1 |
greg |
3.1 |
#ifndef lint |
2 |
greg |
3.2 |
static const char RCSid[] = "$Id: erf.c,v 3.1 2003/02/22 02:07:22 greg Exp $"; |
3 |
greg |
3.1 |
#endif |
4 |
|
|
#ifndef lint |
5 |
|
|
static char sccsid[] = "@(#)erf.c 1.1 87/12/21 SMI"; /* from UCB 4.1 12/25/82 */ |
6 |
|
|
#endif |
7 |
|
|
|
8 |
greg |
3.2 |
#include "rtmath.h" |
9 |
|
|
|
10 |
greg |
3.1 |
/* |
11 |
|
|
C program for floating point error function |
12 |
|
|
|
13 |
|
|
erf(x) returns the error function of its argument |
14 |
|
|
erfc(x) returns 1.0-erf(x) |
15 |
|
|
|
16 |
|
|
erf(x) is defined by |
17 |
|
|
${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$ |
18 |
|
|
|
19 |
|
|
the entry for erfc is provided because of the |
20 |
|
|
extreme loss of relative accuracy if erf(x) is |
21 |
|
|
called for large x and the result subtracted |
22 |
|
|
from 1. (e.g. for x= 10, 12 places are lost). |
23 |
|
|
|
24 |
|
|
There are no error returns. |
25 |
|
|
|
26 |
|
|
Calls exp. |
27 |
|
|
|
28 |
|
|
Coefficients for large x are #5667 from Hart & Cheney (18.72D). |
29 |
|
|
*/ |
30 |
|
|
|
31 |
|
|
#define M 7 |
32 |
|
|
#define N 9 |
33 |
|
|
static double torp = 1.1283791670955125738961589031; |
34 |
|
|
static double p1[] = { |
35 |
|
|
0.804373630960840172832162e5, |
36 |
|
|
0.740407142710151470082064e4, |
37 |
|
|
0.301782788536507577809226e4, |
38 |
|
|
0.380140318123903008244444e2, |
39 |
|
|
0.143383842191748205576712e2, |
40 |
|
|
-.288805137207594084924010e0, |
41 |
|
|
0.007547728033418631287834e0, |
42 |
|
|
}; |
43 |
|
|
static double q1[] = { |
44 |
|
|
0.804373630960840172826266e5, |
45 |
|
|
0.342165257924628539769006e5, |
46 |
|
|
0.637960017324428279487120e4, |
47 |
|
|
0.658070155459240506326937e3, |
48 |
|
|
0.380190713951939403753468e2, |
49 |
|
|
0.100000000000000000000000e1, |
50 |
|
|
0.0, |
51 |
|
|
}; |
52 |
|
|
static double p2[] = { |
53 |
|
|
0.18263348842295112592168999e4, |
54 |
|
|
0.28980293292167655611275846e4, |
55 |
|
|
0.2320439590251635247384768711e4, |
56 |
|
|
0.1143262070703886173606073338e4, |
57 |
|
|
0.3685196154710010637133875746e3, |
58 |
|
|
0.7708161730368428609781633646e2, |
59 |
|
|
0.9675807882987265400604202961e1, |
60 |
|
|
0.5641877825507397413087057563e0, |
61 |
|
|
0.0, |
62 |
|
|
}; |
63 |
|
|
static double q2[] = { |
64 |
|
|
0.18263348842295112595576438e4, |
65 |
|
|
0.495882756472114071495438422e4, |
66 |
|
|
0.60895424232724435504633068e4, |
67 |
|
|
0.4429612803883682726711528526e4, |
68 |
|
|
0.2094384367789539593790281779e4, |
69 |
|
|
0.6617361207107653469211984771e3, |
70 |
|
|
0.1371255960500622202878443578e3, |
71 |
|
|
0.1714980943627607849376131193e2, |
72 |
|
|
1.0, |
73 |
|
|
}; |
74 |
|
|
|
75 |
|
|
double |
76 |
greg |
3.2 |
erf(double arg) { |
77 |
greg |
3.1 |
int sign; |
78 |
|
|
double argsq; |
79 |
|
|
double d, n; |
80 |
|
|
int i; |
81 |
|
|
|
82 |
|
|
sign = 1; |
83 |
|
|
if(arg < 0.){ |
84 |
|
|
arg = -arg; |
85 |
|
|
sign = -1; |
86 |
|
|
} |
87 |
|
|
if(arg < 0.5){ |
88 |
|
|
argsq = arg*arg; |
89 |
|
|
for(n=0,d=0,i=M-1; i>=0; i--){ |
90 |
|
|
n = n*argsq + p1[i]; |
91 |
|
|
d = d*argsq + q1[i]; |
92 |
|
|
} |
93 |
|
|
return(sign*torp*arg*n/d); |
94 |
|
|
} |
95 |
|
|
if(arg >= 10.) |
96 |
|
|
return(sign*1.); |
97 |
|
|
return(sign*(1. - erfc(arg))); |
98 |
|
|
} |
99 |
|
|
|
100 |
|
|
double |
101 |
greg |
3.2 |
erfc(double arg) { |
102 |
greg |
3.1 |
double n, d; |
103 |
|
|
int i; |
104 |
|
|
|
105 |
|
|
if(arg < 0.) |
106 |
|
|
return(2. - erfc(-arg)); |
107 |
|
|
/* |
108 |
|
|
if(arg < 0.5) |
109 |
|
|
return(1. - erf(arg)); |
110 |
|
|
*/ |
111 |
|
|
if(arg >= 10.) |
112 |
|
|
return(0.); |
113 |
|
|
|
114 |
|
|
for(n=0,d=0,i=N-1; i>=0; i--){ |
115 |
|
|
n = n*arg + p2[i]; |
116 |
|
|
d = d*arg + q2[i]; |
117 |
|
|
} |
118 |
|
|
return(exp(-arg*arg)*n/d); |
119 |
|
|
} |