Back to index

glibc  2.9
Functions | Variables
s_erf.c File Reference
#include "math.h"
#include "math_private.h"
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

double __erf (double x)
 weak_alias (__erf, erf)

Variables

static double tiny = 1e-300
static double half = 5.00000000000000000000e-01
static double one = 1.00000000000000000000e+00
static double two = 2.00000000000000000000e+00
static double erx = 8.45062911510467529297e-01
static double efx = 1.28379167095512586316e-01
static double efx8 = 1.02703333676410069053e+00
static double pp []
static double qq []
static double pa []
static double qa []
static double ra []
static double sa []
static double rb []
static double sb []

Function Documentation

double __erf ( double  x)

Definition at line 200 of file s_erf.c.

{
       int32_t hx,ix,i;
       double R,S,P,Q,s,y,z,r;
       GET_HIGH_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix>=0x7ff00000) {        /* erf(nan)=nan */
           i = ((u_int32_t)hx>>31)<<1;
           return (double)(1-i)+one/x;    /* erf(+-inf)=+-1 */
       }

       if(ix < 0x3feb0000) {              /* |x|<0.84375 */
           double r1,r2,s1,s2,s3,z2,z4;
           if(ix < 0x3e300000) {   /* |x|<2**-28 */
               if (ix < 0x00800000)
                  return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
              return x + efx*x;
           }
           z = x*x;
#ifdef DO_NOT_USE_THIS
           r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
           s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
#else
           r1 = pp[0]+z*pp[1]; z2=z*z;
           r2 = pp[2]+z*pp[3]; z4=z2*z2;
           s1 = one+z*qq[1];
           s2 = qq[2]+z*qq[3];
           s3 = qq[4]+z*qq[5];
            r = r1 + z2*r2 + z4*pp[4];
           s  = s1 + z2*s2 + z4*s3;
#endif
           y = r/s;
           return x + x*y;
       }
       if(ix < 0x3ff40000) {              /* 0.84375 <= |x| < 1.25 */
           double s2,s4,s6,P1,P2,P3,P4,Q1,Q2,Q3,Q4;
           s = fabs(x)-one;
#ifdef DO_NOT_USE_THIS
           P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
           Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
#else
           P1 = pa[0]+s*pa[1]; s2=s*s;
           Q1 = one+s*qa[1];   s4=s2*s2;
           P2 = pa[2]+s*pa[3]; s6=s4*s2;
           Q2 = qa[2]+s*qa[3];
           P3 = pa[4]+s*pa[5];
           Q3 = qa[4]+s*qa[5];
           P4 = pa[6];
           Q4 = qa[6];
           P = P1 + s2*P2 + s4*P3 + s6*P4;
           Q = Q1 + s2*Q2 + s4*Q3 + s6*Q4;
#endif
           if(hx>=0) return erx + P/Q; else return -erx - P/Q;
       }
       if (ix >= 0x40180000) {            /* inf>|x|>=6 */
           if(hx>=0) return one-tiny; else return tiny-one;
       }
       x = fabs(x);
       s = one/(x*x);
       if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */
#ifdef DO_NOT_USE_THIS
           R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
                            ra5+s*(ra6+s*ra7))))));
           S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
                            sa5+s*(sa6+s*(sa7+s*sa8)))))));
#else
           double R1,R2,R3,R4,S1,S2,S3,S4,s2,s4,s6,s8;
           R1 = ra[0]+s*ra[1];s2 = s*s;
           S1 = one+s*sa[1];  s4 = s2*s2;
           R2 = ra[2]+s*ra[3];s6 = s4*s2;
           S2 = sa[2]+s*sa[3];s8 = s4*s4;
           R3 = ra[4]+s*ra[5];
           S3 = sa[4]+s*sa[5];
           R4 = ra[6]+s*ra[7];
           S4 = sa[6]+s*sa[7];
           R = R1 + s2*R2 + s4*R3 + s6*R4;
           S = S1 + s2*S2 + s4*S3 + s6*S4 + s8*sa[8];
#endif
       } else {      /* |x| >= 1/0.35 */
#ifdef DO_NOT_USE_THIS
           R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
                            rb5+s*rb6)))));
           S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
                            sb5+s*(sb6+s*sb7))))));
