00001 #include "cdflib.h"
00002 void cumnor(double *arg,double *result,double *ccum)
00003 /*
00004 **********************************************************************
00005  
00006      void cumnor(double *arg,double *result,double *ccum)
00007  
00008  
00009                               Function
00010  
00011  
00012      Computes the cumulative  of    the  normal   distribution,   i.e.,
00013      the integral from -infinity to x of
00014           (1/sqrt(2*pi)) exp(-u*u/2) du
00015  
00016      X --> Upper limit of integration.
00017                                         X is DOUBLE PRECISION
00018  
00019      RESULT <-- Cumulative normal distribution.
00020                                         RESULT is DOUBLE PRECISION
00021  
00022      CCUM <-- Compliment of Cumulative normal distribution.
00023                                         CCUM is DOUBLE PRECISION
00024  
00025      Renaming of function ANORM from:
00026 
00027      Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
00028      Package of Special Function Routines and Test Drivers"
00029      acm Transactions on Mathematical Software. 19, 22-32.
00030 
00031      with slight modifications to return ccum and to deal with
00032      machine constants.
00033  
00034 **********************************************************************
00035   Original Comments:
00036 ------------------------------------------------------------------
00037  
00038  This function evaluates the normal distribution function:
00039  
00040                               / x
00041                      1       |       -t*t/2
00042           P(x) = ----------- |      e       dt
00043                  sqrt(2 pi)  |
00044                              /-oo
00045  
00046    The main computation evaluates near-minimax approximations
00047    derived from those in "Rational Chebyshev approximations for
00048    the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
00049    This transportable program uses rational functions that
00050    theoretically approximate the normal distribution function to
00051    at least 18 significant decimal digits.  The accuracy achieved
00052    depends on the arithmetic system, the compiler, the intrinsic
00053    functions, and proper selection of the machine-dependent
00054    constants.
00055  
00056 *******************************************************************
00057 *******************************************************************
00058  
00059  Explanation of machine-dependent constants.
00060  
00061    MIN   = smallest machine representable number.
00062  
00063    EPS   = argument below which anorm(x) may be represented by
00064            0.5  and above which  x*x  will not underflow.
00065            A conservative value is the largest machine number X
00066            such that   1.0 + X = 1.0   to machine precision.
00067 *******************************************************************
00068 *******************************************************************
00069  
00070  Error returns
00071  
00072   The program returns  ANORM = 0     for  ARG .LE. XLOW.
00073  
00074  
00075  Intrinsic functions required are:
00076  
00077      ABS, AINT, EXP
00078  
00079  
00080   Author: W. J. Cody
00081           Mathematics and Computer Science Division
00082           Argonne National Laboratory
00083           Argonne, IL 60439
00084  
00085   Latest modification: March 15, 1992
00086  
00087 ------------------------------------------------------------------
00088 */
00089 {
00090 static double a[5] = {
00091     2.2352520354606839287e00,1.6102823106855587881e02,1.0676894854603709582e03,
00092     1.8154981253343561249e04,6.5682337918207449113e-2
00093 };
00094 static double b[4] = {
00095     4.7202581904688241870e01,9.7609855173777669322e02,1.0260932208618978205e04,
00096     4.5507789335026729956e04
00097 };
00098 static double c[9] = {
00099     3.9894151208813466764e-1,8.8831497943883759412e00,9.3506656132177855979e01,
00100     5.9727027639480026226e02,2.4945375852903726711e03,6.8481904505362823326e03,
00101     1.1602651437647350124e04,9.8427148383839780218e03,1.0765576773720192317e-8
00102 };
00103 static double d[8] = {
00104     2.2266688044328115691e01,2.3538790178262499861e02,1.5193775994075548050e03,
00105     6.4855582982667607550e03,1.8615571640885098091e04,3.4900952721145977266e04,
00106     3.8912003286093271411e04,1.9685429676859990727e04
00107 };
00108 static double half = 0.5e0;
00109 static double p[6] = {
00110     2.1589853405795699e-1,1.274011611602473639e-1,2.2235277870649807e-2,
00111     1.421619193227893466e-3,2.9112874951168792e-5,2.307344176494017303e-2
00112 };
00113 static double one = 1.0e0;
00114 static double q[5] = {
00115     1.28426009614491121e00,4.68238212480865118e-1,6.59881378689285515e-2,
00116     3.78239633202758244e-3,7.29751555083966205e-5
00117 };
00118 static double sixten = 1.60e0;
00119 static double sqrpi = 3.9894228040143267794e-1;
00120 static double thrsh = 0.66291e0;
00121 static double root32 = 5.656854248e0;
00122 static double zero = 0.0e0;
00123 static int K1 = 1;
00124 static int K2 = 2;
00125 static int i;
00126 static double del,eps,temp,x,xden,xnum,y,xsq,min;
00127 /*
00128 ------------------------------------------------------------------
00129   Machine dependent constants
00130 ------------------------------------------------------------------
00131 */
00132     eps = spmpar(&K1)*0.5e0;
00133     min = spmpar(&K2);
00134     x = *arg;
00135     y = fabs(x);
00136     if(y <= thrsh) {
00137 /*
00138 ------------------------------------------------------------------
00139   Evaluate  anorm  for  |X| <= 0.66291
00140 ------------------------------------------------------------------
00141 */
00142         xsq = zero;
00143         if(y > eps) xsq = x*x;
00144         xnum = a[4]*xsq;
00145         xden = xsq;
00146         for(i=0; i<3; i++) {
00147             xnum = (xnum+a[i])*xsq;
00148             xden = (xden+b[i])*xsq;
00149         }
00150         *result = x*(xnum+a[3])/(xden+b[3]);
00151         temp = *result;
00152         *result = half+temp;
00153         *ccum = half-temp;
00154     }
00155 /*
00156 ------------------------------------------------------------------
00157   Evaluate  anorm  for 0.66291 <= |X| <= sqrt(32)
00158 ------------------------------------------------------------------
00159 */
00160     else if(y <= root32) {
00161         xnum = c[8]*y;
00162         xden = y;
00163         for(i=0; i<7; i++) {
00164             xnum = (xnum+c[i])*y;
00165             xden = (xden+d[i])*y;
00166         }
00167         *result = (xnum+c[7])/(xden+d[7]);
00168         xsq = fifdint(y*sixten)/sixten;
00169         del = (y-xsq)*(y+xsq);
00170         *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
00171         *ccum = one-*result;
00172         if(x > zero) {
00173             temp = *result;
00174             *result = *ccum;
00175             *ccum = temp;
00176         }
00177     }
00178 /*
00179 ------------------------------------------------------------------
00180   Evaluate  anorm  for |X| > sqrt(32)
00181 ------------------------------------------------------------------
00182 */
00183     else  {
00184         *result = zero;
00185         xsq = one/(x*x);
00186         xnum = p[5]*xsq;
00187         xden = xsq;
00188         for(i=0; i<4; i++) {
00189             xnum = (xnum+p[i])*xsq;
00190             xden = (xden+q[i])*xsq;
00191         }
00192         *result = xsq*(xnum+p[4])/(xden+q[4]);
00193         *result = (sqrpi-*result)/y;
00194         xsq = fifdint(x*sixten)/sixten;
00195         del = (x-xsq)*(x+xsq);
00196         *result = exp(-(xsq*xsq*half))*exp(-(del*half))**result;
00197         *ccum = one-*result;
00198         if(x > zero) {
00199             temp = *result;
00200             *result = *ccum;
00201             *ccum = temp;
00202         }
00203     }
00204     if(*result < min) *result = 0.0e0;
00205 /*
00206 ------------------------------------------------------------------
00207   Fix up for negative argument, erf, etc.
00208 ------------------------------------------------------------------
00209 ----------Last card of ANORM ----------
00210 */
00211     if(*ccum < min) *ccum = 0.0e0;
00212 } /* END */