← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/cephes/ndtr.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/ndtr.c
 
1
/*                                                      ndtr.c
 
2
 *
 
3
 *      Normal distribution function
 
4
 *
 
5
 *
 
6
 *
 
7
 * SYNOPSIS:
 
8
 *
 
9
 * double x, y, ndtr();
 
10
 *
 
11
 * y = ndtr( x );
 
12
 *
 
13
 *
 
14
 *
 
15
 * DESCRIPTION:
 
16
 *
 
17
 * Returns the area under the Gaussian probability density
 
18
 * function, integrated from minus infinity to x:
 
19
 *
 
20
 *                            x
 
21
 *                             -
 
22
 *                   1        | |          2
 
23
 *    ndtr(x)  = ---------    |    exp( - t /2 ) dt
 
24
 *               sqrt(2pi)  | |
 
25
 *                           -
 
26
 *                          -inf.
 
27
 *
 
28
 *             =  ( 1 + erf(z) ) / 2
 
29
 *             =  erfc(z) / 2
 
30
 *
 
31
 * where z = x/sqrt(2). Computation is via the functions
 
32
 * erf and erfc.
 
33
 *
 
34
 *
 
35
 * ACCURACY:
 
36
 *
 
37
 *                      Relative error:
 
38
 * arithmetic   domain     # trials      peak         rms
 
39
 *    DEC      -13,0         8000       2.1e-15     4.8e-16
 
40
 *    IEEE     -13,0        30000       3.4e-14     6.7e-15
 
41
 *
 
42
 *
 
43
 * ERROR MESSAGES:
 
44
 *
 
45
 *   message         condition         value returned
 
46
 * erfc underflow    x > 37.519379347       0.0
 
47
 *
 
48
 */
 
49
/*                                                     erf.c
 
50
 *
 
51
 *      Error function
 
52
 *
 
53
 *
 
54
 *
 
55
 * SYNOPSIS:
 
56
 *
 
57
 * double x, y, erf();
 
58
 *
 
59
 * y = erf( x );
 
60
 *
 
61
 *
 
62
 *
 
63
 * DESCRIPTION:
 
64
 *
 
65
 * The integral is
 
66
 *
 
67
 *                           x 
 
68
 *                            -
 
69
 *                 2         | |          2
 
70
 *   erf(x)  =  --------     |    exp( - t  ) dt.
 
71
 *              sqrt(pi)   | |
 
72
 *                          -
 
73
 *                           0
 
74
 *
 
75
 * The magnitude of x is limited to 9.231948545 for DEC
 
76
 * arithmetic; 1 or -1 is returned outside this range.
 
77
 *
 
78
 * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
 
79
 * erf(x) = 1 - erfc(x).
 
80
 *
 
81
 *
 
82
 *
 
83
 * ACCURACY:
 
84
 *
 
85
 *                      Relative error:
 
86
 * arithmetic   domain     # trials      peak         rms
 
87
 *    DEC       0,1         14000       4.7e-17     1.5e-17
 
88
 *    IEEE      0,1         30000       3.7e-16     1.0e-16
 
89
 *
 
90
 */
 
91
/*                                                     erfc.c
 
92
 *
 
93
 *      Complementary error function
 
94
 *
 
95
 *
 
96
 *
 
97
 * SYNOPSIS:
 
98
 *
 
99
 * double x, y, erfc();
 
100
 *
 
101
 * y = erfc( x );
 
102
 *
 
103
 *
 
104
 *
 
105
 * DESCRIPTION:
 
106
 *
 
107
 *
 
108
 *  1 - erf(x) =
 
109
 *
 
110
 *                           inf. 
 
111
 *                             -
 
112
 *                  2         | |          2
 
113
 *   erfc(x)  =  --------     |    exp( - t  ) dt
 
114
 *               sqrt(pi)   | |
 
115
 *                           -
 
116
 *                            x
 
117
 *
 
118
 *
 
119
 * For small x, erfc(x) = 1 - erf(x); otherwise rational
 
120
 * approximations are computed.
 
121
 *
 
122
 *
 
123
 *
 
124
 * ACCURACY:
 
125
 *
 
126
 *                      Relative error:
 
127
 * arithmetic   domain     # trials      peak         rms
 
128
 *    DEC       0, 9.2319   12000       5.1e-16     1.2e-16
 
129
 *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
 
130
 *
 
131
 *
 
132
 * ERROR MESSAGES:
 
133
 *
 
134
 *   message         condition              value returned
 
135
 * erfc underflow    x > 9.231948545 (DEC)       0.0
 
136
 *
 
137
 *
 
138
 */
 
139
 
 
140
 
 
141
/*
 
142
Cephes Math Library Release 2.2:  June, 1992
 
143
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
 
144
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
145
*/
 
146
 
 
147
 
 
148
#include "mconf.h"
 