#else
           double R1,R2,R3,S1,S2,S3,S4,s2,s4,s6;
           R1 = rb[0]+s*rb[1];s2 = s*s;
           S1 = one+s*sb[1];  s4 = s2*s2;
           R2 = rb[2]+s*rb[3];s6 = s4*s2;
           S2 = sb[2]+s*sb[3];
           R3 = rb[4]+s*rb[5];
           S3 = sb[4]+s*sb[5];
           S4 = sb[6]+s*sb[7];
           R = R1 + s2*R2 + s4*R3 + s6*rb[6];
           S = S1 + s2*S2 + s4*S3 + s6*S4;
#endif
       }
       z  = x;
       SET_LOW_WORD(z,0);
       r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
       if(hx>=0) return one-r/x; else return  r/x-one;
}

Here is the call graph for this function:

weak_alias ( __erf  ,
erf   
)

Definition at line 304 of file s_erf.c.

{
       int32_t hx,ix;
       double R,S,P,Q,s,y,z,r;
       GET_HIGH_WORD(hx,x);
       ix = hx&0x7fffffff;
       if(ix>=0x7ff00000) {               /* erfc(nan)=nan */
                                          /* erfc(+-inf)=0,2 */
           return (double)(((u_int32_t)hx>>31)<<1)+one/x;
       }

       if(ix < 0x3feb0000) {              /* |x|<0.84375 */
           double r1,r2,s1,s2,s3,z2,z4;
           if(ix < 0x3c700000)     /* |x|<2**-56 */
              return one-x;
           z = x*x;
#ifdef DO_NOT_USE_THIS
           r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
           s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
#else
           r1 = pp[0]+z*pp[1]; z2=z*z;
           r2 = pp[2]+z*pp[3]; z4=z2*z2;
           s1 = one+z*qq[1];
           s2 = qq[2]+z*qq[3];
           s3 = qq[4]+z*qq[5];
            r = r1 + z2*r2 + z4*pp[4];
           s  = s1 + z2*s2 + z4*s3;
#endif
           y = r/s;
           if(hx < 0x3fd00000) {   /* x<1/4 */
              return one-(x+x*y);
           } else {
              r = x*y;
              r += (x-half);
               return half - r ;
           }
       }
       if(ix < 0x3ff40000) {              /* 0.84375 <= |x| < 1.25 */
           double s2,s4,s6,P1,P2,P3,P4,Q1,Q2,Q3,Q4;
           s = fabs(x)-one;
#ifdef DO_NOT_USE_THIS
           P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
           Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
#else
           P1 = pa[0]+s*pa[1]; s2=s*s;
           Q1 = one+s*qa[1];   s4=s2*s2;
           P2 = pa[2]+s*pa[3]; s6=s4*s2;
           Q2 = qa[2]+s*qa[3];
           P3 = pa[4]+s*pa[5];
           Q3 = qa[4]+s*qa[5];
           P4 = pa[6];
           Q4 = qa[6];
           P = P1 + s2*P2 + s4*P3 + s6*P4;
           Q = Q1 + s2*Q2 + s4*Q3 + s6*Q4;
#endif
           if(hx>=0) {
               z  = one-erx; return z - P/Q;
           } else {
              z = erx+P/Q; return one+z;
           }
       }
       if (ix < 0x403c0000) {             /* |x|<28 */
           x = fabs(x);
           s = one/(x*x);
           if(ix< 0x4006DB6D) {    /* |x| < 1/.35 ~ 2.857143*/
#ifdef DO_NOT_USE_THIS
               R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
                            ra5+s*(ra6+s*ra7))))));
               S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
                            sa5+s*(sa6+s*(sa7+s*sa8)))))));
#else
              double R1,R2,R3,R4,S1,S2,S3,S4,s2,s4,s6,s8;
           R1 = ra[0]+s*ra[1];s2 = s*s;
           S1 = one+s*sa[1];  s4 = s2*s2;
           R2 = ra[2]+s*ra[3];s6 = s4*s2;
           S2 = sa[2]+s*sa[3];s8 = s4*s4;
           R3 = ra[4]+s*ra[5];
           S3 = sa[4]+s*sa[5];
           R4 = ra[6]+s*ra[7];
           S4 = sa[6]+s*sa[7];
           R = R1 + s2*R2 + s4*R3 + s6*R4;
           S = S1 + s2*S2 + s4*S3 + s6*S4 + s8*sa[8];
#endif
           } else {                /* |x| >= 1/.35 ~ 2.857143 */
              double R1,R2,R3,S1,S2,S3,S4,s2,s4,s6;
              if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
#ifdef DO_NOT_USE_THIS
               R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
                            rb5+s*rb6)))));
               S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
                            sb5+s*(sb6+s*sb7))))));
