lib/libm/ndtrif.c

   1 /*                                                      ndtrif.c
   2  *
   3  *      Inverse of Normal distribution function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, ndtrif();
  10  *
  11  * x = ndtrif( y );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns the argument, x, for which the area under the
  18  * Gaussian probability density function (integrated from
  19  * minus infinity to x) is equal to y.
  20  *
  21  *
  22  * For small arguments 0 < y < exp(-2), the program computes
  23  * z = sqrt( -2.0 * log(y) );  then the approximation is
  24  * x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).
  25  * There are two rational functions P/Q, one for 0 < y < exp(-32)
  26  * and the other for y up to exp(-2).  For larger arguments,
  27  * w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  *                      Relative error:
  33  * arithmetic   domain        # trials      peak         rms
  34  *    IEEE     1e-38, 1        30000       3.6e-7      5.0e-8
  35  *
  36  *
  37  * ERROR MESSAGES:
  38  *
  39  *   message         condition    value returned
  40  * ndtrif domain      x <= 0        -MAXNUM
  41  * ndtrif domain      x >= 1         MAXNUM
  42  *
  43  */
  44 \f
  45
  46 /*
  47 Cephes Math Library Release 2.2:  July, 1992
  48 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  49 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  50 */
  51
  52 #include "mconf.h"
  53 extern float MAXNUMF;
  54
  55 /* sqrt(2pi) */
  56 static float s2pi = 2.50662827463100050242;
  57
  58 /* approximation for 0 <= |y - 0.5| <= 3/8 */
  59 static float P0[5] = {
  60 -5.99633501014107895267E1,
  61  9.80010754185999661536E1,
  62 -5.66762857469070293439E1,
  63  1.39312609387279679503E1,
  64 -1.23916583867381258016E0,
  65 };
  66 static float Q0[8] = {
  67 /* 1.00000000000000000000E0,*/
  68  1.95448858338141759834E0,
  69  4.67627912898881538453E0,
  70  8.63602421390890590575E1,
  71 -2.25462687854119370527E2,
  72  2.00260212380060660359E2,
  73 -8.20372256168333339912E1,
  74  1.59056225126211695515E1,
  75 -1.18331621121330003142E0,
  76 };
  77
  78 /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
  79  * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
  80  */
  81 static float P1[9] = {
  82  4.05544892305962419923E0,
  83  3.15251094599893866154E1,
  84  5.71628192246421288162E1,
  85  4.40805073893200834700E1,
  86  1.46849561928858024014E1,
  87  2.18663306850790267539E0,
  88 -1.40256079171354495875E-1,
  89 -3.50424626827848203418E-2,
  90 -8.57456785154685413611E-4,
  91 };
  92 static float Q1[8] = {
  93 /*  1.00000000000000000000E0,*/
  94  1.57799883256466749731E1,
  95  4.53907635128879210584E1,
  96  4.13172038254672030440E1,
  97  1.50425385692907503408E1,
  98  2.50464946208309415979E0,
  99 -1.42182922854787788574E-1,
 100 -3.80806407691578277194E-2,
 101 -9.33259480895457427372E-4,
 102 };
 103
 104
 105 /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
 106  * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
 107  */
 108
 109 static float P2[9] = {
 110   3.23774891776946035970E0,
 111   6.91522889068984211695E0,
 112   3.93881025292474443415E0,
 113   1.33303460815807542389E0,
 114   2.01485389549179081538E-1,
 115   1.23716634817820021358E-2,
 116   3.01581553508235416007E-4,
 117   2.65806974686737550832E-6,
 118   6.23974539184983293730E-9,
 119 };
 120 static float Q2[8] = {
 121 /*  1.00000000000000000000E0,*/
 122   6.02427039364742014255E0,
 123   3.67983563856160859403E0,
 124   1.37702099489081330271E0,
 125   2.16236993594496635890E-1,
 126   1.34204006088543189037E-2,
 127   3.28014464682127739104E-4,
 128   2.89247864745380683936E-6,
 129   6.79019408009981274425E-9,
 130 };
 131
 132 #ifdef ANSIC
 133 float polevlf(float, float *, int);
 134 float p1evlf(float, float *, int);
 135 float logf(float), sqrtf(float);
 136 #else
 137 float polevlf(), p1evlf(), logf(), sqrtf();
 138 #endif
 139
 140
 141 #ifdef ANSIC
 142 float ndtrif(float yy0)
 143 #else
 144 float ndtrif(yy0)
 145 double yy0;
 146 #endif
 147 {
 148 float y0, x, y, z, y2, x0, x1;
 149 int code;
 150
 151 y0 = yy0;
 152 if( y0 <= 0.0 )
 153         {
 154         mtherr( "ndtrif", DOMAIN );
 155         return( -MAXNUMF );
 156         }
 157 if( y0 >= 1.0 )
 158         {
 159         mtherr( "ndtrif", DOMAIN );
 160         return( MAXNUMF );
 161         }
 162 code = 1;
 163 y = y0;
 164 if( y > (1.0 - 0.13533528323661269189) ) /* 0.135... = exp(-2) */
 165         {
 166         y = 1.0 - y;
 167         code = 0;
 168         }
 169
 170 if( y > 0.13533528323661269189 )
 171         {
 172         y = y - 0.5;
 173         y2 = y * y;
 174         x = y + y * (y2 * polevlf( y2, P0, 4)/p1evlf( y2, Q0, 8 ));
 175         x = x * s2pi;
 176         return(x);
 177         }
 178
 179 x = sqrtf( -2.0 * logf(y) );
 180 x0 = x - logf(x)/x;
 181
 182 z = 1.0/x;
 183 if( x < 8.0 ) /* y > exp(-32) = 1.2664165549e-14 */
 184         x1 = z * polevlf( z, P1, 8 )/p1evlf( z, Q1, 8 );
 185 else
 186         x1 = z * polevlf( z, P2, 8 )/p1evlf( z, Q2, 8 );
 187 x = x0 - x1;
 188 if( code != 0 )
 189         x = -x;
 190 return( x );
 191 }