Doxygen Source Code Documentation
cdf_33.c File Reference
#include "cdflib.h"
Go to the source code of this file.
Functions | |
void | cumnor (double *arg, double *result, double *ccum) |
Function Documentation
|
Definition at line 2 of file cdf_33.c. References a, arg, c, fifdint(), i, p, q, and spmpar(). Referenced by cdfnor(), and dinvnr().
00025 : 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; } /* END */ |