# include # include # include # include "asa091.h" /******************************************************************************/ double alnorm ( double x, int upper ) /******************************************************************************/ /* Purpose: alnorm() computes the cumulative density of the standard normal distribution. Licensing: This code is distributed under the MIT license. Modified: 03 November 2010 Author: Original FORTRAN77 version by David Hill. C version by John Burkardt. Reference: David Hill, Algorithm AS 66: The Normal Integral, Applied Statistics, Volume 22, Number 3, 1973, pages 424-427. Parameters: Input, double X, is one endpoint of the semi-infinite interval over which the integration takes place. Input, int UPPER, determines whether the upper or lower interval is to be integrated: 1 => integrate from X to + Infinity; 0 => integrate from - Infinity to X. Output, double ALNORM, the integral of the standard normal distribution over the desired interval. */ { double a1 = 5.75885480458; double a2 = 2.62433121679; double a3 = 5.92885724438; double b1 = -29.8213557807; double b2 = 48.6959930692; double c1 = -0.000000038052; double c2 = 0.000398064794; double c3 = -0.151679116635; double c4 = 4.8385912808; double c5 = 0.742380924027; double c6 = 3.99019417011; double con = 1.28; double d1 = 1.00000615302; double d2 = 1.98615381364; double d3 = 5.29330324926; double d4 = -15.1508972451; double d5 = 30.789933034; double ltone = 7.0; double p = 0.398942280444; double q = 0.39990348504; double r = 0.398942280385; int up; double utzero = 18.66; double value; double y; double z; up = upper; z = x; if ( z < 0.0 ) { up = !up; z = - z; } if ( ltone < z && ( ( !up ) || utzero < z ) ) { if ( up ) { value = 0.0; } else { value = 1.0; } return value; } y = 0.5 * z * z; if ( z <= con ) { value = 0.5 - z * ( p - q * y / ( y + a1 + b1 / ( y + a2 + b2 / ( y + a3 )))); } else { value = r * exp ( - y ) / ( z + c1 + d1 / ( z + c2 + d2 / ( z + c3 + d3 / ( z + c4 + d4 / ( z + c5 + d5 / ( z + c6 )))))); } if ( !up ) { value = 1.0 - value; } return value; }