/*	File:	CDF.cpp

			Christer Karlsson


This code implements a function that calculates the 
standard normal CDF (x), using an approximation from 
Abromowitz and Stegun Handbook of Mathematical Functions.

http://www.math.sfu.ca/~cbm/aands/page_932.htm 

   CDF(x) = 1-Z(x)*(b1*t+b2*t^2+b3*t^3+b4*t^4+b5*t^5)+err
            (Where |err| < 7.5x10^-8)    
	    
        t =  1/(1+p*x)

       b1 =  0.319381530
       b2 = -0.356563782
       b3 =  1.781477937
       b4 = -1.821255978
       b5 =  1.330274429
       p  =  0.2316419

Allowed inputs are doubles. The program breaks with ctrl-C   */

#include <iostream>
#include <cmath>
using namespace std;

// Constants and global variables
   const double PI = 3.141592654;

// Functions
// The Gaussian p.d.f with mean = 0 and stddev = 1
double Z(const double x)
{
	return (1/sqrt(2*PI))*exp(-x*x/2.0 );
}

// Start of the  Abromowitz and Stegun approximation function
double CDF(const double x)
{
  const double b1 =  0.319381530;
  const double b2 = -0.356563782;
  const double b3 =  1.781477937;
  const double b4 = -1.821255978;
  const double b5 =  1.330274429;
  const double p  =  0.2316419;

  if(x >= 0.0) {
      double t = 1.0 / (1.0 + p*x);
      return (1.0 - Z(x)*t* 
      (t*(t*(t*(t*b5 + b4) + b3) + b2) + b1));
  } 
  else { 
      double t = 1.0 / ( 1.0 - p * x );
      return ( Z(x)*t* 
      (t*(t*(t*(t*b5 + b4) + b3) + b2) + b1));
  }
}

int main()
{
	double x=0;

	do
	{
		cout << "Enter an X: ";
		cin >> x;
		cout << CDF(x) << '\n';
	} while(1);

	return 0;
}