149
 
 
150
extern double SQRTH, NAN;
 
151
extern double MAXLOG;
 
152
 
 
153
#ifdef UNK
 
154
static double P[] = {
 
155
 2.46196981473530512524E-10,
 
156
 5.64189564831068821977E-1,
 
157
 7.46321056442269912687E0,
 
158
 4.86371970985681366614E1,
 
159
 1.96520832956077098242E2,
 
160
 5.26445194995477358631E2,
 
161
 9.34528527171957607540E2,
 
162
 1.02755188689515710272E3,
 
163
 5.57535335369399327526E2
 
164
};
 
165
static double Q[] = {
 
166
/* 1.00000000000000000000E0,*/
 
167
 1.32281951154744992508E1,
 
168
 8.67072140885989742329E1,
 
169
 3.54937778887819891062E2,
 
170
 9.75708501743205489753E2,
 
171
 1.82390916687909736289E3,
 
172
 2.24633760818710981792E3,
 
173
 1.65666309194161350182E3,
 
174
 5.57535340817727675546E2
 
175
};
 
176
static double R[] = {
 
177
 5.64189583547755073984E-1,
 
178
 1.27536670759978104416E0,
 
179
 5.01905042251180477414E0,
 
180
 6.16021097993053585195E0,
 
181
 7.40974269950448939160E0,
 
182
 2.97886665372100240670E0
 
183
};
 
184
static double S[] = {
 
185
/* 1.00000000000000000000E0,*/
 
186
 2.26052863220117276590E0,
 
187
 9.39603524938001434673E0,
 
188
 1.20489539808096656605E1,
 
189
 1.70814450747565897222E1,
 
190
 9.60896809063285878198E0,
 
191
 3.36907645100081516050E0
 
192
};
 
193
static double T[] = {
 
194
 9.60497373987051638749E0,
 
195
 9.00260197203842689217E1,
 
196
 2.23200534594684319226E3,
 
197
 7.00332514112805075473E3,
 
198
 5.55923013010394962768E4
 
199
};
 
200
static double U[] = {
 
201
/* 1.00000000000000000000E0,*/
 
202
 3.35617141647503099647E1,
 
203
 5.21357949780152679795E2,
 
204
 4.59432382970980127987E3,
 
205
 2.26290000613890934246E4,
 
206
 4.92673942608635921086E4
 
207
};
 
208
 
 
209
#define UTHRESH 37.519379347
 
210
#endif
 
211
 
 
212
#ifdef DEC
 
213
static unsigned short P[] = {
 
214
0030207,0054445,0011173,0021706,
 
215
0040020,0067272,0030661,0122075,
 
216
0040756,0151236,0173053,0067042,
 
217
0041502,0106175,0062555,0151457,
 
218
0042104,0102525,0047401,0003667,
 
219
0042403,0116176,0011446,0075303,
 
220
0042551,0120723,0061641,0123275,
 
221
0042600,0070651,0007264,0134516,
 
222
0042413,0061102,0167507,0176625
 
223
};
 
224
static unsigned short Q[] = {
 
225
/*0040200,0000000,0000000,0000000,*/
 
226
0041123,0123257,0165741,0017142,
 
227
0041655,0065027,0173413,0115450,
 
228
0042261,0074011,0021573,0004150,
 
229
0042563,0166530,0013662,0007200,
 
230
0042743,0176427,0162443,0105214,
 
231
0043014,0062546,0153727,0123772,
 
232
0042717,0012470,0006227,0067424,
 
233
0042413,0061103,0003042,0013254
 
234
};
 
235
static unsigned short R[] = {
 
236
0040020,0067272,0101024,0155421,
 
237
0040243,0037467,0056706,0026462,
 
238
0040640,0116017,0120665,0034315,
 
239
0040705,0020162,0143350,0060137,
 
240
0040755,0016234,0134304,0130157,
 
241
0040476,0122700,0051070,0015473
 
242
};
 
243
static unsigned short S[] = {
 
244
/*0040200,0000000,0000000,0000000,*/
 
245
0040420,0126200,0044276,0070413,
 
246
0041026,0053051,0007302,0063746,
 
247
0041100,0144203,0174051,0061151,
 
248
0041210,0123314,0126343,0177646,
 
249
0041031,0137125,0051431,0033011,
 
250
0040527,0117362,0152661,0066201
 
251
};
 
252
static unsigned short T[] = {
 
253
0041031,0126770,0170672,0166101,
 
254
0041664,0006522,0072360,0031770,
 
255
0043013,0100025,0162641,0126671,
 
256
0043332,0155231,0161627,0076200,
 
257
0044131,0024115,0021020,0117343
 
258
};
 
