diff options
Diffstat (limited to 'libm/float/ndtrf.c')
| -rw-r--r-- | libm/float/ndtrf.c | 281 |
1 files changed, 281 insertions, 0 deletions
diff --git a/libm/float/ndtrf.c b/libm/float/ndtrf.c new file mode 100644 index 000000000..c08d69eca --- /dev/null +++ b/libm/float/ndtrf.c @@ -0,0 +1,281 @@ +/* ndtrf.c + * + * Normal distribution function + * + * + * + * SYNOPSIS: + * + * float x, y, ndtrf(); + * + * y = ndtrf( x ); + * + * + * + * DESCRIPTION: + * + * Returns the area under the Gaussian probability density + * function, integrated from minus infinity to x: + * + * x + * - + * 1 | | 2 + * ndtr(x) = --------- | exp( - t /2 ) dt + * sqrt(2pi) | | + * - + * -inf. + * + * = ( 1 + erf(z) ) / 2 + * = erfc(z) / 2 + * + * where z = x/sqrt(2). Computation is via the functions + * erf and erfc. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -13,0 50000 1.5e-5 2.6e-6 + * + * + * ERROR MESSAGES: + * + * See erfcf(). + * + */ + /* erff.c + * + * Error function + * + * + * + * SYNOPSIS: + * + * float x, y, erff(); + * + * y = erff( x ); + * + * + * + * DESCRIPTION: + * + * The integral is + * + * x + * - + * 2 | | 2 + * erf(x) = -------- | exp( - t ) dt. + * sqrt(pi) | | + * - + * 0 + * + * The magnitude of x is limited to 9.231948545 for DEC + * arithmetic; 1 or -1 is returned outside this range. + * + * For 0 <= |x| < 1, erf(x) = x * P(x**2); otherwise + * erf(x) = 1 - erfc(x). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -9.3,9.3 50000 1.7e-7 2.8e-8 + * + */ + /* erfcf.c + * + * Complementary error function + * + * + * + * SYNOPSIS: + * + * float x, y, erfcf(); + * + * y = erfcf( x ); + * + * + * + * DESCRIPTION: + * + * + * 1 - erf(x) = + * + * inf. + * - + * 2 | | 2 + * erfc(x) = -------- | exp( - t ) dt + * sqrt(pi) | | + * - + * x + * + * + * For small x, erfc(x) = 1 - erf(x); otherwise polynomial + * approximations 1/x P(1/x**2) are computed. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -9.3,9.3 50000 3.9e-6 7.2e-7 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * erfcf underflow x**2 > MAXLOGF 0.0 + * + * + */ + + +/* +Cephes Math Library Release 2.2: June, 1992 +Copyright 1984, 1987, 1988 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +#include <math.h> + +extern float MAXLOGF, SQRTHF; + + +/* erfc(x) = exp(-x^2) P(1/x), 1 < x < 2 */ +static float P[] = { + 2.326819970068386E-002, +-1.387039388740657E-001, + 3.687424674597105E-001, +-5.824733027278666E-001, + 6.210004621745983E-001, +-4.944515323274145E-001, + 3.404879937665872E-001, +-2.741127028184656E-001, + 5.638259427386472E-001 +}; + +/* erfc(x) = exp(-x^2) 1/x P(1/x^2), 2 < x < 14 */ +static float R[] = { +-1.047766399936249E+001, + 1.297719955372516E+001, +-7.495518717768503E+000, + 2.921019019210786E+000, +-1.015265279202700E+000, + 4.218463358204948E-001, +-2.820767439740514E-001, + 5.641895067754075E-001 +}; + +/* erf(x) = x P(x^2), 0 < x < 1 */ +static float T[] = { + 7.853861353153693E-005, +-8.010193625184903E-004, + 5.188327685732524E-003, +-2.685381193529856E-002, + 1.128358514861418E-001, +-3.761262582423300E-001, + 1.128379165726710E+000 +}; + +/*#define UTHRESH 37.519379347*/ + +#define UTHRESH 14.0 + +#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) + +#ifdef ANSIC +float polevlf(float, float *, int); +float expf(float), logf(float), erff(float), erfcf(float); +#else +float polevlf(), expf(), logf(), erff(), erfcf(); +#endif + + + +float ndtrf(float aa) +{ +float x, y, z; + +x = aa; +x *= SQRTHF; +z = fabsf(x); + +if( z < SQRTHF ) + y = 0.5 + 0.5 * erff(x); +else + { + y = 0.5 * erfcf(z); + + if( x > 0 ) + y = 1.0 - y; + } + +return(y); +} + + +float erfcf(float aa) +{ +float a, p,q,x,y,z; + + +a = aa; +x = fabsf(a); + +if( x < 1.0 ) + return( 1.0 - erff(a) ); + +z = -a * a; + +if( z < -MAXLOGF ) + { +under: + mtherr( "erfcf", UNDERFLOW ); + if( a < 0 ) + return( 2.0 ); + else + return( 0.0 ); + } + +z = expf(z); +q = 1.0/x; +y = q * q; +if( x < 2.0 ) + { + p = polevlf( y, P, 8 ); + } +else + { + p = polevlf( y, R, 7 ); + } + +y = z * q * p; + +if( a < 0 ) + y = 2.0 - y; + +if( y == 0.0 ) + goto under; + +return(y); +} + + +float erff(float xx) +{ +float x, y, z; + +x = xx; +if( fabsf(x) > 1.0 ) + return( 1.0 - erfcf(x) ); + +z = x * x; +y = x * polevlf( z, T, 6 ); +return( y ); + +} |