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