Back to index

radiance  4R0+20100331
erf.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: erf.c,v 3.1 2003/02/22 02:07:22 greg Exp $";
00003 #endif
00004 #ifndef lint
00005 static char sccsid[] = "@(#)erf.c 1.1 87/12/21 SMI"; /* from UCB 4.1 12/25/82 */
00006 #endif
00007 
00008 /*
00009        C program for floating point error function
00010 
00011        erf(x) returns the error function of its argument
00012        erfc(x) returns 1.0-erf(x)
00013 
00014        erf(x) is defined by
00015        ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
00016 
00017        the entry for erfc is provided because of the
00018        extreme loss of relative accuracy if erf(x) is
00019        called for large x and the result subtracted
00020        from 1. (e.g. for x= 10, 12 places are lost).
00021 
00022        There are no error returns.
00023 
00024        Calls exp.
00025 
00026        Coefficients for large x are #5667 from Hart & Cheney (18.72D).
00027 */
00028 
00029 #define M 7
00030 #define N 9
00031 static double torp = 1.1283791670955125738961589031;
00032 static double p1[] = {
00033        0.804373630960840172832162e5,
00034        0.740407142710151470082064e4,
00035        0.301782788536507577809226e4,
00036        0.380140318123903008244444e2,
00037        0.143383842191748205576712e2,
00038        -.288805137207594084924010e0,
00039        0.007547728033418631287834e0,
00040 };
00041 static double q1[]  = {
00042        0.804373630960840172826266e5,
00043        0.342165257924628539769006e5,
00044        0.637960017324428279487120e4,
00045        0.658070155459240506326937e3,
00046        0.380190713951939403753468e2,
00047        0.100000000000000000000000e1,
00048        0.0,
00049 };
00050 static double p2[]  = {
00051        0.18263348842295112592168999e4,
00052        0.28980293292167655611275846e4,
00053        0.2320439590251635247384768711e4,
00054        0.1143262070703886173606073338e4,
00055        0.3685196154710010637133875746e3,
00056        0.7708161730368428609781633646e2,
00057        0.9675807882987265400604202961e1,
00058        0.5641877825507397413087057563e0,
00059        0.0,
00060 };
00061 static double q2[]  = {
00062        0.18263348842295112595576438e4,
00063        0.495882756472114071495438422e4,
00064        0.60895424232724435504633068e4,
00065        0.4429612803883682726711528526e4,
00066        0.2094384367789539593790281779e4,
00067        0.6617361207107653469211984771e3,
00068        0.1371255960500622202878443578e3,
00069        0.1714980943627607849376131193e2,
00070        1.0,
00071 };
00072 
00073 double
00074 erf(arg) double arg;{
00075        double erfc();
00076        int sign;
00077        double argsq;
00078        double d, n;
00079        int i;
00080 
00081        sign = 1;
00082        if(arg < 0.){
00083               arg = -arg;
00084               sign = -1;
00085        }
00086        if(arg < 0.5){
00087               argsq = arg*arg;
00088               for(n=0,d=0,i=M-1; i>=0; i--){
00089                      n = n*argsq + p1[i];
00090                      d = d*argsq + q1[i];
00091               }
00092               return(sign*torp*arg*n/d);
00093        }
00094        if(arg >= 10.)
00095               return(sign*1.);
00096        return(sign*(1. - erfc(arg)));
00097 }
00098 
00099 double
00100 erfc(arg) double arg;{
00101        double erf();
00102        double exp();
00103        double n, d;
00104        int i;
00105 
00106        if(arg < 0.)
00107               return(2. - erfc(-arg));
00108 /*
00109        if(arg < 0.5)
00110               return(1. - erf(arg));
00111 */
00112        if(arg >= 10.)
00113               return(0.);
00114 
00115        for(n=0,d=0,i=N-1; i>=0; i--){
00116               n = n*arg + p2[i];
00117               d = d*arg + q2[i];
00118        }
00119        return(exp(-arg*arg)*n/d);
00120 }