3 * Normal distribution function
17 * Returns the area under the Gaussian probability density
18 * function, integrated from minus infinity to x:
23 * ndtr(x) = --------- | exp( - t /2 ) dt
28 * = ( 1 + erf(z) ) / 2
31 * where z = x/sqrt(2). Computation is via the functions
38 * arithmetic domain # trials peak rms
39 * IEEE -13,0 50000 1.5e-5 2.6e-6
68 * erf(x) = -------- | exp( - t ) dt.
73 * The magnitude of x is limited to 9.231948545 for DEC
74 * arithmetic; 1 or -1 is returned outside this range.
76 * For 0 <= |x| < 1, erf(x) = x * P(x**2); otherwise
77 * erf(x) = 1 - erfc(x).
84 * arithmetic domain # trials peak rms
85 * IEEE -9.3,9.3 50000 1.7e-7 2.8e-8
90 * Complementary error function
96 * float x, y, erfcf();
110 * erfc(x) = -------- | exp( - t ) dt
116 * For small x, erfc(x) = 1 - erf(x); otherwise polynomial
117 * approximations 1/x P(1/x**2) are computed.
124 * arithmetic domain # trials peak rms
125 * IEEE -9.3,9.3 50000 3.9e-6 7.2e-7
130 * message condition value returned
131 * erfcf underflow x**2 > MAXLOGF 0.0
138 Cephes Math Library Release 2.2: June, 1992
139 Copyright 1984, 1987, 1988 by Stephen L. Moshier
140 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
146 extern float MAXLOGF, SQRTHF;
149 /* erfc(x) = exp(-x^2) P(1/x), 1 < x < 2 */
151 2.326819970068386E-002,
152 -1.387039388740657E-001,
153 3.687424674597105E-001,
154 -5.824733027278666E-001,
155 6.210004621745983E-001,
156 -4.944515323274145E-001,
157 3.404879937665872E-001,
158 -2.741127028184656E-001,
159 5.638259427386472E-001
162 /* erfc(x) = exp(-x^2) 1/x P(1/x^2), 2 < x < 14 */
164 -1.047766399936249E+001,
165 1.297719955372516E+001,
166 -7.495518717768503E+000,
167 2.921019019210786E+000,
168 -1.015265279202700E+000,
169 4.218463358204948E-001,
170 -2.820767439740514E-001,
171 5.641895067754075E-001
174 /* erf(x) = x P(x^2), 0 < x < 1 */
176 7.853861353153693E-005,
177 -8.010193625184903E-004,
178 5.188327685732524E-003,
179 -2.685381193529856E-002,
180 1.128358514861418E-001,
181 -3.761262582423300E-001,
182 1.128379165726710E+000
185 /*#define UTHRESH 37.519379347*/
189 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
192 float polevlf(float, float *, int);
193 float expf(float), logf(float), erff(float), erfcf(float);
195 float polevlf(), expf(), logf(), erff(), erfcf();
200 float ndtrf(float aa)
209 y = 0.5 + 0.5 * erff(x);
222 float erfcf(float aa)
231 return( 1.0 - erff(a) );
238 mtherr( "erfcf", UNDERFLOW );
250 p = polevlf( y, P, 8 );
254 p = polevlf( y, R, 7 );
275 return( 1.0 - erfcf(x) );
278 y = x * polevlf( z, T, 6 );