Owen's T function

tfn.py

def tfn(x, fx):
    """
    TFN calculates the T-function of Owen.
    
    Adapted from http://people.sc.fsu.edu/~jburkardt/c_src/asa076/asa076.html
    
    Author:

    Original FORTRAN77 version by JC Young, Christoph Minder.
    C version by John Burkardt.

    Reference:

    MA Porter, DJ Winstanley,
    Remark AS R30:
    A Remark on Algorithm AS76:
    An Integral Useful in Calculating Noncentral T and Bivariate
    Normal Probabilities,
    Applied Statistics,
    Volume 28, Number 1, 1979, page 113.

    JC Young, Christoph Minder,
    Algorithm AS 76: 
    An Algorithm Useful in Calculating Non-Central T and 
    Bivariate Normal Distributions,
    Applied Statistics,
    Volume 23, Number 3, 1974, pages 455-457.    
    """
    r = [
        0.1477621, 
        0.1346334, 
        0.1095432, 
        0.0747257, 
        0.0333357,
    ]
    tp = 0.159155
    tv1 = 1.0e-35
    tv2 = 15.0
    tv3 = 15.0
    tv4 = 1.0e-05
    u = [
        0.0744372, 
        0.2166977, 
        0.3397048,
        0.4325317, 
        0.4869533,
    ]
    # Test for X near zero.
    if np.fabs(x) < tv1:
        return tp * np.arctan(fx)
    # Test for large values of abs(X).
    if tv2 < np.fabs(x):
        return 0.0
    # Test for FX near zero.
    if np.fabs(fx) < tv1:
        return 0.0
    # Test whether abs ( FX ) is so large that it must be truncated.
    xs = - 0.5 * x * x
    x2 = fx
    fxs = fx * fx
    # Computation of truncation point by Newton iteration.
    if tv3 <= np.log1p(fxs) - xs * fxs:
        x1 = 0.5 * fx
        fxs = 0.25 * fxs
        while np.fabs(x2 - x1) >= tv4:
            rt = fxs + 1.0
            x2 = x1 + ( ( xs * fxs + tv3 - np.log(rt) ) 
                        / ( 2.0 * x1 * ( 1.0 / rt - xs ) ) )
            fxs = x2 * x2
            x1 = x2
    # Gaussian quadrature.
    rt = 0.0
    for rr, uu in zip(r, u):
        r1 = 1.0 + fxs * (0.5 + uu)**2
        r2 = 1.0 + fxs * (0.5 - uu)**2
        rt += rr * ( np.exp(xs * r1) / r1 + np.exp(xs * r2) / r2 )
    return rt * x2 * tp
owen_t = np.vectorize(tfn)