source: sasview/src/sas/models/c_extension/cephes/ndtr.c @ 79492222

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 79492222 was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.

  • Property mode set to 100644
File size: 10.7 KB
Line 
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 with care to avoid error amplification in computing exp(-x^2).
33 *
34 *
35 * ACCURACY:
36 *
37 *                      Relative error:
38 * arithmetic   domain     # trials      peak         rms
39 *    IEEE     -13,0        30000       1.3e-15     2.2e-16
40 *
41 *
42 * ERROR MESSAGES:
43 *
44 *   message         condition         value returned
45 * erfc underflow    x > 37.519379347       0.0
46 *
47 */
48/*                                                     erf.c
49 *
50 *      Error function
51 *
52 *
53 *
54 * SYNOPSIS:
55 *
56 * double x, y, erf();
57 *
58 * y = erf( x );
59 *
60 *
61 *
62 * DESCRIPTION:
63 *
64 * The integral is
65 *
66 *                           x
67 *                            -
68 *                 2         | |          2
69 *   erf(x)  =  --------     |    exp( - t  ) dt.
70 *              sqrt(pi)   | |
71 *                          -
72 *                           0
73 *
74 * The magnitude of x is limited to 9.231948545 for DEC
75 * arithmetic; 1 or -1 is returned outside this range.
76 *
77 * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
78 * erf(x) = 1 - erfc(x).
79 *
80 *
81 *
82 * ACCURACY:
83 *
84 *                      Relative error:
85 * arithmetic   domain     # trials      peak         rms
86 *    DEC       0,1         14000       4.7e-17     1.5e-17
87 *    IEEE      0,1         30000       3.7e-16     1.0e-16
88 *
89 */
90/*                                                     erfc.c
91 *
92 *      Complementary error function
93 *
94 *
95 *
96 * SYNOPSIS:
97 *
98 * double x, y, erfc();
99 *
100 * y = erfc( x );
101 *
102 *
103 *
104 * DESCRIPTION:
105 *
106 *
107 *  1 - erf(x) =
108 *
109 *                           inf.
110 *                             -
111 *                  2         | |          2
112 *   erfc(x)  =  --------     |    exp( - t  ) dt
113 *               sqrt(pi)   | |
114 *                           -
115 *                            x
116 *
117 *
118 * For small x, erfc(x) = 1 - erf(x); otherwise rational
119 * approximations are computed.
120 *
121 * A special function expx2.c is used to suppress error amplification
122 * in computing exp(-x^2).
123 *
124 *
125 * ACCURACY:
126 *
127 *                      Relative error:
128 * arithmetic   domain     # trials      peak         rms
129 *    IEEE      0,26.6417   30000       1.3e-15     2.2e-16
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/*
142Cephes Math Library Release 2.9:  November, 2000
143Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier
144*/
145
146
147#include "mconf.h"
148
149extern double SQRTH;
150extern double MAXLOG;
151
152/* Define this macro to suppress error propagation in exp(x^2)
153   by using the expx2 function.  The tradeoff is that doing so
154   generates two calls to the exponential function instead of one.  */
155#define USE_EXPXSQ 1
156
157#ifdef UNK
158static double P[] = {
159 2.46196981473530512524E-10,
160 5.64189564831068821977E-1,
161 7.46321056442269912687E0,
162 4.86371970985681366614E1,
163 1.96520832956077098242E2,
164 5.26445194995477358631E2,
165 9.34528527171957607540E2,
166 1.02755188689515710272E3,
167 5.57535335369399327526E2
168};
169static double Q[] = {
170/* 1.00000000000000000000E0,*/
171 1.32281951154744992508E1,
172 8.67072140885989742329E1,
173 3.54937778887819891062E2,
174 9.75708501743205489753E2,
175 1.82390916687909736289E3,
176 2.24633760818710981792E3,
177 1.65666309194161350182E3,
178 5.57535340817727675546E2
179};
180static double R[] = {
181 5.64189583547755073984E-1,
182 1.27536670759978104416E0,
183 5.01905042251180477414E0,
184 6.16021097993053585195E0,
185 7.40974269950448939160E0,
186 2.97886665372100240670E0
187};
188static double S[] = {
189/* 1.00000000000000000000E0,*/
190 2.26052863220117276590E0,
191 9.39603524938001434673E0,
192 1.20489539808096656605E1,
193 1.70814450747565897222E1,
194 9.60896809063285878198E0,
195 3.36907645100081516050E0
196};
197static double T[] = {
198 9.60497373987051638749E0,
199 9.00260197203842689217E1,
200 2.23200534594684319226E3,
201 7.00332514112805075473E3,
202 5.55923013010394962768E4
203};
204static double U[] = {
205/* 1.00000000000000000000E0,*/
206 3.35617141647503099647E1,
207 5.21357949780152679795E2,
208 4.59432382970980127987E3,
209 2.26290000613890934246E4,
210 4.92673942608635921086E4
211};
212
213#define UTHRESH 37.519379347
214#endif
215
216#ifdef DEC
217static unsigned short P[] = {
2180030207,0054445,0011173,0021706,
2190040020,0067272,0030661,0122075,
2200040756,0151236,0173053,0067042,
2210041502,0106175,0062555,0151457,
2220042104,0102525,0047401,0003667,
2230042403,0116176,0011446,0075303,
2240042551,0120723,0061641,0123275,
2250042600,0070651,0007264,0134516,
2260042413,0061102,0167507,0176625
227};
228static unsigned short Q[] = {
229/*0040200,0000000,0000000,0000000,*/
2300041123,0123257,0165741,0017142,
2310041655,0065027,0173413,0115450,
2320042261,0074011,0021573,0004150,
2330042563,0166530,0013662,0007200,
2340042743,0176427,0162443,0105214,
2350043014,0062546,0153727,0123772,
2360042717,0012470,0006227,0067424,
2370042413,0061103,0003042,0013254
238};
239static unsigned short R[] = {
2400040020,0067272,0101024,0155421,
2410040243,0037467,0056706,0026462,
2420040640,0116017,0120665,0034315,
2430040705,0020162,0143350,0060137,
2440040755,0016234,0134304,0130157,
2450040476,0122700,0051070,0015473
246};
247static unsigned short S[] = {
248/*0040200,0000000,0000000,0000000,*/
2490040420,0126200,0044276,0070413,
2500041026,0053051,0007302,0063746,
2510041100,0144203,0174051,0061151,
2520041210,0123314,0126343,0177646,
2530041031,0137125,0051431,0033011,
2540040527,0117362,0152661,0066201
255};
256static unsigned short T[] = {
2570041031,0126770,0170672,0166101,
2580041664,0006522,0072360,0031770,
2590043013,0100025,0162641,0126671,
2600043332,0155231,0161627,0076200,
2610044131,0024115,0021020,0117343
262};
263static unsigned short U[] = {
264/*0040200,0000000,0000000,0000000,*/
2650041406,0037461,0177575,0032714,
2660042402,0053350,0123061,0153557,
2670043217,0111227,0032007,0164217,
2680043660,0145000,0004013,0160114,
2690044100,0071544,0167107,0125471
270};
271#define UTHRESH 14.0
272#endif
273
274#ifdef IBMPC
275static unsigned short P[] = {
2760x6479,0xa24f,0xeb24,0x3df0,
2770x3488,0x4636,0x0dd7,0x3fe2,
2780x6dc4,0xdec5,0xda53,0x401d,
2790xba66,0xacad,0x518f,0x4048,
2800x20f7,0xa9e0,0x90aa,0x4068,
2810xcf58,0xc264,0x738f,0x4080,
2820x34d8,0x6c74,0x343a,0x408d,
2830x972a,0x21d6,0x0e35,0x4090,
2840xffb3,0x5de8,0x6c48,0x4081
285};
286static unsigned short Q[] = {
287/*0x0000,0x0000,0x0000,0x3ff0,*/
2880x23cc,0xfd7c,0x74d5,0x402a,
2890x7365,0xfee1,0xad42,0x4055,
2900x610d,0x246f,0x2f01,0x4076,
2910x41d0,0x02f6,0x7dab,0x408e,
2920x7151,0xfca4,0x7fa2,0x409c,
2930xf4ff,0xdafa,0x8cac,0x40a1,
2940xede2,0x0192,0xe2a7,0x4099,
2950x42d6,0x60c4,0x6c48,0x4081
296};
297static unsigned short R[] = {
2980x9b62,0x5042,0x0dd7,0x3fe2,
2990xc5a6,0xebb8,0x67e6,0x3ff4,
3000xa71a,0xf436,0x1381,0x4014,
3010x0c0c,0x58dd,0xa40e,0x4018,
3020x960e,0x9718,0xa393,0x401d,
3030x0367,0x0a47,0xd4b8,0x4007
304};
305static unsigned short S[] = {
306/*0x0000,0x0000,0x0000,0x3ff0,*/
3070xce21,0x0917,0x1590,0x4002,
3080x4cfd,0x21d8,0xcac5,0x4022,
3090x2c4d,0x7f05,0x1910,0x4028,
3100x7ff5,0x959c,0x14d9,0x4031,
3110x26c1,0xaa63,0x37ca,0x4023,
3120x2d90,0x5ab6,0xf3de,0x400a
313};
314static unsigned short T[] = {
3150x5d88,0x1e37,0x35bf,0x4023,
3160x067f,0x4e9e,0x81aa,0x4056,
3170x35b7,0xbcb4,0x7002,0x40a1,
3180xef90,0x3c72,0x5b53,0x40bb,
3190x13dc,0xa442,0x2509,0x40eb
320};
321static unsigned short U[] = {
322/*0x0000,0x0000,0x0000,0x3ff0,*/
3230xa6ba,0x3fef,0xc7e6,0x4040,
3240x3aee,0x14c6,0x4add,0x4080,
3250xfd12,0xe680,0xf252,0x40b1,
3260x7c0a,0x0101,0x1940,0x40d6,
3270xf567,0x9dc8,0x0e6c,0x40e8
328};
329#define UTHRESH 37.519379347
330#endif
331
332#ifdef MIEEE
333static unsigned short P[] = {
3340x3df0,0xeb24,0xa24f,0x6479,
3350x3fe2,0x0dd7,0x4636,0x3488,
3360x401d,0xda53,0xdec5,0x6dc4,
3370x4048,0x518f,0xacad,0xba66,
3380x4068,0x90aa,0xa9e0,0x20f7,
3390x4080,0x738f,0xc264,0xcf58,
3400x408d,0x343a,0x6c74,0x34d8,
3410x4090,0x0e35,0x21d6,0x972a,
3420x4081,0x6c48,0x5de8,0xffb3
343};
344static unsigned short Q[] = {
3450x402a,0x74d5,0xfd7c,0x23cc,
3460x4055,0xad42,0xfee1,0x7365,
3470x4076,0x2f01,0x246f,0x610d,
3480x408e,0x7dab,0x02f6,0x41d0,
3490x409c,0x7fa2,0xfca4,0x7151,
3500x40a1,0x8cac,0xdafa,0xf4ff,
3510x4099,0xe2a7,0x0192,0xede2,
3520x4081,0x6c48,0x60c4,0x42d6
353};
354static unsigned short R[] = {
3550x3fe2,0x0dd7,0x5042,0x9b62,
3560x3ff4,0x67e6,0xebb8,0xc5a6,
3570x4014,0x1381,0xf436,0xa71a,
3580x4018,0xa40e,0x58dd,0x0c0c,
3590x401d,0xa393,0x9718,0x960e,
3600x4007,0xd4b8,0x0a47,0x0367
361};
362static unsigned short S[] = {
3630x4002,0x1590,0x0917,0xce21,
3640x4022,0xcac5,0x21d8,0x4cfd,
3650x4028,0x1910,0x7f05,0x2c4d,
3660x4031,0x14d9,0x959c,0x7ff5,
3670x4023,0x37ca,0xaa63,0x26c1,
3680x400a,0xf3de,0x5ab6,0x2d90
369};
370static unsigned short T[] = {
3710x4023,0x35bf,0x1e37,0x5d88,
3720x4056,0x81aa,0x4e9e,0x067f,
3730x40a1,0x7002,0xbcb4,0x35b7,
3740x40bb,0x5b53,0x3c72,0xef90,
3750x40eb,0x2509,0xa442,0x13dc
376};
377static unsigned short U[] = {
3780x4040,0xc7e6,0x3fef,0xa6ba,
3790x4080,0x4add,0x14c6,0x3aee,
3800x40b1,0xf252,0xe680,0xfd12,
3810x40d6,0x1940,0x0101,0x7c0a,
3820x40e8,0x0e6c,0x9dc8,0xf567
383};
384#define UTHRESH 37.519379347
385#endif
386
387#ifdef ANSIPROT
388extern double polevl ( double, void *, int );
389extern double p1evl ( double, void *, int );
390extern double exp ( double );
391extern double log ( double );
392extern double fabs ( double );
393extern double sqrt ( double );
394extern double expx2 ( double, int );
395double erf ( double );
396double erfc ( double );
397static double erfce ( double );
398#else
399double polevl(), p1evl(), exp(), log(), fabs();
400double erf(), erfc(), expx2(), sqrt();
401static double erfce();
402#endif
403
404double ndtr(a)
405double a;
406{
407double x, y, z;
408
409x = a * SQRTH;
410z = fabs(x);
411
412/* if( z < SQRTH ) */
413if( z < 1.0 )
414        y = 0.5 + 0.5 * erf(x);
415
416else
417        {
418#ifdef USE_EXPXSQ
419        /* See below for erfce. */
420        y = 0.5 * erfce(z);
421        /* Multiply by exp(-x^2 / 2)  */
422        z = expx2(a, -1);
423        y = y * sqrt(z);
424#else
425        y = 0.5 * erfc(z);
426#endif
427        if( x > 0 )
428                y = 1.0 - y;
429        }
430
431return(y);
432}
433
434
435double erfc(a)
436double a;
437{
438double p,q,x,y,z;
439
440
441if( a < 0.0 )
442        x = -a;
443else
444        x = a;
445
446if( x < 1.0 )
447        return( 1.0 - erf(a) );
448
449z = -a * a;
450
451if( z < -MAXLOG )
452        {
453under:
454        mtherr( "erfc", UNDERFLOW );
455        if( a < 0 )
456                return( 2.0 );
457        else
458                return( 0.0 );
459        }
460
461#ifdef USE_EXPXSQ
462/* Compute z = exp(z).  */
463z = expx2(a, -1);
464#else
465z = exp(z);
466#endif
467if( x < 8.0 )
468        {
469        p = polevl( x, P, 8 );
470        q = p1evl( x, Q, 8 );
471        }
472else
473        {
474        p = polevl( x, R, 5 );
475        q = p1evl( x, S, 6 );
476        }
477y = (z * p)/q;
478
479if( a < 0 )
480        y = 2.0 - y;
481
482if( y == 0.0 )
483        goto under;
484
485return(y);
486}
487
488
489/* Exponentially scaled erfc function
490   exp(x^2) erfc(x)
491   valid for x > 1.
492   Use with ndtr and expx2.  */
493static double erfce(x)
494double x;
495{
496double p,q;
497
498if( x < 8.0 )
499        {
500        p = polevl( x, P, 8 );
501        q = p1evl( x, Q, 8 );
502        }
503else
504        {
505        p = polevl( x, R, 5 );
506        q = p1evl( x, S, 6 );
507        }
508return (p/q);
509}
510
511
512
513double erf(x)
514double x;
515{
516double y, z;
517
518if( fabs(x) > 1.0 )
519        return( 1.0 - erfc(x) );
520z = x * x;
521y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
522return( y );
523
524}
Note: See TracBrowser for help on using the repository browser.