← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cephes/ndtri.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

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