259
static unsigned short U[] = {
 
260
/*0040200,0000000,0000000,0000000,*/
 
261
0041406,0037461,0177575,0032714,
 
262
0042402,0053350,0123061,0153557,
 
263
0043217,0111227,0032007,0164217,
 
264
0043660,0145000,0004013,0160114,
 
265
0044100,0071544,0167107,0125471
 
266
};
 
267
#define UTHRESH 14.0
 
268
#endif
 
269
 
 
270
#ifdef IBMPC
 
271
static unsigned short P[] = {
 
272
0x6479,0xa24f,0xeb24,0x3df0,
 
273
0x3488,0x4636,0x0dd7,0x3fe2,
 
274
0x6dc4,0xdec5,0xda53,0x401d,
 
275
0xba66,0xacad,0x518f,0x4048,
 
276
0x20f7,0xa9e0,0x90aa,0x4068,
 
277
0xcf58,0xc264,0x738f,0x4080,
 
278
0x34d8,0x6c74,0x343a,0x408d,
 
279
0x972a,0x21d6,0x0e35,0x4090,
 
280
0xffb3,0x5de8,0x6c48,0x4081
 
281
};
 
282
static unsigned short Q[] = {
 
283
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
284
0x23cc,0xfd7c,0x74d5,0x402a,
 
285
0x7365,0xfee1,0xad42,0x4055,
 
286
0x610d,0x246f,0x2f01,0x4076,
 
287
0x41d0,0x02f6,0x7dab,0x408e,
 
288
0x7151,0xfca4,0x7fa2,0x409c,
 
289
0xf4ff,0xdafa,0x8cac,0x40a1,
 
290
0xede2,0x0192,0xe2a7,0x4099,
 
291
0x42d6,0x60c4,0x6c48,0x4081
 
292
};
 
293
static unsigned short R[] = {
 
294
0x9b62,0x5042,0x0dd7,0x3fe2,
 
295
0xc5a6,0xebb8,0x67e6,0x3ff4,
 
296
0xa71a,0xf436,0x1381,0x4014,
 
297
0x0c0c,0x58dd,0xa40e,0x4018,
 
298
0x960e,0x9718,0xa393,0x401d,
 
299
0x0367,0x0a47,0xd4b8,0x4007
 
300
};
 
301
static unsigned short S[] = {
 
302
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
303
0xce21,0x0917,0x1590,0x4002,
 
304
0x4cfd,0x21d8,0xcac5,0x4022,
 
305
0x2c4d,0x7f05,0x1910,0x4028,
 
306
0x7ff5,0x959c,0x14d9,0x4031,
 
307
0x26c1,0xaa63,0x37ca,0x4023,
 
308
0x2d90,0x5ab6,0xf3de,0x400a
 
309
};
 
310
static unsigned short T[] = {
 
311
0x5d88,0x1e37,0x35bf,0x4023,
 
312
0x067f,0x4e9e,0x81aa,0x4056,
 
313
0x35b7,0xbcb4,0x7002,0x40a1,
 
314
0xef90,0x3c72,0x5b53,0x40bb,
 
315
0x13dc,0xa442,0x2509,0x40eb
 
316
};
 
317
static unsigned short U[] = {
 
318
/*0x0000,0x0000,0x0000,0x3ff0,*/
 
319
0xa6ba,0x3fef,0xc7e6,0x4040,
 
320
0x3aee,0x14c6,0x4add,0x4080,
 
321
0xfd12,0xe680,0xf252,0x40b1,
 
322
0x7c0a,0x0101,0x1940,0x40d6,
 
323
0xf567,0x9dc8,0x0e6c,0x40e8
 
324
};
 
325
#define UTHRESH 37.519379347
 
326
#endif
 
327
 
 
328
#ifdef MIEEE
 
329
static unsigned short P[] = {
 
330
0x3df0,0xeb24,0xa24f,0x6479,
 
331
0x3fe2,0x0dd7,0x4636,0x3488,
 
332
0x401d,0xda53,0xdec5,0x6dc4,
 
333
0x4048,0x518f,0xacad,0xba66,
 
334
0x4068,0x90aa,0xa9e0,0x20f7,
 
335
0x4080,0x738f,0xc264,0xcf58,
 
336
0x408d,0x343a,0x6c74,0x34d8,
 
337
0x4090,0x0e35,0x21d6,0x972a,
 
338
0x4081,0x6c48,0x5de8,0xffb3
 
339
};
 
340
static unsigned short Q[] = {
 
341
0x402a,0x74d5,0xfd7c,0x23cc,
 
342
0x4055,0xad42,0xfee1,0x7365,
 
343
0x4076,0x2f01,0x246f,0x610d,
 
344
0x408e,0x7dab,0x02f6,0x41d0,
 
345
0x409c,0x7fa2,0xfca4,0x7151,
 
346
0x40a1,0x8cac,0xdafa,0xf4ff,
 
347
0x4099,0xe2a7,0x0192,0xede2,
 
348
0x4081,0x6c48,0x60c4,0x42d6
 
349
};
 
