Numerical Estimate of the Inverse Error Function

The function that numerically approximates the inverse error (erf) function $[\large \frac{2}{\sqrt{\pi}}\int_0^xe^{-t^2}dt]^{-1}$ is captured in the C function below - it's a fragment of the source code of my (now free) iOS App SciStatCalc. Note the one iteration of Halley's method to refine the result.

Feel free to use it, if you wish to do so.

double inverf(double x)

{   

    double w,p;

    

    w = -log((1.00-x)*(1.00+x));

    

    if ( w < 5.000000 ) {

        w = w - 2.500000;

        p= 2.81022636e-08;

        p= 3.43273939e-07 + p*w;

        p= -3.5233877e-06 + p*w;

        p = -4.39150654e-06 + p*w;

        p = 0.00021858087 + p*w;

        p = -0.00125372503 + p*w;

        p = -0.00417768164 + p*w;

        p = 0.246640727 + p*w;

        p = 1.50140941 + p*w;

    } else

    {

        w = sqrtf(w) - 3.000000;

        p =  -0.000200214257;

        p = 0.000100950558 + p*w;

        p = 0.00134934322 + p*w;

        p = -0.00367342844 + p*w;

        p = 0.00573950773 + p*w;

        p = -0.0076224613 + p*w;

        p = 0.00943887047 + p*w;

        p = 1.00167406 + p*w;

        p = 2.83297682 + p*w;

    }

    

    double res_ra = p*x; // assign to rational estimate variable

    

    // Halley's method to refine estimate of inverse erf

    double res_hm = 0;

    double fx = (erf(res_ra) - x);

    double df = (2.0/sqrt(M_PI))*exp(-(res_ra * res_ra));

    double d2f = -2 * res_ra * df;

    

    res_hm = res_ra - (2*fx*df)/((2*df*df) - (fx*d2f));

    return res_hm;

}

The above function can be used to estimate the inverse cumulative distribution function (quantile) of the Gaussian distribution (mu is the mean, and sig the standard deviation of the distribution), as below:-

double gauss_icdf(double p, double mu, double sig)
{
    double res;
    res = (inverf((2*p) - 1) *sig*sqrt(2.0)) + mu;
    return res;
}