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