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 #include "cephes.h"
57
58 extern double MAXNUM;
59
60 #ifdef UNK
61 /* sqrt(2pi) */
62 static double s2pi = 2.50662827463100050242E0;
63 #endif
64
65 #ifdef DEC
66 static union us2d_t s2p = {{0040440,0066230,0177661,0034055}};
67 #define s2pi s2p.d
68 #endif
69
70 #ifdef IBMPC
71 static union us2d_t s2p = {{0x2706,0x1ff6,0x0d93,0x4004}};
72 #define s2pi s2p.d
73 #endif
74
75 #ifdef MIEEE
76 static union us2d_t s2p = { {0x4004,0x0d93,0x1ff6,0x2706} };
77 #define s2pi s2p.d
78 #endif
79
80
81 /* approximation for 0 <= |y - 0.5| <= 3/8 */
82 #ifdef UNK
83 static double P0[5] = {
84 -5.99633501014107895267E1,
85 9.80010754185999661536E1,
86 -5.66762857469070293439E1,
87 1.39312609387279679503E1,
88 -1.23916583867381258016E0,
89 };
90 static double Q0[8] = {
91 /* 1.00000000000000000000E0,*/
92 1.95448858338141759834E0,
93 4.67627912898881538453E0,
94 8.63602421390890590575E1,
95 -2.25462687854119370527E2,
96 2.00260212380060660359E2,
97 -8.20372256168333339912E1,
98 1.59056225126211695515E1,
99 -1.18331621121330003142E0,
100 };
101 #endif
102 #ifdef DEC
103 static unsigned short P0[20] = {
104 0141557,0155170,0071360,0120550,
105 0041704,0000214,0172417,0067307,
106 0141542,0132204,0040066,0156723,
107 0041136,0163161,0157276,0007747,
108 0140236,0116374,0073666,0051764,
109 };
110 static unsigned short Q0[32] = {
111 /*0040200,0000000,0000000,0000000,*/
112 0040372,0026256,0110403,0123707,
113 0040625,0122024,0020277,0026661,
114 0041654,0134161,0124134,0007244,
115 0142141,0073162,0133021,0131371,
116 0042110,0041235,0043516,0057767,
117 0141644,0011417,0036155,0137305,
118 0041176,0076556,0004043,0125430,
119 0140227,0073347,0152776,0067251,
120 };
121 #endif
122 #ifdef IBMPC
123 static unsigned short P0[20] = {
124 0x142d,0x0e5e,0xfb4f,0xc04d,
125 0xedd9,0x9ea1,0x8011,0x4058,
126 0xdbba,0x8806,0x5690,0xc04c,
127 0xc1fd,0x3bd7,0xdcce,0x402b,
128 0xca7e,0x8ef6,0xd39f,0xbff3,
129 };
130 static unsigned short Q0[36] = {
131 /*0x0000,0x0000,0x0000,0x3ff0,*/
132 0x74f9,0xd220,0x4595,0x3fff,
133 0xe5b6,0x8417,0xb482,0x4012,
134 0x81d4,0x350b,0x970e,0x4055,
135 0x365f,0x56c2,0x2ece,0xc06c,
136 0xcbff,0xa8e9,0x0853,0x4069,
137 0xb7d9,0xe78d,0x8261,0xc054,
138 0x7563,0xc104,0xcfad,0x402f,
139 0xcdd5,0xfabf,0xeedc,0xbff2,
140 };
141 #endif
142 #ifdef MIEEE
143 static unsigned short P0[20] = {
144 0xc04d,0xfb4f,0x0e5e,0x142d,
145 0x4058,0x8011,0x9ea1,0xedd9,
146 0xc04c,0x5690,0x8806,0xdbba,
147 0x402b,0xdcce,0x3bd7,0xc1fd,
148 0xbff3,0xd39f,0x8ef6,0xca7e,
149 };
150 static unsigned short Q0[32] = {
151 /*0x3ff0,0x0000,0x0000,0x0000,*/
152 0x3fff,0x4595,0xd220,0x74f9,
153 0x4012,0xb482,0x8417,0xe5b6,
154 0x4055,0x970e,0x350b,0x81d4,
155 0xc06c,0x2ece,0x56c2,0x365f,
156 0x4069,0x0853,0xa8e9,0xcbff,
157 0xc054,0x8261,0xe78d,0xb7d9,
158 0x402f,0xcfad,0xc104,0x7563,
159 0xbff2,0xeedc,0xfabf,0xcdd5,
160 };
161 #endif
162
163
164 /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
165 * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
166 */
167 #ifdef UNK
168 static double P1[9] = {
169 4.05544892305962419923E0,
170 3.15251094599893866154E1,
171 5.71628192246421288162E1,
172 4.40805073893200834700E1,
173 1.46849561928858024014E1,
174 2.18663306850790267539E0,
175 -1.40256079171354495875E-1,
176 -3.50424626827848203418E-2,
177 -8.57456785154685413611E-4,
178 };
179 static double Q1[8] = {
180 /* 1.00000000000000000000E0,*/
181 1.57799883256466749731E1,
182 4.53907635128879210584E1,
183 4.13172038254672030440E1,
184 1.50425385692907503408E1,
185 2.50464946208309415979E0,
186 -1.42182922854787788574E-1,
187 -3.80806407691578277194E-2,
188 -9.33259480895457427372E-4,
189 };
190 #endif
191 #ifdef DEC
192 static unsigned short P1[36] = {
193 0040601,0143074,0150744,0073326,
194 0041374,0031554,0113253,0146016,
195 0041544,0123272,0012463,0176771,
196 0041460,0051160,0103560,0156511,
197 0041152,0172624,0117772,0030755,
198 0040413,0170713,0151545,0176413,
199 0137417,0117512,0022154,0131671,
200 0137017,0104257,0071432,0007072,
201 0135540,0143363,0063137,0036166,
202 };
203 static unsigned short Q1[32] = {
204 /*0040200,0000000,0000000,0000000,*/
205 0041174,0075325,0004736,0120326,
206 0041465,0110044,0047561,0045567,
207 0041445,0042321,0012142,0030340,
208 0041160,0127074,0166076,0141051,
209 0040440,0046055,0040745,0150400,
210 0137421,0114146,0067330,0010621,
211 0137033,0175162,0025555,0114351,
212 0135564,0122773,0145750,0030357,
213 };
214 #endif
215 #ifdef IBMPC
216 static unsigned short P1[36] = {
217 0x8edb,0x9a3c,0x38c7,0x4010,
218 0x7982,0x92d5,0x866d,0x403f,
219 0x7fbf,0x42a6,0x94d7,0x404c,
220 0x1ba9,0x10ee,0x0a4e,0x4046,
221 0x463e,0x93ff,0x5eb2,0x402d,
222 0xbfa1,0x7a6c,0x7e39,0x4001,
223 0x9677,0x448d,0xf3e9,0xbfc1,
224 0x41c7,0xee63,0xf115,0xbfa1,
225 0xe78f,0x6ccb,0x18de,0xbf4c,
226 };
227 static unsigned short Q1[32] = {
228 /*0x0000,0x0000,0x0000,0x3ff0,*/
229 0xd41b,0xa13b,0x8f5a,0x402f,
230 0x296f,0x89ee,0xb204,0x4046,
231 0x461c,0x228c,0xa89a,0x4044,
232 0xd845,0x9d87,0x15c7,0x402e,
233 0xba20,0xa83c,0x0985,0x4004,
234 0x0232,0xcddb,0x330c,0xbfc2,
235 0xb31d,0x456d,0x7f4e,0xbfa3,
236 0x061e,0x797d,0x94bf,0xbf4e,
237 };
238 #endif
239 #ifdef MIEEE
240 static unsigned short P1[36] = {
241 0x4010,0x38c7,0x9a3c,0x8edb,
242 0x403f,0x866d,0x92d5,0x7982,
243 0x404c,0x94d7,0x42a6,0x7fbf,
244 0x4046,0x0a4e,0x10ee,0x1ba9,
245 0x402d,0x5eb2,0x93ff,0x463e,
246 0x4001,0x7e39,0x7a6c,0xbfa1,
247 0xbfc1,0xf3e9,0x448d,0x9677,
248 0xbfa1,0xf115,0xee63,0x41c7,
249 0xbf4c,0x18de,0x6ccb,0xe78f,
250 };
251 static unsigned short Q1[32] = {
252 /*0x3ff0,0x0000,0x0000,0x0000,*/
253 0x402f,0x8f5a,0xa13b,0xd41b,
254 0x4046,0xb204,0x89ee,0x296f,
255 0x4044,0xa89a,0x228c,0x461c,
256 0x402e,0x15c7,0x9d87,0xd845,
257 0x4004,0x0985,0xa83c,0xba20,
258 0xbfc2,0x330c,0xcddb,0x0232,
259 0xbfa3,0x7f4e,0x456d,0xb31d,
260 0xbf4e,0x94bf,0x797d,0x061e,
261 };
262 #endif
263
264 /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
265 * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
266 */
267
268 #ifdef UNK
269 static double P2[9] = {
270 3.23774891776946035970E0,
271 6.91522889068984211695E0,
272 3.93881025292474443415E0,
273 1.33303460815807542389E0,
274 2.01485389549179081538E-1,
275 1.23716634817820021358E-2,
276 3.01581553508235416007E-4,
277 2.65806974686737550832E-6,
278 6.23974539184983293730E-9,
279 };
280 static double Q2[8] = {
281 /* 1.00000000000000000000E0,*/
282 6.02427039364742014255E0,
283 3.67983563856160859403E0,
284 1.37702099489081330271E0,
285 2.16236993594496635890E-1,
286 1.34204006088543189037E-2,
287 3.28014464682127739104E-4,
288 2.89247864745380683936E-6,
289 6.79019408009981274425E-9,
290 };
291 #endif
292 #ifdef DEC
293 static unsigned short P2[36] = {
294 0040517,0033507,0036236,0125641,
295 0040735,0044616,0014473,0140133,
296 0040574,0012567,0114535,0102541,
297 0040252,0120340,0143474,0150135,
298 0037516,0051057,0115361,0031211,
299 0036512,0131204,0101511,0125144,
300 0035236,0016627,0043160,0140216,
301 0033462,0060512,0060141,0010641,
302 0031326,0062541,0101304,0077706,
303 };
304 static unsigned short Q2[32] = {
305 /*0040200,0000000,0000000,0000000,*/
306 0040700,0143322,0132137,0040501,
307 0040553,0101155,0053221,0140257,
308 0040260,0041071,0052573,0010004,
309 0037535,0066472,0177261,0162330,
310 0036533,0160475,0066666,0036132,
311 0035253,0174533,0027771,0044027,
312 0033502,0016147,0117666,0063671,
313 0031351,0047455,0141663,0054751,
314 };
315 #endif
316 #ifdef IBMPC
317 static unsigned short P2[36] = {
318 0xd574,0xe793,0xe6e8,0x4009,
319 0x780b,0xc327,0xa931,0x401b,
320 0xb0ac,0xf32b,0x82ae,0x400f,
321 0x9a0c,0x18e7,0x541c,0x3ff5,
322 0x2651,0xf35e,0xca45,0x3fc9,
323 0x354d,0x9069,0x5650,0x3f89,
324 0x1812,0xe8ce,0xc3b2,0x3f33,
325 0x2234,0x4c0c,0x4c29,0x3ec6,
326 0x8ff9,0x3058,0xccac,0x3e3a,
327 };
328 static unsigned short Q2[32] = {
329 /*0x0000,0x0000,0x0000,0x3ff0,*/
330 0xe828,0x568b,0x18da,0x4018,
331 0x3816,0xaad2,0x704d,0x400d,
332 0x6200,0x2aaf,0x0847,0x3ff6,
333 0x3c9b,0x5fd6,0xada7,0x3fcb,
334 0xc78b,0xadb6,0x7c27,0x3f8b,
335 0x2903,0x65ff,0x7f2b,0x3f35,
336 0xccf7,0xf3f6,0x438c,0x3ec8,
337 0x6b3d,0xb876,0x29e5,0x3e3d,
338 };
339 #endif
340 #ifdef MIEEE
341 static unsigned short P2[36] = {
342 0x4009,0xe6e8,0xe793,0xd574,
343 0x401b,0xa931,0xc327,0x780b,
344 0x400f,0x82ae,0xf32b,0xb0ac,
345 0x3ff5,0x541c,0x18e7,0x9a0c,
346 0x3fc9,0xca45,0xf35e,0x2651,
347 0x3f89,0x5650,0x9069,0x354d,
348 0x3f33,0xc3b2,0xe8ce,0x1812,
349 0x3ec6,0x4c29,0x4c0c,0x2234,
350 0x3e3a,0xccac,0x3058,0x8ff9,
351 };
352 static unsigned short Q2[32] = {
353 /*0x3ff0,0x0000,0x0000,0x0000,*/
354 0x4018,0x18da,0x568b,0xe828,
355 0x400d,0x704d,0xaad2,0x3816,
356 0x3ff6,0x0847,0x2aaf,0x6200,
357 0x3fcb,0xada7,0x5fd6,0x3c9b,
358 0x3f8b,0x7c27,0xadb6,0xc78b,
359 0x3f35,0x7f2b,0x65ff,0x2903,
360 0x3ec8,0x438c,0xf3f6,0xccf7,
361 0x3e3d,0x29e5,0xb876,0x6b3d,
362 };
363 #endif
364
ndtri(y0)365 double ndtri(y0)
366 double y0;
367 {
368 double x, y, z, y2, x0, x1;
369 int code;
370
371 if( y0 <= 0.0 )
372 {
373 mtherr( "ndtri", DOMAIN );
374 return( -MAXNUM );
375 }
376 if( y0 >= 1.0 )
377 {
378 mtherr( "ndtri", DOMAIN );
379 return( MAXNUM );
380 }
381 code = 1;
382 y = y0;
383 if( y > (1.0 - 0.13533528323661269189) ) /* 0.135... = exp(-2) */
384 {
385 y = 1.0 - y;
386 code = 0;
387 }
388
389 if( y > 0.13533528323661269189 )
390 {
391 y = y - 0.5;
392 y2 = y * y;
393 x = y + y * (y2 * polevl( y2, P0, 4)/p1evl( y2, Q0, 8 ));
394 x = x * s2pi;
395 return(x);
396 }
397
398 x = sqrt( -2.0 * log(y) );
399 x0 = x - log(x)/x;
400
401 z = 1.0/x;
402 if( x < 8.0 ) /* y > exp(-32) = 1.2664165549e-14 */
403 x1 = z * polevl( z, P1, 8 )/p1evl( z, Q1, 8 );
404 else
405 x1 = z * polevl( z, P2, 8 )/p1evl( z, Q2, 8 );
406 x = x0 - x1;
407 if( code != 0 )
408 x = -x;
409 return( x );
410 }
411