lib/libm/ndtr.c

   1 /*                                                      ndtr.c
   2  *
   3  *      Normal distribution function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, ndtr();
  10  *
  11  * y = ndtr( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the area under the Gaussian probability density
  18  * function, integrated from minus infinity to x:
  19  *
  20  *                            x
  21  *                             -
  22  *                   1        | |          2
  23  *    ndtr(x)  = ---------    |    exp( - t /2 ) dt
  24  *               sqrt(2pi)  | |
  25  *                           -
  26  *                          -inf.
  27  *
  28  *             =  ( 1 + erf(z) ) / 2
  29  *             =  erfc(z) / 2
  30  *
  31  * where z = x/sqrt(2). Computation is via the functions
  32  * erf and erfc with care to avoid error amplification in computing exp(-x^2).
  33  *
  34  *
  35  * ACCURACY:
  36  *
  37  *                      Relative error:
  38  * arithmetic   domain     # trials      peak         rms
  39  *    IEEE     -13,0        30000       1.3e-15     2.2e-16
  40  *
  41  *
  42  * ERROR MESSAGES:
  43  *
  44  *   message         condition         value returned
  45  * erfc underflow    x > 37.519379347       0.0
  46  *
  47  */
  48 \f/*                                                     erf.c
  49  *
  50  *      Error function
  51  *
  52  *
  53  *
  54  * SYNOPSIS:
  55  *
  56  * double x, y, erf();
  57  *
  58  * y = erf( x );
  59  *
  60  *
  61  *
  62  * DESCRIPTION:
  63  *
  64  * The integral is
  65  *
  66  *                           x
  67  *                            -
  68  *                 2         | |          2
  69  *   erf(x)  =  --------     |    exp( - t  ) dt.
  70  *              sqrt(pi)   | |
  71  *                          -
  72  *                           0
  73  *
  74  * The magnitude of x is limited to 9.231948545 for DEC
  75  * arithmetic; 1 or -1 is returned outside this range.
  76  *
  77  * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
  78  * erf(x) = 1 - erfc(x).
  79  *
  80  *
  81  *
  82  * ACCURACY:
  83  *
  84  *                      Relative error:
  85  * arithmetic   domain     # trials      peak         rms
  86  *    DEC       0,1         14000       4.7e-17     1.5e-17
  87  *    IEEE      0,1         30000       3.7e-16     1.0e-16
  88  *
  89  */
  90 \f/*                                                     erfc.c
  91  *
  92  *      Complementary error function
  93  *
  94  *
  95  *
  96  * SYNOPSIS:
  97  *
  98  * double x, y, erfc();
  99  *
 100  * y = erfc( x );
 101  *
 102  *
 103  *
 104  * DESCRIPTION:
 105  *
 106  *
 107  *  1 - erf(x) =
 108  *
 109  *                           inf.
 110  *                             -
 111  *                  2         | |          2
 112  *   erfc(x)  =  --------     |    exp( - t  ) dt
 113  *               sqrt(pi)   | |
 114  *                           -
 115  *                            x
 116  *
 117  *
 118  * For small x, erfc(x) = 1 - erf(x); otherwise rational
 119  * approximations are computed.
 120  *
 121  * A special function expx2l.c is used to suppress error amplification
 122  * in computing exp(-x^2).
 123  *
 124  *
 125  * ACCURACY:
 126  *
 127  *                      Relative error:
 128  * arithmetic   domain     # trials      peak         rms
 129  *    IEEE      0,26.6417   30000       1.3e-15     2.2e-16
 130  *
 131  *
 132  * ERROR MESSAGES:
 133  *
 134  *   message         condition              value returned
 135  * erfc underflow    x > 9.231948545 (DEC)       0.0
 136  *
 137  *
 138  */
 139 \f
 140
 141 /*
 142 Cephes Math Library Release 2.9:  November, 2000
 143 Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
 144 */
 145
 146
 147 #include "mconf.h"
 148
 149 extern double SQRTH;
 150 extern double MAXLOG;
 151
 152
 153 #ifdef UNK
 154 static double P[] = {
 155  2.46196981473530512524E-10,
 156  5.64189564831068821977E-1,
 157  7.46321056442269912687E0,
 158  4.86371970985681366614E1,
 159  1.96520832956077098242E2,
 160  5.26445194995477358631E2,
 161  9.34528527171957607540E2,
 162  1.02755188689515710272E3,
 163  5.57535335369399327526E2
 164 };
 165 static double Q[] = {
 166 /* 1.00000000000000000000E0,*/
 167  1.32281951154744992508E1,
 168  8.67072140885989742329E1,
 169  3.54937778887819891062E2,
 170  9.75708501743205489753E2,
 171  1.82390916687909736289E3,
 172  2.24633760818710981792E3,
 173  1.65666309194161350182E3,
 174  5.57535340817727675546E2
 175 };
 176 static double R[] = {
 177  5.64189583547755073984E-1,
 178  1.27536670759978104416E0,
 179  5.01905042251180477414E0,
 180  6.16021097993053585195E0,
 181  7.40974269950448939160E0,
 182  2.97886665372100240670E0
 183 };
 184 static double S[] = {
 185 /* 1.00000000000000000000E0,*/
 186  2.26052863220117276590E0,
 187  9.39603524938001434673E0,
 188  1.20489539808096656605E1,
 189  1.70814450747565897222E1,
 190  9.60896809063285878198E0,
 191  3.36907645100081516050E0
 192 };
 193 static double T[] = {
 194  9.60497373987051638749E0,
 195  9.00260197203842689217E1,
 196  2.23200534594684319226E3,
 197  7.00332514112805075473E3,
 198  5.55923013010394962768E4
 199 };
 200 static double U[] = {
 201 /* 1.00000000000000000000E0,*/
 202  3.35617141647503099647E1,
 203  5.21357949780152679795E2,
 204  4.59432382970980127987E3,
 205  2.26290000613890934246E4,
 206  4.92673942608635921086E4
 207 };
 208
 209 #define UTHRESH 37.519379347
 210 #endif
 211
 212 #ifdef DEC
 213 static unsigned short P[] = {
 214 0030207,0054445,0011173,0021706,
 215 0040020,0067272,0030661,0122075,
 216 0040756,0151236,0173053,0067042,
 217 0041502,0106175,0062555,0151457,
 218 0042104,0102525,0047401,0003667,
 219 0042403,0116176,0011446,0075303,
 220 0042551,0120723,0061641,0123275,
 221 0042600,0070651,0007264,0134516,
 222 0042413,0061102,0167507,0176625
 223 };
 224 static unsigned short Q[] = {
 225 /*0040200,0000000,0000000,0000000,*/
 226 0041123,0123257,0165741,0017142,
 227 0041655,0065027,0173413,0115450,
 228 0042261,0074011,0021573,0004150,
 229 0042563,0166530,0013662,0007200,
 230 0042743,0176427,0162443,0105214,
 231 0043014,0062546,0153727,0123772,
 232 0042717,0012470,0006227,0067424,
 233 0042413,0061103,0003042,0013254
 234 };
 235 static unsigned short R[] = {
 236 0040020,0067272,0101024,0155421,
 237 0040243,0037467,0056706,0026462,
 238 0040640,0116017,0120665,0034315,
 239 0040705,0020162,0143350,0060137,
 240 0040755,0016234,0134304,0130157,
 241 0040476,0122700,0051070,0015473
 242 };
 243 static unsigned short S[] = {
 244 /*0040200,0000000,0000000,0000000,*/
 245 0040420,0126200,0044276,0070413,
 246 0041026,0053051,0007302,0063746,
 247 0041100,0144203,0174051,0061151,
 248 0041210,0123314,0126343,0177646,
 249 0041031,0137125,0051431,0033011,
 250 0040527,0117362,0152661,0066201
 251 };
 252 static unsigned short T[] = {
 253 0041031,0126770,0170672,0166101,
 254 0041664,0006522,0072360,0031770,
 255 0043013,0100025,0162641,0126671,
 256 0043332,0155231,0161627,0076200,
 257 0044131,0024115,0021020,0117343
 258 };
 259 static unsigned short U[] = {
 260 /*0040200,0000000,0000000,0000000,*/
 261 0041406,0037461,0177575,0032714,
 262 0042402,0053350,0123061,0153557,
 263 0043217,0111227,0032007,0164217,
 264 0043660,0145000,0004013,0160114,
 265 0044100,0071544,0167107,0125471
 266 };
 267 #define UTHRESH 14.0
 268 #endif
 269
 270 #ifdef IBMPC
 271 static unsigned short P[] = {
 272 0x6479,0xa24f,0xeb24,0x3df0,
 273 0x3488,0x4636,0x0dd7,0x3fe2,
 274 0x6dc4,0xdec5,0xda53,0x401d,
 275 0xba66,0xacad,0x518f,0x4048,
 276 0x20f7,0xa9e0,0x90aa,0x4068,
 277 0xcf58,0xc264,0x738f,0x4080,
 278 0x34d8,0x6c74,0x343a,0x408d,
 279 0x972a,0x21d6,0x0e35,0x4090,
 280 0xffb3,0x5de8,0x6c48,0x4081
 281 };
 282 static unsigned short Q[] = {
 283 /*0x0000,0x0000,0x0000,0x3ff0,*/
 284 0x23cc,0xfd7c,0x74d5,0x402a,
 285 0x7365,0xfee1,0xad42,0x4055,
 286 0x610d,0x246f,0x2f01,0x4076,
 287 0x41d0,0x02f6,0x7dab,0x408e,
 288 0x7151,0xfca4,0x7fa2,0x409c,
 289 0xf4ff,0xdafa,0x8cac,0x40a1,
 290 0xede2,0x0192,0xe2a7,0x4099,
 291 0x42d6,0x60c4,0x6c48,0x4081
 292 };
 293 static unsigned short R[] = {
 294 0x9b62,0x5042,0x0dd7,0x3fe2,
 295 0xc5a6,0xebb8,0x67e6,0x3ff4,
 296 0xa71a,0xf436,0x1381,0x4014,
 297 0x0c0c,0x58dd,0xa40e,0x4018,
 298 0x960e,0x9718,0xa393,0x401d,
 299 0x0367,0x0a47,0xd4b8,0x4007
 300 };
 301 static unsigned short S[] = {
 302 /*0x0000,0x0000,0x0000,0x3ff0,*/
 303 0xce21,0x0917,0x1590,0x4002,
 304 0x4cfd,0x21d8,0xcac5,0x4022,
 305 0x2c4d,0x7f05,0x1910,0x4028,
 306 0x7ff5,0x959c,0x14d9,0x4031,
 307 0x26c1,0xaa63,0x37ca,0x4023,
 308 0x2d90,0x5ab6,0xf3de,0x400a
 309 };
 310 static unsigned short T[] = {
 311 0x5d88,0x1e37,0x35bf,0x4023,
 312 0x067f,0x4e9e,0x81aa,0x4056,
 313 0x35b7,0xbcb4,0x7002,0x40a1,
 314 0xef90,0x3c72,0x5b53,0x40bb,
 315 0x13dc,0xa442,0x2509,0x40eb
 316 };
 317 static unsigned short U[] = {
 318 /*0x0000,0x0000,0x0000,0x3ff0,*/
 319 0xa6ba,0x3fef,0xc7e6,0x4040,
 320 0x3aee,0x14c6,0x4add,0x4080,
 321 0xfd12,0xe680,0xf252,0x40b1,
 322 0x7c0a,0x0101,0x1940,0x40d6,
 323 0xf567,0x9dc8,0x0e6c,0x40e8
 324 };
 325 #define UTHRESH 37.519379347
 326 #endif
 327
 328 #ifdef MIEEE
 329 static unsigned short P[] = {
 330 0x3df0,0xeb24,0xa24f,0x6479,
 331 0x3fe2,0x0dd7,0x4636,0x3488,
 332 0x401d,0xda53,0xdec5,0x6dc4,
 333 0x4048,0x518f,0xacad,0xba66,
 334 0x4068,0x90aa,0xa9e0,0x20f7,
 335 0x4080,0x738f,0xc264,0xcf58,
 336 0x408d,0x343a,0x6c74,0x34d8,
 337 0x4090,0x0e35,0x21d6,0x972a,
 338 0x4081,0x6c48,0x5de8,0xffb3
 339 };
 340 static unsigned short Q[] = {
 341 0x402a,0x74d5,0xfd7c,0x23cc,
 342 0x4055,0xad42,0xfee1,0x7365,
 343 0x4076,0x2f01,0x246f,0x610d,
 344 0x408e,0x7dab,0x02f6,0x41d0,
 345 0x409c,0x7fa2,0xfca4,0x7151,
 346 0x40a1,0x8cac,0xdafa,0xf4ff,
 347 0x4099,0xe2a7,0x0192,0xede2,
 348 0x4081,0x6c48,0x60c4,0x42d6
 349 };
 350 static unsigned short R[] = {
 351 0x3fe2,0x0dd7,0x5042,0x9b62,
 352 0x3ff4,0x67e6,0xebb8,0xc5a6,
 353 0x4014,0x1381,0xf436,0xa71a,
 354 0x4018,0xa40e,0x58dd,0x0c0c,
 355 0x401d,0xa393,0x9718,0x960e,
 356 0x4007,0xd4b8,0x0a47,0x0367
 357 };
 358 static unsigned short S[] = {
 359 0x4002,0x1590,0x0917,0xce21,
 360 0x4022,0xcac5,0x21d8,0x4cfd,
 361 0x4028,0x1910,0x7f05,0x2c4d,
 362 0x4031,0x14d9,0x959c,0x7ff5,
 363 0x4023,0x37ca,0xaa63,0x26c1,
 364 0x400a,0xf3de,0x5ab6,0x2d90
 365 };
 366 static unsigned short T[] = {
 367 0x4023,0x35bf,0x1e37,0x5d88,
 368 0x4056,0x81aa,0x4e9e,0x067f,
 369 0x40a1,0x7002,0xbcb4,0x35b7,
 370 0x40bb,0x5b53,0x3c72,0xef90,
 371 0x40eb,0x2509,0xa442,0x13dc
 372 };
 373 static unsigned short U[] = {
 374 0x4040,0xc7e6,0x3fef,0xa6ba,
 375 0x4080,0x4add,0x14c6,0x3aee,
 376 0x40b1,0xf252,0xe680,0xfd12,
 377 0x40d6,0x1940,0x0101,0x7c0a,
 378 0x40e8,0x0e6c,0x9dc8,0xf567
 379 };
 380 #define UTHRESH 37.519379347
 381 #endif
 382
 383 #ifdef ANSIPROT
 384 extern double polevl ( double, void *, int );
 385 extern double p1evl ( double, void *, int );
 386 extern double exp ( double );
 387 extern double log ( double );
 388 extern double fabs ( double );
 389 extern double sqrt ( double );
 390 extern double expx2 ( double, int );
 391 double erf ( double );
 392 double erfc ( double );
 393 static double erfce ( double );
 394 #else
 395 double polevl(), p1evl(), exp(), log(), fabs();
 396 double erf(), erfc(), expx2(), sqrt();
 397 static double erfce();
 398 #endif
 399
 400 double ndtr(a)
 401 double a;
 402 {
 403 double x, y, z;
 404
 405 x = a * SQRTH;
 406 z = fabs(x);
 407
 408 /* if( z < SQRTH ) */
 409 if( z < 1.0 )
 410         y = 0.5 + 0.5 * erf(x);
 411
 412 else
 413         {
 414         /* See below for erfce. */
 415         y = 0.5 * erfce(z);
 416         /* Multiply by exp(-x^2 / 2)  */
 417         z = expx2(a, -1);
 418         y = y * sqrt(z);
 419         if( x > 0 )
 420                 y = 1.0 - y;
 421         }
 422
 423 return(y);
 424 }
 425
 426
 427 double erfc(a)
 428 double a;
 429 {
 430 double p,q,x,y,z;
 431
 432
 433 if( a < 0.0 )
 434         x = -a;
 435 else
 436         x = a;
 437
 438 if( x < 1.0 )
 439         return( 1.0 - erf(a) );
 440
 441 z = -a * a;
 442
 443 if( z < -MAXLOG )
 444         {
 445 under:
 446         mtherr( "erfc", UNDERFLOW );
 447         if( a < 0 )
 448                 return( 2.0 );
 449         else
 450                 return( 0.0 );
 451         }
 452
 453 /* Compute z = exp(z).  */
 454 /* z = exp(z); */
 455 z = expx2(a, -1);
 456
 457 if( x < 8.0 )
 458         {
 459         p = polevl( x, P, 8 );
 460         q = p1evl( x, Q, 8 );
 461         }
 462 else
 463         {
 464         p = polevl( x, R, 5 );
 465         q = p1evl( x, S, 6 );
 466         }
 467 y = (z * p)/q;
 468
 469 if( a < 0 )
 470         y = 2.0 - y;
 471
 472 if( y == 0.0 )
 473         goto under;
 474
 475 return(y);
 476 }
 477
 478
 479 /* Exponentially scaled erfc function
 480    exp(x^2) erfc(x)
 481    valid for x > 1.
 482    Use with ndtr and expx2.  */
 483 static double erfce(x)
 484 double x;
 485 {
 486 double p,q;
 487
 488 if( x < 8.0 )
 489         {
 490         p = polevl( x, P, 8 );
 491         q = p1evl( x, Q, 8 );
 492         }
 493 else
 494         {
 495         p = polevl( x, R, 5 );
 496         q = p1evl( x, S, 6 );
 497         }
 498 return (p/q);
 499 }
 500
 501
 502
 503 double erf(x)
 504 double x;
 505 {
 506 double y, z;
 507
 508 if( fabs(x) > 1.0 )
 509         return( 1.0 - erfc(x) );
 510 z = x * x;
 511 y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
 512 return( y );
 513
 514 }