normsinv
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

//@version=5
// @description Returns the inverse of the standard normal cumulative distribution. 
//      The distribution has a mean of zero and a standard deviation of one; i.e., 
//      normsinv seeks that value z such that a normal distribtuio of mean of zero 
//      and standard deviation one is equal to the input probability.
// Reference:
//		https://github.com/jeremybarnes/cephes/blob/master/cprob/ndtri.c
library("normsinv")
//
//     Inverse of Normal distribution function
//
//
//
// SYNOPSIS:
//
// double x, y, ndtri();
//
// x = ndtri( y );
//
//
//
// DESCRIPTION:
//
// Returns the argument, x, for which the area under the
// Gaussian probability density function (integrated from
// minus infinity to x) is equal to y.
//
//
// For small arguments 0 < y < exp(-2), the program computes
// z = sqrt( -2.0 // log(y) );  then the approximation is
// x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
// There are two rational functions P/Q, one for 0 < y < exp(-32)
// and the other for y up to exp(-2).  For larger arguments,
// w = y - 0.5, and  x/sqrt(2pi) = w + w////3 R(w////2)/S(w////2)).
//
//
// ACCURACY:
//
//                      Relative error:
// arithmetic   domain        # trials      peak         rms
//    IEEE     0.125, 1        20000       7.2e-16     1.3e-16
//    IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17
//
//
// ERROR MESSAGES:
//
//   message         condition    value returned
// ndtri domain       x <= 0        -1E10
// ndtri domain       x >= 1         1E10
//
 
// Evaluates polynomial of degree N
polevl(float x, float[] coef, int n)=>
    float ans = array.get(coef, 0)
    for i = 1 to n
        ans := ans * x + array.get(coef, i)
    ans

// Evaluates polynomial of degree N with assumption that coef[N] = 1.0
p1evl(float x, float[] coef, int n)=>
    float ans = x + array.get(coef, 0)
    for i = 1 to n - 1
        ans := ans * x + array.get(coef, i)
    ans

// @function Returns the inverse of the standard normal cumulative distribution. The distribution has a mean of zero and a standard deviation of one.
// @param y0 float, probability corresponding to the normal distribution.
// @returns float, z-score
export normsinv(float y0)=>

    // Cephes Math Library Release 2.1:  January, 1989
    // Copyright 1984, 1987, 1989 by Stephen L. Moshier
    // Direct inquiries to 30 Frost Street, Cambridge, MA 02140

    // sqrt(2pi) 
    float s2pi = 2.50662827463100050242E0

    // approximation for 0 <= |y - 0.5| <= 3/8 
    var array<float> P0 = array.new_float(5, 0)
    array.set(P0, 0, -5.99633501014107895267E1)
    array.set(P0, 1, 9.80010754185999661536E1)
    array.set(P0, 2, -5.66762857469070293439E1)
    array.set(P0, 3, 1.39312609387279679503E1)
    array.set(P0, 4, -1.23916583867381258016E0)

    var array<float> Q0 = array.new_float(8, 0)
    array.set(Q0, 0, 1.95448858338141759834E0)
    array.set(Q0, 1, 4.67627912898881538453E0)
    array.set(Q0, 2, 8.63602421390890590575E1)
    array.set(Q0, 3, -2.25462687854119370527E2)
    array.set(Q0, 4, 2.00260212380060660359E2)
    array.set(Q0, 5, -8.20372256168333339912E1)
    array.set(Q0, 6, 1.59056225126211695515E1)
    array.set(Q0, 7, -1.18331621121330003142E0)

    // Approximation for interval z = sqrt(-2 log y ) between 2 and 8
    // i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
    var array<float> P1 = array.new_float(9, 0)
    array.set(P1, 0, 4.05544892305962419923E0)
    array.set(P1, 1, 3.15251094599893866154E1)
    array.set(P1, 2, 5.71628192246421288162E1)
    array.set(P1, 3, 4.40805073893200834700E1)
    array.set(P1, 4, 1.46849561928858024014E1)
    array.set(P1, 5, 2.18663306850790267539E0)
    array.set(P1, 6, -1.40256079171354495875E-1)
    array.set(P1, 7, -3.50424626827848203418E-2)
    array.set(P1, 8, -8.57456785154685413611E-4)

    var array<float> Q1 = array.new_float(8, 0)
    array.set(Q1, 0, 1.57799883256466749731E1)
    array.set(Q1, 1, 4.53907635128879210584E1)
    array.set(Q1, 2, 4.13172038254672030440E1)
    array.set(Q1, 3, 1.50425385692907503408E1)
    array.set(Q1, 4, 2.50464946208309415979E0)
    array.set(Q1, 5, -1.42182922854787788574E-1)
    array.set(Q1, 6, -3.80806407691578277194E-2)
    array.set(Q1, 7, -9.33259480895457427372E-4)

    // Approximation for interval z = sqrt(-2 log y ) between 8 and 64
    // i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    var array<float> P2 = array.new_float(9, 0)
    array.set(P2, 0, 3.23774891776946035970E0)
    array.set(P2, 1, 6.91522889068984211695E0)
    array.set(P2, 2, 3.93881025292474443415E0)
    array.set(P2, 3, 1.33303460815807542389E0)
    array.set(P2, 4, 2.01485389549179081538E-1)
    array.set(P2, 5, 1.23716634817820021358E-2)
    array.set(P2, 6, 3.01581553508235416007E-4)
    array.set(P2, 7, 2.65806974686737550832E-6)
    array.set(P2, 8, 6.23974539184983293730E-9)

    var array<float> Q2 = array.new_float(8, 0)
    array.set(Q2, 0, 6.02427039364742014255E0)
    array.set(Q2, 1, 3.67983563856160859403E0)
    array.set(Q2, 2, 1.37702099489081330271E0)
    array.set(Q2, 3, 2.16236993594496635890E-1)
    array.set(Q2, 4, 1.34204006088543189037E-2)
    array.set(Q2, 5, 3.28014464682127739104E-4)
    array.set(Q2, 6, 2.89247864745380683936E-6)
    array.set(Q2, 7, 6.79019408009981274425E-9)

    float x = 0
    float y = 0
    float z = 0
    float y2 = 0
    float x0 = 0
    float x1 = 0
    int code = 0
    if not(y0 <= 0. or y0 >= 1.)
        code := 1
        y := y0
        // 0.135... = exp(-2) 
        if y > (1.0 - 0.13533528323661269189)
            y := 1.0 - y
            code := 0
        if y > 0.13533528323661269189
            y := y - 0.5
            y2 := y * y
            x := y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8))
            x := x * s2pi
        else 
            x := math.sqrt(-2.0 *  math.log(y))
            x0 := x - math.log(x) / x
            z := 1.0 / x
            // y > exp(-32) = 1.2664165549e-14 
            if x < 8.
                x1 := z * polevl(z, P1, 8) / p1evl(z, Q1, 8)
            else
                x1 := z * polevl(z, P2, 8) / p1evl(z, Q2, 8)
            x := x0 - x1
            if (code != 0)
                x := -x
    else 
        if y0 <= 0.
            x := -1E10
        if y0 >= 1
            x := 1E10

    x
    
// example usage
out = normsinv(0.995)
plot(out)