/*
  Lower tail quantile for standard normal distribution function.
  This function returns an approximation of the inverse cumulative
  standard normal distribution function.  I.e., given P, it returns
  an approximation to the X satisfying P = Pr{Z <= X} where Z is a
  random variable from the standard normal distribution.

  The algorithm uses a minimax approximation by rational functions
  and the result has a relative error whose absolute value is less
  than 1.15e-9.

  The shell syntax is (or should be, anyway) POSIX bc.

  Author:      Peter John Acklam
  Time-stamp:  2005-03-10 14:13:52 +01:00
  E-mail:      pjacklam@online.no
  WWW URL:     http://home.online.no/~pjacklam
*/

define z (p) {
    q = p - 0.5;

    a = q;
    if (a < 0) {
        a = -a;
    }

    /*
       Rational approximation for central region
    */

    if (a <= .47575) {

        r = q*q;
        z =  ((((( -39.69683028665376  * r + 220.9460984245205) * r \
                  - 275.9285104469687) * r + 138.3577518672690) * r \
                  - 30.66479806614716) * r + 2.506628277459239) * q \
           / ((((( -54.47609879822406  * r + 161.5858368580409) * r \
                  - 155.6989798598866) * r + 66.80131188771972) * r \
                  - 13.28068155288572) * r + 1);

        return(z);
    }

    /* Rational approximation for tails */

    if (a > .47575) {

        /* If in upper tail, map to lower tail */

        if (q > 0) {
            p = 1 - p;
        }

        r  = sqrt(-2 * l(p));
        z =  (((((-0.007784894002430293 * r - 0.3223964580411365)  * r \
                  -2.400758277161838)   * r - 2.549732539343734)   * r \
                  +4.374664141464968)   * r + 2.938163982698783)       \
           / ((((  0.007784695709041462 * r + 0.3224671290700398)  * r \
                 + 2.445134137142996)   * r + 3.754408661907416)   * r \
                 + 1);

        /* If in upper tail, swap sign */

        if (q > 0) {
            z = -z;
        }

        return(z);
    }
}

/*
  Emacs configuration:
  ====================
  Local Variables:
  mode: c
  eval: (c-set-style "stroustrup")
  End:
*/