libm/float/ndtrf.c

   1 /*                                                      ndtrf.c
   2  *
   3  *      Normal distribution function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, ndtrf();
  10  *
  11  * y = ndtrf( 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.
  33  *
  34  *
  35  * ACCURACY:
  36  *
  37  *                      Relative error:
  38  * arithmetic   domain     # trials      peak         rms
  39  *    IEEE     -13,0        50000       1.5e-5      2.6e-6
  40  *
  41  *
  42  * ERROR MESSAGES:
  43  *
  44  * See erfcf().
  45  *
  46  */
  47 \f/*                                                     erff.c
  48  *
  49  *      Error function
  50  *
  51  *
  52  *
  53  * SYNOPSIS:
  54  *
  55  * float x, y, erff();
  56  *
  57  * y = erff( x );
  58  *
  59  *
  60  *
  61  * DESCRIPTION:
  62  *
  63  * The integral is
  64  *
  65  *                           x
  66  *                            -
  67  *                 2         | |          2
  68  *   erf(x)  =  --------     |    exp( - t  ) dt.
  69  *              sqrt(pi)   | |
  70  *                          -
  71  *                           0
  72  *
  73  * The magnitude of x is limited to 9.231948545 for DEC
  74  * arithmetic; 1 or -1 is returned outside this range.
  75  *
  76  * For 0 <= |x| < 1, erf(x) = x * P(x**2); otherwise
  77  * erf(x) = 1 - erfc(x).
  78  *
  79  *
  80  *
  81  * ACCURACY:
  82  *
  83  *                      Relative error:
  84  * arithmetic   domain     # trials      peak         rms
  85  *    IEEE      -9.3,9.3    50000       1.7e-7      2.8e-8
  86  *
  87  */
  88 \f/*                                                     erfcf.c
  89  *
  90  *      Complementary error function
  91  *
  92  *
  93  *
  94  * SYNOPSIS:
  95  *
  96  * float x, y, erfcf();
  97  *
  98  * y = erfcf( x );
  99  *
 100  *
 101  *
 102  * DESCRIPTION:
 103  *
 104  *
 105  *  1 - erf(x) =
 106  *
 107  *                           inf.
 108  *                             -
 109  *                  2         | |          2
 110  *   erfc(x)  =  --------     |    exp( - t  ) dt
 111  *               sqrt(pi)   | |
 112  *                           -
 113  *                            x
 114  *
 115  *
 116  * For small x, erfc(x) = 1 - erf(x); otherwise polynomial
 117  * approximations 1/x P(1/x**2) are computed.
 118  *
 119  *
 120  *
 121  * ACCURACY:
 122  *
 123  *                      Relative error:
 124  * arithmetic   domain     # trials      peak         rms
 125  *    IEEE      -9.3,9.3    50000       3.9e-6      7.2e-7
 126  *
 127  *
 128  * ERROR MESSAGES:
 129  *
 130  *   message           condition              value returned
 131  * erfcf underflow    x**2 > MAXLOGF              0.0
 132  *
 133  *
 134  */
 135 \f
 136
 137 /*
 138 Cephes Math Library Release 2.2:  June, 1992
 139 Copyright 1984, 1987, 1988 by Stephen L. Moshier
 140 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 141 */
 142
 143
 144 #include <math.h>
 145
 146 extern float MAXLOGF, SQRTHF;
 147
 148
 149 /* erfc(x) = exp(-x^2) P(1/x), 1 < x < 2 */
 150 static float P[] = {
 151  2.326819970068386E-002,
 152 -1.387039388740657E-001,
 153  3.687424674597105E-001,
 154 -5.824733027278666E-001,
 155  6.210004621745983E-001,
 156 -4.944515323274145E-001,
 157  3.404879937665872E-001,
 158 -2.741127028184656E-001,
 159  5.638259427386472E-001
 160 };
 161
 162 /* erfc(x) = exp(-x^2) 1/x P(1/x^2), 2 < x < 14 */
 163 static float R[] = {
 164 -1.047766399936249E+001,
 165  1.297719955372516E+001,
 166 -7.495518717768503E+000,
 167  2.921019019210786E+000,
 168 -1.015265279202700E+000,
 169  4.218463358204948E-001,
 170 -2.820767439740514E-001,
 171  5.641895067754075E-001
 172 };
 173
 174 /* erf(x) = x P(x^2), 0 < x < 1 */
 175 static float T[] = {
 176  7.853861353153693E-005,
 177 -8.010193625184903E-004,
 178  5.188327685732524E-003,
 179 -2.685381193529856E-002,
 180  1.128358514861418E-001,
 181 -3.761262582423300E-001,
 182  1.128379165726710E+000
 183 };
 184
 185 /*#define UTHRESH 37.519379347*/
 186
 187 #define UTHRESH 14.0
 188
 189 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
 190
 191 #ifdef ANSIC
 192 float polevlf(float, float *, int);
 193 float expf(float), logf(float), erff(float), erfcf(float);
 194 #else
 195 float polevlf(), expf(), logf(), erff(), erfcf();
 196 #endif
 197
 198
 199
 200 float ndtrf(float aa)
 201 {
 202 float x, y, z;
 203
 204 x = aa;
 205 x *= SQRTHF;
 206 z = fabsf(x);
 207
 208 if( z < SQRTHF )
 209         y = 0.5 + 0.5 * erff(x);
 210 else
 211         {
 212         y = 0.5 * erfcf(z);
 213
 214         if( x > 0 )
 215                 y = 1.0 - y;
 216         }
 217
 218 return(y);
 219 }
 220
 221
 222 float erfcf(float aa)
 223 {
 224 float a, p,q,x,y,z;
 225
 226
 227 a = aa;
 228 x = fabsf(a);
 229
 230 if( x < 1.0 )
 231         return( 1.0 - erff(a) );
 232
 233 z = -a * a;
 234
 235 if( z < -MAXLOGF )
 236         {
 237 under:
 238         mtherr( "erfcf", UNDERFLOW );
 239         if( a < 0 )
 240                 return( 2.0 );
 241         else
 242                 return( 0.0 );
 243         }
 244
 245 z = expf(z);
 246 q = 1.0/x;
 247 y = q * q;
 248 if( x < 2.0 )
 249         {
 250         p = polevlf( y, P, 8 );
 251         }
 252 else
 253         {
 254         p = polevlf( y, R, 7 );
 255         }
 256
 257 y = z * q * p;
 258
 259 if( a < 0 )
 260         y = 2.0 - y;
 261
 262 if( y == 0.0 )
 263         goto under;
 264
 265 return(y);
 266 }
 267
 268
 269 float erff(float xx)
 270 {
 271 float x, y, z;
 272
 273 x = xx;
 274 if( fabsf(x) > 1.0 )
 275         return( 1.0 - erfcf(x) );
 276
 277 z = x * x;
 278 y = x * polevlf( z, T, 6 );
 279 return( y );
 280
 281 }