1*48402Sbostic /*- 2*48402Sbostic * Copyright (c) 1985 The Regents of the University of California. 3*48402Sbostic * All rights reserved. 4*48402Sbostic * 5*48402Sbostic * %sccs.include.proprietary.c% 634121Sbostic */ 724592Szliu 834121Sbostic #ifndef lint 9*48402Sbostic static char sccsid[] = "@(#)erf.c 5.3 (Berkeley) 04/20/91"; 1034121Sbostic #endif /* not lint */ 1134121Sbostic 1224592Szliu /* 1324592Szliu C program for floating point error function 1424592Szliu 1524592Szliu erf(x) returns the error function of its argument 1624592Szliu erfc(x) returns 1.0-erf(x) 1724592Szliu 1824592Szliu erf(x) is defined by 1924592Szliu ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$ 2024592Szliu 2124592Szliu the entry for erfc is provided because of the 2224592Szliu extreme loss of relative accuracy if erf(x) is 2324592Szliu called for large x and the result subtracted 2424592Szliu from 1. (e.g. for x= 10, 12 places are lost). 2524592Szliu 2624592Szliu There are no error returns. 2724592Szliu 2824592Szliu Calls exp. 2924592Szliu 3024592Szliu Coefficients for large x are #5667 from Hart & Cheney (18.72D). 3124592Szliu */ 3224592Szliu 3324592Szliu #define M 7 3424592Szliu #define N 9 3524592Szliu static double torp = 1.1283791670955125738961589031; 3624592Szliu static double p1[] = { 3724592Szliu 0.804373630960840172832162e5, 3824592Szliu 0.740407142710151470082064e4, 3924592Szliu 0.301782788536507577809226e4, 4024592Szliu 0.380140318123903008244444e2, 4124592Szliu 0.143383842191748205576712e2, 4224592Szliu -.288805137207594084924010e0, 4324592Szliu 0.007547728033418631287834e0, 4424592Szliu }; 4524592Szliu static double q1[] = { 4624592Szliu 0.804373630960840172826266e5, 4724592Szliu 0.342165257924628539769006e5, 4824592Szliu 0.637960017324428279487120e4, 4924592Szliu 0.658070155459240506326937e3, 5024592Szliu 0.380190713951939403753468e2, 5124592Szliu 0.100000000000000000000000e1, 5224592Szliu 0.0, 5324592Szliu }; 5424592Szliu static double p2[] = { 5524592Szliu 0.18263348842295112592168999e4, 5624592Szliu 0.28980293292167655611275846e4, 5724592Szliu 0.2320439590251635247384768711e4, 5824592Szliu 0.1143262070703886173606073338e4, 5924592Szliu 0.3685196154710010637133875746e3, 6024592Szliu 0.7708161730368428609781633646e2, 6124592Szliu 0.9675807882987265400604202961e1, 6224592Szliu 0.5641877825507397413087057563e0, 6324592Szliu 0.0, 6424592Szliu }; 6524592Szliu static double q2[] = { 6624592Szliu 0.18263348842295112595576438e4, 6724592Szliu 0.495882756472114071495438422e4, 6824592Szliu 0.60895424232724435504633068e4, 6924592Szliu 0.4429612803883682726711528526e4, 7024592Szliu 0.2094384367789539593790281779e4, 7124592Szliu 0.6617361207107653469211984771e3, 7224592Szliu 0.1371255960500622202878443578e3, 7324592Szliu 0.1714980943627607849376131193e2, 7424592Szliu 1.0, 7524592Szliu }; 7624592Szliu 7724592Szliu double 7824592Szliu erf(arg) double arg;{ 7924592Szliu double erfc(); 8024592Szliu int sign; 8124592Szliu double argsq; 8224592Szliu double d, n; 8324592Szliu int i; 8424592Szliu 8524592Szliu sign = 1; 8624592Szliu if(arg < 0.){ 8724592Szliu arg = -arg; 8824592Szliu sign = -1; 8924592Szliu } 9024592Szliu if(arg < 0.5){ 9124592Szliu argsq = arg*arg; 9224592Szliu for(n=0,d=0,i=M-1; i>=0; i--){ 9324592Szliu n = n*argsq + p1[i]; 9424592Szliu d = d*argsq + q1[i]; 9524592Szliu } 9624592Szliu return(sign*torp*arg*n/d); 9724592Szliu } 9824592Szliu if(arg >= 10.) 9924592Szliu return(sign*1.); 10024592Szliu return(sign*(1. - erfc(arg))); 10124592Szliu } 10224592Szliu 10324592Szliu double 10424592Szliu erfc(arg) double arg;{ 10524592Szliu double erf(); 10624592Szliu double exp(); 10724592Szliu double n, d; 10824592Szliu int i; 10924592Szliu 11024592Szliu if(arg < 0.) 11124592Szliu return(2. - erfc(-arg)); 11224592Szliu /* 11324592Szliu if(arg < 0.5) 11424592Szliu return(1. - erf(arg)); 11524592Szliu */ 11624592Szliu if(arg >= 10.) 11724592Szliu return(0.); 11824592Szliu 11924592Szliu for(n=0,d=0,i=N-1; i>=0; i--){ 12024592Szliu n = n*arg + p2[i]; 12124592Szliu d = d*arg + q2[i]; 12224592Szliu } 12324592Szliu return(exp(-arg*arg)*n/d); 12424592Szliu } 125