libpok-refman-current-html/html/erff_8c_source.html

00001 /*
00002  *                               POK header
00003  *
00004  * The following file is a part of the POK project. Any modification should
00005  * made according to the POK licence. You CANNOT use this file or a part of
00006  * this file is this part of a file for your own project
00007  *
00008  * For more information on the POK licence, please see our LICENCE FILE
00009  *
00010  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
00011  *
00012  *                                      Copyright (c) 2007-2009 POK team
00013  *
00014  * Created by julien on Fri Jan 30 14:41:34 2009
00015  */
00016
00017 /* s_erff.c -- float version of s_erf.c.
00018  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00019  */
00020
00021 /*
00022  * ====================================================
00023  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00024  *
00025  * Developed at SunPro, a Sun Microsystems, Inc. business.
00026  * Permission to use, copy, modify, and distribute this
00027  * software is freely granted, provided that this notice
00028  * is preserved.
00029  * ====================================================
00030  */
00031
00032 #ifdef POK_NEEDS_LIBMATH
00033
00034 #include <libm.h>
00035 #include "math_private.h"
00036
00037 static const float
00038 tiny        = 1e-30,
00039 half=  5.0000000000e-01, /* 0x3F000000 */
00040 one =  1.0000000000e+00, /* 0x3F800000 */
00041 two =  2.0000000000e+00, /* 0x40000000 */
00042         /* c = (subfloat)0.84506291151 */
00043 erx =  8.4506291151e-01, /* 0x3f58560b */
00044 /*
00045  * Coefficients for approximation to  erf on [0,0.84375]
00046  */
00047 efx =  1.2837916613e-01, /* 0x3e0375d4 */
00048 efx8=  1.0270333290e+00, /* 0x3f8375d4 */
00049 pp0  =  1.2837916613e-01, /* 0x3e0375d4 */
00050 pp1  = -3.2504209876e-01, /* 0xbea66beb */
00051 pp2  = -2.8481749818e-02, /* 0xbce9528f */
00052 pp3  = -5.7702702470e-03, /* 0xbbbd1489 */
00053 pp4  = -2.3763017452e-05, /* 0xb7c756b1 */
00054 qq1  =  3.9791721106e-01, /* 0x3ecbbbce */
00055 qq2  =  6.5022252500e-02, /* 0x3d852a63 */
00056 qq3  =  5.0813062117e-03, /* 0x3ba68116 */
00057 qq4  =  1.3249473704e-04, /* 0x390aee49 */
00058 qq5  = -3.9602282413e-06, /* 0xb684e21a */
00059 /*
00060  * Coefficients for approximation to  erf  in [0.84375,1.25]
00061  */
00062 pa0  = -2.3621185683e-03, /* 0xbb1acdc6 */
00063 pa1  =  4.1485610604e-01, /* 0x3ed46805 */
00064 pa2  = -3.7220788002e-01, /* 0xbebe9208 */
00065 pa3  =  3.1834661961e-01, /* 0x3ea2fe54 */
00066 pa4  = -1.1089469492e-01, /* 0xbde31cc2 */
00067 pa5  =  3.5478305072e-02, /* 0x3d1151b3 */
00068 pa6  = -2.1663755178e-03, /* 0xbb0df9c0 */
00069 qa1  =  1.0642088205e-01, /* 0x3dd9f331 */
00070 qa2  =  5.4039794207e-01, /* 0x3f0a5785 */
00071 qa3  =  7.1828655899e-02, /* 0x3d931ae7 */
00072 qa4  =  1.2617121637e-01, /* 0x3e013307 */
00073 qa5  =  1.3637083583e-02, /* 0x3c5f6e13 */
00074 qa6  =  1.1984500103e-02, /* 0x3c445aa3 */
00075 /*
00076  * Coefficients for approximation to  erfc in [1.25,1/0.35]
00077  */
00078 ra0  = -9.8649440333e-03, /* 0xbc21a093 */
00079 ra1  = -6.9385856390e-01, /* 0xbf31a0b7 */
00080 ra2  = -1.0558626175e+01, /* 0xc128f022 */
00081 ra3  = -6.2375331879e+01, /* 0xc2798057 */
00082 ra4  = -1.6239666748e+02, /* 0xc322658c */
00083 ra5  = -1.8460508728e+02, /* 0xc3389ae7 */
00084 ra6  = -8.1287437439e+01, /* 0xc2a2932b */
00085 ra7  = -9.8143291473e+00, /* 0xc11d077e */
00086 sa1  =  1.9651271820e+01, /* 0x419d35ce */
00087 sa2  =  1.3765776062e+02, /* 0x4309a863 */
00088 sa3  =  4.3456588745e+02, /* 0x43d9486f */
00089 sa4  =  6.4538726807e+02, /* 0x442158c9 */
00090 sa5  =  4.2900814819e+02, /* 0x43d6810b */
00091 sa6  =  1.0863500214e+02, /* 0x42d9451f */
00092 sa7  =  6.5702495575e+00, /* 0x40d23f7c */
00093 sa8  = -6.0424413532e-02, /* 0xbd777f97 */
00094 /*
00095  * Coefficients for approximation to  erfc in [1/.35,28]
00096  */
00097 rb0  = -9.8649431020e-03, /* 0xbc21a092 */
00098 rb1  = -7.9928326607e-01, /* 0xbf4c9dd4 */
00099 rb2  = -1.7757955551e+01, /* 0xc18e104b */
00100 rb3  = -1.6063638306e+02, /* 0xc320a2ea */
00101 rb4  = -6.3756646729e+02, /* 0xc41f6441 */
00102 rb5  = -1.0250950928e+03, /* 0xc480230b */
00103 rb6  = -4.8351919556e+02, /* 0xc3f1c275 */
00104 sb1  =  3.0338060379e+01, /* 0x41f2b459 */
00105 sb2  =  3.2579251099e+02, /* 0x43a2e571 */
00106 sb3  =  1.5367296143e+03, /* 0x44c01759 */
00107 sb4  =  3.1998581543e+03, /* 0x4547fdbb */
00108 sb5  =  2.5530502930e+03, /* 0x451f90ce */
00109 sb6  =  4.7452853394e+02, /* 0x43ed43a7 */
00110 sb7  = -2.2440952301e+01; /* 0xc1b38712 */
00111
00112 float
00113 erff(float x)
00114 {
00115         int32_t hx,ix,i;
00116         float R,S,P,Q,s,y,z,r;
00117         GET_FLOAT_WORD(hx,x);
00118         ix = hx&0x7fffffff;
00119         if(ix>=0x7f800000) {            /* erf(nan)=nan */
00120             i = ((uint32_t)hx>>31)<<1;
00121             return (float)(1-i)+one/x;  /* erf(+-inf)=+-1 */
00122         }
00123
00124         if(ix < 0x3f580000) {           /* |x|<0.84375 */
00125             if(ix < 0x31800000) {       /* |x|<2**-28 */
00126                 if (ix < 0x04000000)
00127                     /*avoid underflow */
00128                     return (float)0.125*((float)8.0*x+efx8*x);
00129                 return x + efx*x;
00130             }
00131             z = x*x;
00132             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
00133             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
00134             y = r/s;
00135             return x + x*y;
00136         }
00137         if(ix < 0x3fa00000) {           /* 0.84375 <= |x| < 1.25 */
00138             s = fabsf(x)-one;
00139             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
00140             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
00141             if(hx>=0) return erx + P/Q; else return -erx - P/Q;
00142         }
00143         if (ix >= 0x40c00000) {         /* inf>|x|>=6 */
00144             if(hx>=0) return one-tiny; else return tiny-one;
00145         }
00146         x = fabsf(x);
00147         s = one/(x*x);
00148         if(ix< 0x4036DB6E) {    /* |x| < 1/0.35 */
00149             R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
00150                                 ra5+s*(ra6+s*ra7))))));
00151             S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
00152                                 sa5+s*(sa6+s*(sa7+s*sa8)))))));
00153         } else {        /* |x| >= 1/0.35 */
00154             R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
00155                                 rb5+s*rb6)))));
00156             S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
00157                                 sb5+s*(sb6+s*sb7))))));
00158         }
00159         GET_FLOAT_WORD(ix,x);
00160         SET_FLOAT_WORD(z,ix&0xfffff000);
00161         r  =  __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S);
00162         if(hx>=0) return one-r/x; else return  r/x-one;
00163 }
00164
00165 float
00166 erfcf(float x)
00167 {
00168         int32_t hx,ix;
00169         float R,S,P,Q,s,y,z,r;
00170         GET_FLOAT_WORD(hx,x);
00171         ix = hx&0x7fffffff;
00172         if(ix>=0x7f800000) {                    /* erfc(nan)=nan */
00173                                                 /* erfc(+-inf)=0,2 */
00174             return (float)(((uint32_t)hx>>31)<<1)+one/x;
00175         }
00176
00177         if(ix < 0x3f580000) {           /* |x|<0.84375 */
00178             if(ix < 0x23800000)         /* |x|<2**-56 */
00179                 return one-x;
00180             z = x*x;
00181             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
00182             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
00183             y = r/s;
00184             if(hx < 0x3e800000) {       /* x<1/4 */
00185                 return one-(x+x*y);
00186             } else {
00187                 r = x*y;
00188                 r += (x-half);
00189                 return half - r ;
00190             }
00191         }
00192         if(ix < 0x3fa00000) {           /* 0.84375 <= |x| < 1.25 */
00193             s = fabsf(x)-one;
00194             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
00195             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
00196             if(hx>=0) {
00197                 z  = one-erx; return z - P/Q;
00198             } else {
00199                 z = erx+P/Q; return one+z;
00200             }
00201         }
00202         if (ix < 0x41e00000) {          /* |x|<28 */
00203             x = fabsf(x);
00204             s = one/(x*x);
00205             if(ix< 0x4036DB6D) {        /* |x| < 1/.35 ~ 2.857143*/
00206                 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
00207                                 ra5+s*(ra6+s*ra7))))));
00208                 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
00209                                 sa5+s*(sa6+s*(sa7+s*sa8)))))));
00210             } else {                    /* |x| >= 1/.35 ~ 2.857143 */
00211                 if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
00212                 R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
00213                                 rb5+s*rb6)))));
00214                 S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
00215                                 sb5+s*(sb6+s*sb7))))));
00216             }
00217             GET_FLOAT_WORD(ix,x);
00218             SET_FLOAT_WORD(z,ix&0xfffff000);
00219             r  =  __ieee754_expf(-z*z-(float)0.5625)*
00220                         __ieee754_expf((z-x)*(z+x)+R/S);
00221             if(hx>0) return r/x; else return two-r/x;
00222         } else {
00223             if(hx>0) return tiny*tiny; else return two-tiny;
00224         }
00225 }
00226
00227 #endif