/***************************************************
* Program to demonstrate inverse normal subroutine *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*              C++ version by J.-P. Moreau, Paris. *
*                      (www.jpmoreau.fr)           *
* ------------------------------------------------ *
* SAMPLE RUN:                                      *
*                                                  *
*  P(Z>X)     X                                    *
* ----------------                                 *
*  0.50    0.0000                                  *
*  0.48    0.0500                                  *
*  0.46    0.1002                                  *
*  0.44    0.1507                                  *
*  0.42    0.2015                                  *
*  0.40    0.2529                                  *
*  0.38    0.3050                                  *
*  0.36    0.3580                                  *
*  0.34    0.4120                                  *
*  0.32    0.4673                                  *
*  0.30    0.5240                                  *
*  0.28    0.5825                                  *
*  0.26    0.6430                                  *
*  0.24    0.7060                                  *
*  0.22    0.7719                                  *
*  0.20    0.8414                                  *
*  0.18    0.9152                                  *
*  0.16    0.9944                                  *
*  0.14    1.0804                                  *
*  0.12    1.1751                                  *
*  0.10    1.2817                                  *
*  0.08    1.4053                                  *
*  0.06    1.5551                                  *
*  0.04    1.7511                                  *
*  0.02    2.0542                                  *
* -0.00    INF.                                    *
*                                                  *
***************************************************/
#include <stdio.h>
#include <math.h>

int    i;
double x,y;


/**********************************************
*   Inverse normal distribution subroutine    *
* ------------------------------------------- *
* This program calculates an approximation to *
* the integral of the normal distribution     *
* function from x to infinity (the tail).     *
* A rational polynomial is used. The input is *
* in y, with the result returned in x. The    *
* accuracy is better then 0.0005 in the range *
* 0 < y < 0.5.                                *
* ------------------------------------------- *
* Reference: Abramowitz and Stegun.           *
**********************************************/
void Inverse_Normal()  {
  double c0,c1,c2,d1,d2,d3,z;
  // Define coefficients
  c0 = 2.515517;
  c1 = 0.802853;
  c2 = 0.010328;
  d1 = 1.432788;
  d2 = 0.189269;
  d3 = 0.001308;
  if (y <= 0)  {
    x = 1e13;
    return;
  }
  z = sqrt(-log(y * y));
  x = 1.0 + d1 * z + d2 * z * z + d3 * z * z * z;
  x = (c0 + c1 * z + c2 * z * z) / x;
  x = z - x;
}

void main()  {
  printf("\n  P(Z>X)     X   \n");
  printf(" ----------------\n");
  i = 1; y = 0.5;
  while (y >= -0.01)  {
    Inverse_Normal();
    if (x < 0.000001) x = 0;
    if (x < 1e10) 
      printf(" %5.2f    %6.4f\n",y,x);
    else
      printf(" %5.2f    INF.\n",y);
    if (i==20) {
      getchar(); i = 1;
    }
    i++;
    y -= 0.02;
  }
  printf("\n");
}


// End of file invnorm.cpp