#else
              R1 = rb[0]+s*rb[1];s2 = s*s;
              S1 = one+s*sb[1];  s4 = s2*s2;
              R2 = rb[2]+s*rb[3];s6 = s4*s2;
              S2 = sb[2]+s*sb[3];
              R3 = rb[4]+s*rb[5];
              S3 = sb[4]+s*sb[5];
              S4 = sb[6]+s*sb[7];
              R = R1 + s2*R2 + s4*R3 + s6*rb[6];
              S = S1 + s2*S2 + s4*S3 + s6*S4;
#endif
           }
           z  = x;
           SET_LOW_WORD(z,0);
           r  =  __ieee754_exp(-z*z-0.5625)*
                     __ieee754_exp((z-x)*(z+x)+R/S);
           if(hx>0) return r/x; else return two-r/x;
       } else {
           if(hx>0) return tiny*tiny; else return two-tiny;
       }
}

Here is the call graph for this function:


Variable Documentation

double efx = 1.28379167095512586316e-01 [static]

Definition at line 132 of file s_erf.c.

double efx8 = 1.02703333676410069053e+00 [static]

Definition at line 133 of file s_erf.c.

double erx = 8.45062911510467529297e-01 [static]

Definition at line 128 of file s_erf.c.

double half = 5.00000000000000000000e-01 [static]

Definition at line 124 of file s_erf.c.

double one = 1.00000000000000000000e+00 [static]

Definition at line 125 of file s_erf.c.

double pa[] [static]
Initial value:
 {-2.36211856075265944077e-03, 
  4.14856118683748331666e-01, 
 -3.72207876035701323847e-01, 
  3.18346619901161753674e-01, 
 -1.10894694282396677476e-01, 
  3.54783043256182359371e-02, 
 -2.16637559486879084300e-03}

Definition at line 147 of file s_erf.c.

double pp[] [static]
Initial value:
  {1.28379167095512558561e-01, 
 -3.25042107247001499370e-01, 
 -2.84817495755985104766e-02, 
 -5.77027029648944159157e-03, 
 -2.37630166566501626084e-05}

Definition at line 134 of file s_erf.c.

double qa[] [static]
Initial value:
  {0.0, 1.06420880400844228286e-01, 
  5.40397917702171048937e-01, 
  7.18286544141962662868e-02, 
  1.26171219808761642112e-01, 
  1.36370839120290507362e-02, 
  1.19844998467991074170e-02}

Definition at line 154 of file s_erf.c.

double qq[] [static]
Initial value:
  {0.0, 3.97917223959155352819e-01, 
  6.50222499887672944485e-02, 
  5.08130628187576562776e-03, 
  1.32494738004321644526e-04, 
 -3.96022827877536812320e-06}

Definition at line 139 of file s_erf.c.

double ra[] [static]
Initial value:
 {-9.86494403484714822705e-03, 
 -6.93858572707181764372e-01, 
 -1.05586262253232909814e+01, 
 -6.23753324503260060396e+01, 
 -1.62396669462573470355e+02, 
 -1.84605092906711035994e+02, 
 -8.12874355063065934246e+01, 
 -9.81432934416914548592e+00}

Definition at line 163 of file s_erf.c.

double rb[] [static]
Initial value:
 {-9.86494292470009928597e-03, 
 -7.99283237680523006574e-01, 
 -1.77579549177547519889e+01, 
 -1.60636384855821916062e+02, 
 -6.37566443368389627722e+02, 
 -1.02509513161107724954e+03, 
 -4.83519191608651397019e+02}

Definition at line 182 of file s_erf.c.

double sa[] [static]
Initial value:
  {0.0,1.96512716674392571292e+01, 
  1.37657754143519042600e+02, 
  4.34565877475229228821e+02, 
  6.45387271733267880336e+02, 
  4.29008140027567833386e+02, 
  1.08635005541779435134e+02, 
  6.57024977031928170135e+00, 
 -6.04244152148580987438e-02}

Definition at line 171 of file s_erf.c.

double sb[] [static]
Initial value:
  {0.0,3.03380607434824582924e+01, 
  3.25792512996573918826e+02, 
  1.53672958608443695994e+03, 
  3.19985821950859553908e+03, 
  2.55305040643316442583e+03, 
  4.74528541206955367215e+02, 
 -2.24409524465858183362e+01}

Definition at line 189 of file s_erf.c.

double tiny = 1e-300 [static]

Definition at line 123 of file s_erf.c.

double two = 2.00000000000000000000e+00 [static]

Definition at line 126 of file s_erf.c.