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