Stand-alone C# code for Phi inverse

The following code first appeared as A literate program to compute the inverse of the normal CDF. See that page for a detailed explanation of the algorithm.

using System;

class Program
{
    static double RationalApproximation(double t)
    {
        // Abramowitz and Stegun formula 26.2.23.
        // The absolute value of the error should be less than 4.5 e-4.
        double[] c = {2.515517, 0.802853, 0.010328};
        double[] d = {1.432788, 0.189269, 0.001308};
        return t - ((c[2]*t + c[1])*t + c[0]) / 
                    (((d[2]*t + d[1])*t + d[0])*t + 1.0);
    }

    static double NormalCDFInverse(double p)
    {
        if (p <= 0.0 || p >= 1.0)
        {
            string msg = String.Format("Invalid input argument: {0}.", p);
            throw new ArgumentOutOfRangeException(msg);
        }

        // See article above for explanation of this section.
        if (p < 0.5)
        {
            // F^-1(p) = - G^-1(p)
            return -RationalApproximation( Math.Sqrt(-2.0*Math.Log(p)) );
        }
        else
        {
            // F^-1(p) = G^-1(1-p)
            return RationalApproximation( Math.Sqrt(-2.0*Math.Log(1.0 - p)) );
        }
    }

    static void demo()
    {
        Console.WriteLine("\nShow NormalCDFInverse is accurate at");
        Console.WriteLine("0.05, 0.15, 0.25, ..., 0.95 ");
        Console.WriteLine("and at a few extreme values.\n");

        double[] p =
        {
            0.0000001,
            0.00001,
            0.001,
            0.05,
            0.15,
            0.25,
            0.35,
            0.45,
            0.55,
            0.65,
            0.75,
            0.85,
            0.95,
            0.999,
            0.99999,
            0.9999999
        };

        // Exact values computed by Mathematica.
        double[] exact =
        {
            -5.199337582187471,
            -4.264890793922602,
            -3.090232306167813,
            -1.6448536269514729,
            -1.0364333894937896,
            -0.6744897501960817,
            -0.38532046640756773,
            -0.12566134685507402,
             0.12566134685507402,
             0.38532046640756773,
             0.6744897501960817,
             1.0364333894937896,
             1.6448536269514729,
             3.090232306167813,
             4.264890793922602,
             5.199337582187471
        };

        double maxerror = 0.0;
    
        Console.WriteLine("p, exact CDF inverse, computed CDF inverse, diff\n");

        for (int i = 0; i < exact.Length; ++i)
        {
            double computed = NormalCDFInverse(p[i]);
            double error = exact[i] - computed;
            Console.WriteLine("{0}, {1}, {2}, {3}", p[i], exact[i], computed, error);
            if (Math.Abs(error) > maxerror)
                maxerror = Math.Abs(error);
        }

        Console.WriteLine("\nMaximum error: {0}\n", maxerror);
    }


    static void Main(string[] args)
    {
        demo();
    }
}

The code is based on an algorithm given in Handbook of Mathematical Functions by Abramowitz and Stegun.

This code is in the public domain. Do whatever you want to with it, no strings attached.

Other versions: C++, C#

More stand-alone numerical code