diff options
Diffstat (limited to 'libm/double/ndtr.c')
-rw-r--r-- | libm/double/ndtr.c | 481 |
1 files changed, 0 insertions, 481 deletions
diff --git a/libm/double/ndtr.c b/libm/double/ndtr.c deleted file mode 100644 index 75d59ab54..000000000 --- a/libm/double/ndtr.c +++ /dev/null @@ -1,481 +0,0 @@ -/* ndtr.c - * - * Normal distribution function - * - * - * - * SYNOPSIS: - * - * double x, y, ndtr(); - * - * y = ndtr( 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 - * DEC -13,0 8000 2.1e-15 4.8e-16 - * IEEE -13,0 30000 3.4e-14 6.7e-15 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * erfc underflow x > 37.519379347 0.0 - * - */ -/* erf.c - * - * Error function - * - * - * - * SYNOPSIS: - * - * double x, y, erf(); - * - * y = erf( 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 * P4(x**2)/Q5(x**2); otherwise - * erf(x) = 1 - erfc(x). - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * DEC 0,1 14000 4.7e-17 1.5e-17 - * IEEE 0,1 30000 3.7e-16 1.0e-16 - * - */ -/* erfc.c - * - * Complementary error function - * - * - * - * SYNOPSIS: - * - * double x, y, erfc(); - * - * y = erfc( 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 rational - * approximations are computed. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * DEC 0, 9.2319 12000 5.1e-16 1.2e-16 - * IEEE 0,26.6417 30000 5.7e-14 1.5e-14 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * erfc underflow x > 9.231948545 (DEC) 0.0 - * - * - */ - - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier -*/ - - -#include <math.h> - -extern double SQRTH; -extern double MAXLOG; - - -#ifdef UNK -static double P[] = { - 2.46196981473530512524E-10, - 5.64189564831068821977E-1, - 7.46321056442269912687E0, - 4.86371970985681366614E1, - 1.96520832956077098242E2, - 5.26445194995477358631E2, - 9.34528527171957607540E2, - 1.02755188689515710272E3, - 5.57535335369399327526E2 -}; -static double Q[] = { -/* 1.00000000000000000000E0,*/ - 1.32281951154744992508E1, - 8.67072140885989742329E1, - 3.54937778887819891062E2, - 9.75708501743205489753E2, - 1.82390916687909736289E3, - 2.24633760818710981792E3, - 1.65666309194161350182E3, - 5.57535340817727675546E2 -}; -static double R[] = { - 5.64189583547755073984E-1, - 1.27536670759978104416E0, - 5.01905042251180477414E0, - 6.16021097993053585195E0, - 7.40974269950448939160E0, - 2.97886665372100240670E0 -}; -static double S[] = { -/* 1.00000000000000000000E0,*/ - 2.26052863220117276590E0, - 9.39603524938001434673E0, - 1.20489539808096656605E1, - 1.70814450747565897222E1, - 9.60896809063285878198E0, - 3.36907645100081516050E0 -}; -static double T[] = { - 9.60497373987051638749E0, - 9.00260197203842689217E1, - 2.23200534594684319226E3, - 7.00332514112805075473E3, - 5.55923013010394962768E4 -}; -static double U[] = { -/* 1.00000000000000000000E0,*/ - 3.35617141647503099647E1, - 5.21357949780152679795E2, - 4.59432382970980127987E3, - 2.26290000613890934246E4, - 4.92673942608635921086E4 -}; - -#define UTHRESH 37.519379347 -#endif - -#ifdef DEC -static unsigned short P[] = { -0030207,0054445,0011173,0021706, -0040020,0067272,0030661,0122075, -0040756,0151236,0173053,0067042, -0041502,0106175,0062555,0151457, -0042104,0102525,0047401,0003667, -0042403,0116176,0011446,0075303, -0042551,0120723,0061641,0123275, -0042600,0070651,0007264,0134516, -0042413,0061102,0167507,0176625 -}; -static unsigned short Q[] = { -/*0040200,0000000,0000000,0000000,*/ -0041123,0123257,0165741,0017142, -0041655,0065027,0173413,0115450, -0042261,0074011,0021573,0004150, -0042563,0166530,0013662,0007200, -0042743,0176427,0162443,0105214, -0043014,0062546,0153727,0123772, -0042717,0012470,0006227,0067424, -0042413,0061103,0003042,0013254 -}; -static unsigned short R[] = { -0040020,0067272,0101024,0155421, -0040243,0037467,0056706,0026462, -0040640,0116017,0120665,0034315, -0040705,0020162,0143350,0060137, -0040755,0016234,0134304,0130157, -0040476,0122700,0051070,0015473 -}; -static unsigned short S[] = { -/*0040200,0000000,0000000,0000000,*/ -0040420,0126200,0044276,0070413, -0041026,0053051,0007302,0063746, -0041100,0144203,0174051,0061151, -0041210,0123314,0126343,0177646, -0041031,0137125,0051431,0033011, -0040527,0117362,0152661,0066201 -}; -static unsigned short T[] = { -0041031,0126770,0170672,0166101, -0041664,0006522,0072360,0031770, -0043013,0100025,0162641,0126671, -0043332,0155231,0161627,0076200, -0044131,0024115,0021020,0117343 -}; -static unsigned short U[] = { -/*0040200,0000000,0000000,0000000,*/ -0041406,0037461,0177575,0032714, -0042402,0053350,0123061,0153557, -0043217,0111227,0032007,0164217, -0043660,0145000,0004013,0160114, -0044100,0071544,0167107,0125471 -}; -#define UTHRESH 14.0 -#endif - -#ifdef IBMPC -static unsigned short P[] = { -0x6479,0xa24f,0xeb24,0x3df0, -0x3488,0x4636,0x0dd7,0x3fe2, -0x6dc4,0xdec5,0xda53,0x401d, -0xba66,0xacad,0x518f,0x4048, -0x20f7,0xa9e0,0x90aa,0x4068, -0xcf58,0xc264,0x738f,0x4080, -0x34d8,0x6c74,0x343a,0x408d, -0x972a,0x21d6,0x0e35,0x4090, -0xffb3,0x5de8,0x6c48,0x4081 -}; -static unsigned short Q[] = { -/*0x0000,0x0000,0x0000,0x3ff0,*/ -0x23cc,0xfd7c,0x74d5,0x402a, -0x7365,0xfee1,0xad42,0x4055, -0x610d,0x246f,0x2f01,0x4076, -0x41d0,0x02f6,0x7dab,0x408e, -0x7151,0xfca4,0x7fa2,0x409c, -0xf4ff,0xdafa,0x8cac,0x40a1, -0xede2,0x0192,0xe2a7,0x4099, -0x42d6,0x60c4,0x6c48,0x4081 -}; -static unsigned short R[] = { -0x9b62,0x5042,0x0dd7,0x3fe2, -0xc5a6,0xebb8,0x67e6,0x3ff4, -0xa71a,0xf436,0x1381,0x4014, -0x0c0c,0x58dd,0xa40e,0x4018, -0x960e,0x9718,0xa393,0x401d, -0x0367,0x0a47,0xd4b8,0x4007 -}; -static unsigned short S[] = { -/*0x0000,0x0000,0x0000,0x3ff0,*/ -0xce21,0x0917,0x1590,0x4002, -0x4cfd,0x21d8,0xcac5,0x4022, -0x2c4d,0x7f05,0x1910,0x4028, -0x7ff5,0x959c,0x14d9,0x4031, -0x26c1,0xaa63,0x37ca,0x4023, -0x2d90,0x5ab6,0xf3de,0x400a -}; -static unsigned short T[] = { -0x5d88,0x1e37,0x35bf,0x4023, -0x067f,0x4e9e,0x81aa,0x4056, -0x35b7,0xbcb4,0x7002,0x40a1, -0xef90,0x3c72,0x5b53,0x40bb, -0x13dc,0xa442,0x2509,0x40eb -}; -static unsigned short U[] = { -/*0x0000,0x0000,0x0000,0x3ff0,*/ -0xa6ba,0x3fef,0xc7e6,0x4040, -0x3aee,0x14c6,0x4add,0x4080, -0xfd12,0xe680,0xf252,0x40b1, -0x7c0a,0x0101,0x1940,0x40d6, -0xf567,0x9dc8,0x0e6c,0x40e8 -}; -#define UTHRESH 37.519379347 -#endif - -#ifdef MIEEE -static unsigned short P[] = { -0x3df0,0xeb24,0xa24f,0x6479, -0x3fe2,0x0dd7,0x4636,0x3488, -0x401d,0xda53,0xdec5,0x6dc4, -0x4048,0x518f,0xacad,0xba66, -0x4068,0x90aa,0xa9e0,0x20f7, -0x4080,0x738f,0xc264,0xcf58, -0x408d,0x343a,0x6c74,0x34d8, -0x4090,0x0e35,0x21d6,0x972a, -0x4081,0x6c48,0x5de8,0xffb3 -}; -static unsigned short Q[] = { -0x402a,0x74d5,0xfd7c,0x23cc, -0x4055,0xad42,0xfee1,0x7365, -0x4076,0x2f01,0x246f,0x610d, -0x408e,0x7dab,0x02f6,0x41d0, -0x409c,0x7fa2,0xfca4,0x7151, -0x40a1,0x8cac,0xdafa,0xf4ff, -0x4099,0xe2a7,0x0192,0xede2, -0x4081,0x6c48,0x60c4,0x42d6 -}; -static unsigned short R[] = { -0x3fe2,0x0dd7,0x5042,0x9b62, -0x3ff4,0x67e6,0xebb8,0xc5a6, -0x4014,0x1381,0xf436,0xa71a, -0x4018,0xa40e,0x58dd,0x0c0c, -0x401d,0xa393,0x9718,0x960e, -0x4007,0xd4b8,0x0a47,0x0367 -}; -static unsigned short S[] = { -0x4002,0x1590,0x0917,0xce21, -0x4022,0xcac5,0x21d8,0x4cfd, -0x4028,0x1910,0x7f05,0x2c4d, -0x4031,0x14d9,0x959c,0x7ff5, -0x4023,0x37ca,0xaa63,0x26c1, -0x400a,0xf3de,0x5ab6,0x2d90 -}; -static unsigned short T[] = { -0x4023,0x35bf,0x1e37,0x5d88, -0x4056,0x81aa,0x4e9e,0x067f, -0x40a1,0x7002,0xbcb4,0x35b7, -0x40bb,0x5b53,0x3c72,0xef90, -0x40eb,0x2509,0xa442,0x13dc -}; -static unsigned short U[] = { -0x4040,0xc7e6,0x3fef,0xa6ba, -0x4080,0x4add,0x14c6,0x3aee, -0x40b1,0xf252,0xe680,0xfd12, -0x40d6,0x1940,0x0101,0x7c0a, -0x40e8,0x0e6c,0x9dc8,0xf567 -}; -#define UTHRESH 37.519379347 -#endif - -#ifdef ANSIPROT -extern double polevl ( double, void *, int ); -extern double p1evl ( double, void *, int ); -extern double exp ( double ); -extern double log ( double ); -extern double fabs ( double ); -double erf ( double ); -double erfc ( double ); -#else -double polevl(), p1evl(), exp(), log(), fabs(); -double erf(), erfc(); -#endif - -double ndtr(a) -double a; -{ -double x, y, z; - -x = a * SQRTH; -z = fabs(x); - -if( z < SQRTH ) - y = 0.5 + 0.5 * erf(x); - -else - { - y = 0.5 * erfc(z); - - if( x > 0 ) - y = 1.0 - y; - } - -return(y); -} - - -double erfc(a) -double a; -{ -double p,q,x,y,z; - - -if( a < 0.0 ) - x = -a; -else - x = a; - -if( x < 1.0 ) - return( 1.0 - erf(a) ); - -z = -a * a; - -if( z < -MAXLOG ) - { -under: - mtherr( "erfc", UNDERFLOW ); - if( a < 0 ) - return( 2.0 ); - else - return( 0.0 ); - } - -z = exp(z); - -if( x < 8.0 ) - { - p = polevl( x, P, 8 ); - q = p1evl( x, Q, 8 ); - } -else - { - p = polevl( x, R, 5 ); - q = p1evl( x, S, 6 ); - } -y = (z * p)/q; - -if( a < 0 ) - y = 2.0 - y; - -if( y == 0.0 ) - goto under; - -return(y); -} - - - -double erf(x) -double x; -{ -double y, z; - -if( fabs(x) > 1.0 ) - return( 1.0 - erfc(x) ); -z = x * x; -y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 ); -return( y ); - -} |