code: def _interpolate_p(p, r, v):
      """
      interpolates p based on the values in the A table for the
      scalar value of r and the scalar value of v
      """
  
      # interpolate p (v should be in table)
      # if .5 < p < .75 use linear interpolation in q
      # if p > .75 use quadratic interpolation in log(y + r/v)
      # by -1. / (1. + 1.5 * _phi((1. + p)/2.))
  
      # find the 3 closest v values
      p0, p1, p2 = _select_ps(p)
      try:
          y0 = _func(A[(p0, v)], p0, r, v) + 1.0
      except:
          print(p, r, v)
          raise
      y1 = _func(A[(p1, v)], p1, r, v) + 1.0
      y2 = _func(A[(p2, v)], p2, r, v) + 1.0
  
      y_log0 = math.log(y0 + float(r) / float(v))
      y_log1 = math.log(y1 + float(r) / float(v))
      y_log2 = math.log(y2 + float(r) / float(v))
  
      # If p < .85 apply only the ordinate transformation
      # if p > .85 apply the ordinate and the abcissa transformation
      # In both cases apply quadratic interpolation
      if p > 0.85:
          p_t = _ptransform(p)
          p0_t = _ptransform(p0)
          p1_t = _ptransform(p1)
          p2_t = _ptransform(p2)
  
          # calculate derivatives for quadratic interpolation
          d2 = (
              2
              * (
                  (y_log2 - y_log1) / (p2_t - p1_t)
                  - (y_log1 - y_log0) / (p1_t - p0_t)
              )
              / (p2_t - p0_t)
          )
          if (p2 + p0) >= (p1 + p1):
              d1 = (y_log2 - y_log1) / (p2_t - p1_t) - 0.5 * d2 * (
                  p2_t - p1_t
              )
          else:
              d1 = (y_log1 - y_log0) / (p1_t - p0_t) + 0.5 * d2 * (
                  p1_t - p0_t
              )
          d0 = y_log1
  
          # interpolate value
          y_log = (
              (d2 / 2.0) * (p_t - p1_t) ** 2.0 + d1 * (p_t - p1_t) + d0
          )
  
          # transform back to y
          y = math.exp(y_log) - float(r) / float(v)
  
      elif p > 0.5:
          # calculate derivatives for quadratic interpolation
          d2 = (
              2
              * (
                  (y_log2 - y_log1) / (p2 - p1)
                  - (y_log1 - y_log0) / (p1 - p0)
              )
              / (p2 - p0)
          )
          if (p2 + p0) >= (p1 + p1):
              d1 = (y_log2 - y_log1) / (p2 - p1) - 0.5 * d2 * (p2 - p1)
          else:
              d1 = (y_log1 - y_log0) / (p1 - p0) + 0.5 * d2 * (p1 - p0)
          d0 = y_log1
  
          # interpolate values
          y_log = (d2 / 2.0) * (p - p1) ** 2.0 + d1 * (p - p1) + d0
  
          # transform back to y
          y = math.exp(y_log) - float(r) / float(v)
  
      else:
          # linear interpolation in q and p
          v = min(v, 1e38)
          q0 = (
              math.sqrt(2)
              * -y0
              * scipy.stats.t.isf((1.0 + p0) / 2.0, v)
          )
          q1 = (
              math.sqrt(2)
              * -y1
              * scipy.stats.t.isf((1.0 + p1) / 2.0, v)
          )
  
          d1 = (q1 - q0) / (p1 - p0)
          d0 = q0
  
          # interpolate values
          q = d1 * (p - p0) + d0
  
          # transform back to y
          y = -q / (
              math.sqrt(2) * scipy.stats.t.isf((1.0 + p) / 2.0, v)
          )
  
      return y
  
file_path: statsmodels/statsmodels/stats/libqsturng/qsturng_.py
imports: from statsmodels.compat.python import lrange
  import math
  import scipy.stats
  import numpy
  from scipy.optimize import fminbound
label_desc: Use image node_id 6 for calling a global function with example usage: _interpolate_p(p, r, v) and returns: y