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