diff options
Diffstat (limited to 'libm/float/ndtrf.c')
-rw-r--r-- | libm/float/ndtrf.c | 281 |
1 files changed, 0 insertions, 281 deletions
diff --git a/libm/float/ndtrf.c b/libm/float/ndtrf.c deleted file mode 100644 index c08d69eca..000000000 --- a/libm/float/ndtrf.c +++ /dev/null @@ -1,281 +0,0 @@ -/* 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 ); - -} |