350
static unsigned short R[] = {
 
351
0x3fe2,0x0dd7,0x5042,0x9b62,
 
352
0x3ff4,0x67e6,0xebb8,0xc5a6,
 
353
0x4014,0x1381,0xf436,0xa71a,
 
354
0x4018,0xa40e,0x58dd,0x0c0c,
 
355
0x401d,0xa393,0x9718,0x960e,
 
356
0x4007,0xd4b8,0x0a47,0x0367
 
357
};
 
358
static unsigned short S[] = {
 
359
0x4002,0x1590,0x0917,0xce21,
 
360
0x4022,0xcac5,0x21d8,0x4cfd,
 
361
0x4028,0x1910,0x7f05,0x2c4d,
 
362
0x4031,0x14d9,0x959c,0x7ff5,
 
363
0x4023,0x37ca,0xaa63,0x26c1,
 
364
0x400a,0xf3de,0x5ab6,0x2d90
 
365
};
 
366
static unsigned short T[] = {
 
367
0x4023,0x35bf,0x1e37,0x5d88,
 
368
0x4056,0x81aa,0x4e9e,0x067f,
 
369
0x40a1,0x7002,0xbcb4,0x35b7,
 
370
0x40bb,0x5b53,0x3c72,0xef90,
 
371
0x40eb,0x2509,0xa442,0x13dc
 
372
};
 
373
static unsigned short U[] = {
 
374
0x4040,0xc7e6,0x3fef,0xa6ba,
 
375
0x4080,0x4add,0x14c6,0x3aee,
 
376
0x40b1,0xf252,0xe680,0xfd12,
 
377
0x40d6,0x1940,0x0101,0x7c0a,
 
378
0x40e8,0x0e6c,0x9dc8,0xf567
 
379
};
 
380
#define UTHRESH 37.519379347
 
381
#endif
 
382
 
 
383
#ifndef ANSIPROT
 
384
double polevl(), p1evl(), exp(), log(), fabs();
 
385
double erf(), erfc();
 
386
#endif
 
387
 
 
388
double ndtr(double a)
 
389
{
 
390
double x, y, z;
 
391
 
 
392
if (isnan(a)) {
 
393
  mtherr("ndtr", DOMAIN);
 
394
  return (NAN);
 
395
}
 
396
 
 
397
x = a * SQRTH;
 
398
z = fabs(x);
 
399
 
 
400
if( z < SQRTH )
 
401
        y = 0.5 + 0.5 * erf(x);
 
402
 
 
403
else
 
404
        {
 
405
        y = 0.5 * erfc(z);
 
406
 
 
407
        if( x > 0 )
 
408
                y = 1.0 - y;
 
409
        }
 
410
 
 
411
return(y);
 
412
}
 
413
 
 
414
 
 
415
double erfc(double a)
 
416
{
 
417
double p,q,x,y,z;
 
418
 
 
419
if (isnan(a)) {
 
420
  mtherr("erfc", DOMAIN);
 
421
  return (NAN);
 
422
}
 
423
 
 
424
if( a < 0.0 )
 
425
        x = -a;
 
426
else
 
427
        x = a;
 
428
 
 
429
if( x < 1.0 )
 
430
        return( 1.0 - erf(a) );
 
431
 
 
432
z = -a * a;
 
433
 
 
434
if( z < -MAXLOG )
 
435
        {
 
436
under:
 
437
        mtherr( "erfc", UNDERFLOW );
 
438
        if( a < 0 )
 
439
                return( 2.0 );
 
440
        else
 
441
                return( 0.0 );
 
442
        }
 
443
 
 
444
z = exp(z);
 
445
 
 
446
if( x < 8.0 )
 
447
        {
 
448
        p = polevl( x, P, 8 );
 
449
        q = p1evl( x, Q, 8 );
 
450
        }
 
451
else
 
452
        {
 
453
        p = polevl( x, R, 5 );
 
454
        q = p1evl( x, S, 6 );
 
455
        }
 
456
y = (z * p)/q;
 
457
 
 
458
if( a < 0 )
 
459
        y = 2.0 - y;
 
460
 
 
461
if( y == 0.0 )
 
462
        goto under;
 
463
 
 
464
return(y);
 
465
}
 
466
 
 
467
 
 
468
 
 
469
double erf(double x)
 
470
{
 
471
double y, z;
 
472
 
 
473
if (isnan(x)) {
 
474
  mtherr("erf", DOMAIN);
 
475
  return (NAN);
 
476
}
 
477
 
 
478
if( fabs(x) > 1.0 )
 
479
        return( 1.0 - erfc(x) );
 
480
z = x * x;
 
481
 
 
482
y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
 
483
return( y );
 
484
 
 
485
}