File Coverage

lib/PDL/Math/ndtri.c
Criterion Covered Total %
statement 23 28 82.1
branch 9 12 75.0
condition n/a
subroutine n/a
pod n/a
total 32 40 80.0


line stmt bran cond sub pod time code
1             /* ndtri.c
2             *
3             * Inverse of Normal distribution function
4             *
5             *
6             *
7             * SYNOPSIS:
8             *
9             * double x, y, ndtri();
10             *
11             * x = ndtri( 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             * DEC 0.125, 1 5500 9.5e-17 2.1e-17
35             * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
36             * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
37             * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
38             *
39             *
40             * ERROR MESSAGES:
41             *
42             * message condition value returned
43             * ndtri domain x <= 0 -MAXNUM
44             * ndtri domain x >= 1 MAXNUM
45             *
46             */
47            
48              
49             /*
50             Cephes Math Library Release 2.1: January, 1989
51             Copyright 1984, 1987, 1989 by Stephen L. Moshier
52             Direct inquiries to 30 Frost Street, Cambridge, MA 02140
53             */
54              
55             #include "mconf.h"
56             extern double MAXNUM;
57              
58             /* sqrt(2pi) */
59             static double s2pi = 2.50662827463100050242E0;
60              
61             /* approximation for 0 <= |y - 0.5| <= 3/8 */
62             static double P0[5] = {
63             -5.99633501014107895267E1,
64             9.80010754185999661536E1,
65             -5.66762857469070293439E1,
66             1.39312609387279679503E1,
67             -1.23916583867381258016E0,
68             };
69             static double Q0[8] = {
70             /* 1.00000000000000000000E0,*/
71             1.95448858338141759834E0,
72             4.67627912898881538453E0,
73             8.63602421390890590575E1,
74             -2.25462687854119370527E2,
75             2.00260212380060660359E2,
76             -8.20372256168333339912E1,
77             1.59056225126211695515E1,
78             -1.18331621121330003142E0,
79             };
80              
81              
82             /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
83             * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
84             */
85             static double P1[9] = {
86             4.05544892305962419923E0,
87             3.15251094599893866154E1,
88             5.71628192246421288162E1,
89             4.40805073893200834700E1,
90             1.46849561928858024014E1,
91             2.18663306850790267539E0,
92             -1.40256079171354495875E-1,
93             -3.50424626827848203418E-2,
94             -8.57456785154685413611E-4,
95             };
96             static double Q1[8] = {
97             /* 1.00000000000000000000E0,*/
98             1.57799883256466749731E1,
99             4.53907635128879210584E1,
100             4.13172038254672030440E1,
101             1.50425385692907503408E1,
102             2.50464946208309415979E0,
103             -1.42182922854787788574E-1,
104             -3.80806407691578277194E-2,
105             -9.33259480895457427372E-4,
106             };
107              
108             /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
109             * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
110             */
111              
112             static double P2[9] = {
113             3.23774891776946035970E0,
114             6.91522889068984211695E0,
115             3.93881025292474443415E0,
116             1.33303460815807542389E0,
117             2.01485389549179081538E-1,
118             1.23716634817820021358E-2,
119             3.01581553508235416007E-4,
120             2.65806974686737550832E-6,
121             6.23974539184983293730E-9,
122             };
123             static double Q2[8] = {
124             /* 1.00000000000000000000E0,*/
125             6.02427039364742014255E0,
126             3.67983563856160859403E0,
127             1.37702099489081330271E0,
128             2.16236993594496635890E-1,
129             1.34204006088543189037E-2,
130             3.28014464682127739104E-4,
131             2.89247864745380683936E-6,
132             6.79019408009981274425E-9,
133             };
134              
135             #ifndef ANSIPROT
136             double polevl(), p1evl(), log(), sqrt();
137             #endif
138              
139 28           double ndtri(double y0)
140             {
141             double x, y, z, y2, x0, x1;
142             int code;
143              
144 28 50         if( y0 <= 0.0 )
145             {
146 0           mtherr( "ndtri", DOMAIN );
147 0           return( -MAXNUM );
148             }
149 28 50         if( y0 >= 1.0 )
150             {
151 0           mtherr( "ndtri", DOMAIN );
152 0           return( MAXNUM );
153             }
154 28           code = 1;
155 28           y = y0;
156 28 100         if( y > (1.0 - 0.13533528323661269189) ) /* 0.135... = exp(-2) */
157             {
158 2           y = 1.0 - y;
159 2           code = 0;
160             }
161              
162 28 100         if( y > 0.13533528323661269189 )
163             {
164 22           y = y - 0.5;
165 22           y2 = y * y;
166 22           x = y + y * (y2 * polevl( y2, P0, 4)/p1evl( y2, Q0, 8 ));
167 22           x = x * s2pi;
168 22           return(x);
169             }
170              
171 6           x = sqrt( -2.0 * log(y) );
172 6           x0 = x - log(x)/x;
173              
174 6           z = 1.0/x;
175 6 50         if( x < 8.0 ) /* y > exp(-32) = 1.2664165549e-14 */
176 6           x1 = z * polevl( z, P1, 8 )/p1evl( z, Q1, 8 );
177             else
178 0           x1 = z * polevl( z, P2, 8 )/p1evl( z, Q2, 8 );
179 6           x = x0 - x1;
180 6 100         if( code != 0 )
181 4           x = -x;
182 6           return( x );
183             }