Code for generating a Padé approximant for the erfinv function.

erfinv_approx.py

"""
erfinv_approx.py

Create a Pade approximant for the erfinv function.
"""

import mpmath
import numpy as np


def create_pade_approx():
    # dps = 400 is way overkill, but there is reason for it.  If I set dps
    # to M, some of the terms in ts and in p and q have magnitude on the
    # order of 10**-M or smaller.  Presumably these should be zero.  By
    # making dps = 400, these terms are sure to converted to 0 when the
    # coefficients are converted to float64.
    with mpmath.workdps(400):

        def erfinvy_over_y(y):
            return mpmath.erfinv(y)/y if y != 0 else mpmath.sqrt(mpmath.pi)/2

        # Compute the Taylor polynomial of erfinv(y)/y at 0.
        ts = mpmath.taylor(erfinvy_over_y, 0, 20)

        # Convert the Taylor series to the Pade rational approximation.
        p, q = mpmath.pade(ts, 10, 10)

        # By construction, every other coefficient in p and q is 0.
        # Drop the zeros to create NumPy polynomials whose argument is
        # expected to be y**2.
        P2 = np.polynomial.Polynomial([float(c) for c in p[::2]])
        Q2 = np.polynomial.Polynomial([float(c) for c in q[::2]])
        return P2, Q2


def check_rel_error(P2, Q2):
    x = np.concatenate((np.geomspace(1e-80, 1e-7, 250),
                        np.linspace(1.1e-7, 0.25, 1500)))
    x = np.concatenate((-x[::-1], x))
    max_relerr = 0
    for k, v in enumerate(x):
        with mpmath.workdps(50):
            exact = float(mpmath.erfinv(v))
        v2 = v**2
        delta = exact - v*(P2(v2)/Q2(v2))
        relerr = abs(delta/exact)
        if relerr > max_relerr:
            max_relerr = relerr
    return max_relerr


def generate_c_code(P2, Q2):
    print('double polevl(double x, const double coef[], int N);')
    print()
    print('/*')
    print(' *  Compute erfinv(y) using a Pade approximation.')
    print(' *  Tests suggest that on the interval [-0.25, 0.25],')
    print(' *  the relative error is roughly 3e-16 or less.')
    print(' */')
    print('double erfinv_pade_approx(double y)')
    print('{')
    print('    // Coefficients are in reverse order; the constant term')
    print('    // is the last term in the array.')
    s1 = '    double P2[] = {'
    print(s1, end='')
    print((',\n' + ' '*len(s1)).join([repr(v) for v in P2.coef[::-1]]), end='')
    print('};')
    s2 = '    double Q2[] = {'
    print(s2, end='')
    print((',\n' + ' '*len(s2)).join([repr(v) for v in Q2.coef[::-1]]), end='')
    print('};')
    print('    double y2 = y*y;')
    print(f'    double num = polevl(y2, P2, {len(P2)-1});')
    print(f'    double den = polevl(y2, Q2, {len(Q2)-1});')
    print('    return y * (num / den);')
    print('}')