source: sasmodels/sasmodels/models/lib/sas_erf.c @ 1557a1e

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 1557a1e was 1557a1e, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

windows MSVC runtime is missing erf

  • Property mode set to 100644
File size: 6.3 KB
Line 
1
2/*
3 * Cephes Math Library Release 2.2:  June, 1992
4 * Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
5 * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
6 */
7/*
8 *
9 *      Error function
10 *
11 *
12 *
13 * SYNOPSIS:
14 *
15 * double x, y, erf();
16 *
17 * y = erf( x );
18 *
19 *
20 *
21 * DESCRIPTION:
22 *
23 * The integral is
24 *
25 *                           x
26 *                            -
27 *                 2         | |          2
28 *   erf(x)  =  --------     |    exp( - t  ) dt.
29 *              sqrt(pi)   | |
30 *                          -
31 *                           0
32 *
33 * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
34 * erf(x) = 1 - erfc(x).
35 *
36 *
37 *
38 * ACCURACY:
39 *
40 *                      Relative error:
41 * arithmetic   domain     # trials      peak         rms
42 *    IEEE      0,1         30000       3.7e-16     1.0e-16
43 *
44 */
45/*
46 *
47 *      Complementary error function
48 *
49 *
50 *
51 * SYNOPSIS:
52 *
53 * double x, y, erfc();
54 *
55 * y = erfc( x );
56 *
57 *
58 *
59 * DESCRIPTION:
60 *
61 *
62 *  1 - erf(x) =
63 *
64 *                           inf.
65 *                             -
66 *                  2         | |          2
67 *   erfc(x)  =  --------     |    exp( - t  ) dt
68 *               sqrt(pi)   | |
69 *                           -
70 *                            x
71 *
72 *
73 * For small x, erfc(x) = 1 - erf(x); otherwise rational
74 * approximations are computed.
75 *
76 *
77 *
78 * ACCURACY:
79 *
80 *                      Relative error:
81 * arithmetic   domain     # trials      peak         rms
82 *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
83 */
84
85#ifdef NEED_ERF
86
87#if FLOAT_SIZE>4  // DOUBLE_PRECISION
88double erf(double x);
89double erfc(double a);
90
91constant double PD[] = {
92    2.46196981473530512524E-10,
93    5.64189564831068821977E-1,
94    7.46321056442269912687E0,
95    4.86371970985681366614E1,
96    1.96520832956077098242E2,
97    5.26445194995477358631E2,
98    9.34528527171957607540E2,
99    1.02755188689515710272E3,
100    5.57535335369399327526E2
101};
102
103constant double QD[] = {
104    /* 1.00000000000000000000E0, */
105    1.32281951154744992508E1,
106    8.67072140885989742329E1,
107    3.54937778887819891062E2,
108    9.75708501743205489753E2,
109    1.82390916687909736289E3,
110    2.24633760818710981792E3,
111    1.65666309194161350182E3,
112    5.57535340817727675546E2
113};
114
115constant double RD[] = {
116    5.64189583547755073984E-1,
117    1.27536670759978104416E0,
118    5.01905042251180477414E0,
119    6.16021097993053585195E0,
120    7.40974269950448939160E0,
121    2.97886665372100240670E0
122};
123
124constant double SD[] = {
125    /* 1.00000000000000000000E0, */
126    2.26052863220117276590E0,
127    9.39603524938001434673E0,
128    1.20489539808096656605E1,
129    1.70814450747565897222E1,
130    9.60896809063285878198E0,
131    3.36907645100081516050E0
132};
133
134constant double TD[] = {
135    9.60497373987051638749E0,
136    9.00260197203842689217E1,
137    2.23200534594684319226E3,
138    7.00332514112805075473E3,
139    5.55923013010394962768E4
140};
141
142constant double UD[] = {
143    /* 1.00000000000000000000E0, */
144    3.35617141647503099647E1,
145    5.21357949780152679795E2,
146    4.59432382970980127987E3,
147    2.26290000613890934246E4,
148    4.92673942608635921086E4
149};
150
151double erfc(double a)
152{
153    double MAXLOG = 88.72283905206835;
154    double p, q, x, y, z;
155
156
157    /*if (a < 0.0)
158        x = -a;
159    else
160        x = a;*/
161
162    x = fabs(a);
163
164
165    if (x < 1.0) {
166        //The line bellow is a troublemaker for GPU, so sas_erf function
167        //is explicit here for the case < 1.0
168        //return (1.0 - sas_erf(a));
169        z = x * x;
170        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
171
172        return y;
173    }
174
175    z = -a * a;
176
177    if (z < -MAXLOG) {
178        if (a < 0)
179            return (2.0);
180        else
181            return (0.0);
182    }
183
184    z = exp(z);
185
186
187    if (x < 8.0) {
188        p = polevl(x, PD, 8);
189        q = p1evl(x, QD, 8);
190    }
191    else {
192        p = polevl(x, RD, 5);
193        q = p1evl(x, SD, 6);
194    }
195    y = (z * p) / q;
196
197    if (a < 0)
198        y = 2.0 - y;
199
200    if (y == 0.0) {
201        if (a < 0)
202            return (2.0);
203        else
204            return (0.0);
205    }
206    return y;
207}
208
209
210double erf(double x)
211{
212    double y, z;
213
214    if (fabs(x) > 1.0)
215        return (1.0 - erfc(x));
216
217    z = x * x;
218    #if FLOAT_SIZE>4
219        y = x * polevl(z, TD, 4) / p1evl(z, UD, 5);
220    #else
221        y = x * polevl( z, TF, 6 );
222    #endif
223
224    return y;
225}
226
227#else // SINGLE PRECISION
228
229double erff(double x);
230double erfcf(double a);
231
232/* erfc(x) = exp(-x^2) P(1/x), 1 < x < 2 */
233constant double PF[] = {
234    2.326819970068386E-002,
235    -1.387039388740657E-001,
236    3.687424674597105E-001,
237    -5.824733027278666E-001,
238    6.210004621745983E-001,
239    -4.944515323274145E-001,
240    3.404879937665872E-001,
241    -2.741127028184656E-001,
242    5.638259427386472E-001
243};
244
245/* erfc(x) = exp(-x^2) 1/x P(1/x^2), 2 < x < 14 */
246constant double RF[] = {
247    -1.047766399936249E+001,
248    1.297719955372516E+001,
249    -7.495518717768503E+000,
250    2.921019019210786E+000,
251    -1.015265279202700E+000,
252    4.218463358204948E-001,
253    -2.820767439740514E-001,
254    5.641895067754075E-001
255};
256
257/* erf(x) = x P(x^2), 0 < x < 1 */
258 constant double TF[] = {
259    7.853861353153693E-005,
260    -8.010193625184903E-004,
261    5.188327685732524E-003,
262    -2.685381193529856E-002,
263    1.128358514861418E-001,
264    -3.761262582423300E-001,
265    1.128379165726710E+000
266};
267
268
269float erfcf(float a)
270{
271    float MAXLOG = 88.72283905206835;
272    float p, q, x, y, z;
273
274
275    /*if (a < 0.0)
276        x = -a;
277    else
278        x = a;*/
279
280    x = fabsf(a);
281
282
283    if (x < 1.0) {
284        //The line below is a troublemaker for GPU, so sas_erf function
285        //is explicit here for the case < 1.0
286        //return (1.0 - sas_erf(a));
287        z = x * x;
288        y = x * polevl( z, TF, 6 );
289
290        return y;
291    }
292
293    z = -a * a;
294
295    if (z < -MAXLOG) {
296        if (a < 0)
297            return (2.0);
298        else
299            return (0.0);
300    }
301
302    z = expf(z);
303
304
305    q=1.0/x;
306    y=q*q;
307    if( x < 2.0 ) {
308        p = polevl( y, PF, 8 );
309    } else {
310        p = polevl( y, RF, 7 );
311    }
312    y = z * q * p;
313
314    if (a < 0)
315        y = 2.0 - y;
316
317    if (y == 0.0) {
318        if (a < 0)
319            return (2.0);
320        else
321            return (0.0);
322    }
323    return y;
324}
325
326
327float erff(float x)
328{
329    float y, z;
330
331    if (fabsf(x) > 1.0)
332        return (1.0 - erfcf(x));
333
334    z = x * x;
335    y = x * polevl( z, TF, 6 );
336
337    return y;
338}
339
340#endif // SINGLE_PRECISION
341#endif // NEED_ERF
342
343#if FLOAT_SIZE>4
344#define sas_erf erf
345#define sas_erfc erfc
346#else
347#define sas_erf erff
348#define sas_erfc erfcf
349#endif
Note: See TracBrowser for help on using the repository browser.