def _interpolate_p(p, r, 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.
    except:
        print(p,r,v)
    y1 = _func(A[(p1, v)], p1, r, v) + 1.
    y2 = _func(A[(p2, v)], p2, r, v) + 1.

    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 > .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.) * (p_t-p1_t)**2. + d1 * (p_t-p1_t) + d0

        # transform back to y
        y = math.exp(y_log) - float(r)/float(v)

    elif p > .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.) * (p-p1)**2. + d1 * (p-p1) + d0

        # transform back to y
        y = math.exp(y_log) - float(r)/float(v)

    else:
        # linear interpolation in q and p
        q0 = math.sqrt(2) * -y0 * \
             scipy.stats.t.isf((1.+p0)/2., (v,1e38)[v>1e38])
        q1 = math.sqrt(2) * -y1 * \
             scipy.stats.t.isf((1.+p1)/2., (v,1e38)[v>1e38])

        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.+p)/2., (v,1e38)[v>1e38]))

    return y
QR factorization.
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/stats/libqsturng/qsturng_.py/_interpolate_p--crossover with--dataset/ETHPy150Open lebedov/scikit-cuda/skcuda/cula.py/culaDeviceZgeqrf