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.8: June, 2000
51 Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
52 */
53
54 #include "mconf.h"
55
56 /* sqrt(2pi) */
57 static double s2pi = 2.50662827463100050242E0;
58
59 /* approximation for 0 <= |y - 0.5| <= 3/8 */
60 static double P0[5] = {
61 -5.99633501014107895267E1,
62 9.80010754185999661536E1,
63 -5.66762857469070293439E1,
64 1.39312609387279679503E1,
65 -1.23916583867381258016E0,
66 };
67
68 static double Q0[8] = {
69 /* 1.00000000000000000000E0,*/
70 1.95448858338141759834E0,
71 4.67627912898881538453E0,
72 8.63602421390890590575E1,
73 -2.25462687854119370527E2,
74 2.00260212380060660359E2,
75 -8.20372256168333339912E1,
76 1.59056225126211695515E1,
77 -1.18331621121330003142E0,
78 };
79
80 /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
81 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
82 */
83
84 static double P1[9] = {
85 4.05544892305962419923E0,
86 3.15251094599893866154E1,
87 5.71628192246421288162E1,
88 4.40805073893200834700E1,
89 1.46849561928858024014E1,
90 2.18663306850790267539E0,
91 -1.40256079171354495875E-1,
92 -3.50424626827848203418E-2,
93 -8.57456785154685413611E-4,
94 };
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
124 static double Q2[8] = {
125 /* 1.00000000000000000000E0,*/
126 6.02427039364742014255E0,
127 3.67983563856160859403E0,
128 1.37702099489081330271E0,
129 2.16236993594496635890E-1,
130 1.34204006088543189037E-2,
131 3.28014464682127739104E-4,
132 2.89247864745380683936E-6,
133 6.79019408009981274425E-9,
134 };
135
ndtri(double y0)136 double ndtri (double y0)
137 {
138 double expm2 = 0.13533528323661269189; /* exp(-2) */
139 double x, y, z, y2, x0, x1;
140 int code;
141
142 if (y0 <= 0.0) {
143 mtherr_with_arg("ndtri", CEPHES_DOMAIN, y0);
144 return -MAXNUM;
145 }
146
147 if (y0 >= 1.0) {
148 mtherr_with_arg("ndtri", CEPHES_DOMAIN, y0);
149 return MAXNUM;
150 }
151
152 code = 1;
153 y = y0;
154
155 if (y > (1.0 - expm2)) {
156 y = 1.0 - y;
157 code = 0;
158 }
159
160 if (y > expm2) {
161 y = y - 0.5;
162 y2 = y * y;
163 x = y + y * (y2 * polevl(y2, P0, 4)/p1evl(y2, Q0, 8));
164 x = x * s2pi;
165 return x;
166 }
167
168 x = sqrt( -2.0 * log(y) );
169 x0 = x - log(x)/x;
170
171 z = 1.0/x;
172
173 if (x < 8.0) {
174 /* y > exp(-32) = 1.2664165549e-14 */
175 x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
176 } else {
177 x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
178 }
179
180 x = x0 - x1;
181
182 if (code != 0) {
183 x = -x;
184 }
185
186 return x;
187 }
188