summaryrefslogtreecommitdiff
path: root/libm/float/ndtrf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/ndtrf.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/float/ndtrf.c')
-rw-r--r--libm/float/ndtrf.c281
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 );
-
-}