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