1 /*- 2 * Copyright (c) 1992 The Regents of the University of California. 3 * All rights reserved. 4 * 5 * %sccs.include.redist.c% 6 */ 7 8 #ifndef lint 9 static char sccsid[] = "@(#)erf.c 5.4 (Berkeley) 12/02/92"; 10 #endif /* not lint */ 11 12 /* 13 * ==================================================== 14 * Copyright (C) 1992 by Sun Microsystems, Inc. 15 * 16 * Developed at SunPro, a Sun Microsystems, Inc. business. 17 * Permission to use, copy, modify, and distribute this 18 * software is freely granted, provided that this notice 19 * is preserved. 20 * ==================================================== 21 * 22 * ******************* WARNING ******************** 23 * This is an alpha version of SunPro's FDLIBM (Freely 24 * Distributable Math Library) for IEEE double precision 25 * arithmetic. FDLIBM is a basic math library written 26 * in C that runs on machines that conform to IEEE 27 * Standard 754/854. This alpha version is distributed 28 * for testing purpose. Those who use this software 29 * should report any bugs to 30 * 31 * fdlibm-comments@sunpro.eng.sun.com 32 * 33 * -- K.C. Ng, Oct 12, 1992 34 * ************************************************ 35 */ 36 37 /* Modified Nov 30, 1992 P. McILROY: 38 * Replaced expansion for x > 6 39 * Add #ifdef's for vax/tahoe. 40 */ 41 42 /* double erf(double x) 43 * double erfc(double x) 44 * x 45 * 2 |\ 46 * erf(x) = --------- | exp(-t*t)dt 47 * sqrt(pi) \| 48 * 0 49 * 50 * erfc(x) = 1-erf(x) 51 * 52 * Method: 53 * 1. Reduce x to |x| by erf(-x) = -erf(x) 54 * 2. For x in [0, 0.84375] 55 * erf(x) = x + x*P(x^2) 56 * erfc(x) = 1 - erf(x) if x<=0.25 57 * = 0.5 + ((0.5-x)-x*P) if x in [0.25,0.84375] 58 * where 59 * 2 2 4 20 60 * P = P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x ) 61 * is an approximation to (erf(x)-x)/x with precision 62 * 63 * -56.45 64 * | P - (erf(x)-x)/x | <= 2 65 * 66 * 67 * Remark. The formula is derived by noting 68 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) 69 * and that 70 * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 71 * is close to one. The interval is chosen because the fixed 72 * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is 73 * near 0.6174), and by some experiment, 0.84375 is chosen to 74 * guarantee the error is less than one ulp for erf. 75 * 76 * 3. For x in [0.84375,1.25], let s = x - 1, and 77 * c = 0.84506291151 rounded to single (24 bits) 78 * erf(x) = c + P1(s)/Q1(s) 79 * erfc(x) = (1-c) - P1(s)/Q1(s) 80 * |P1/Q1 - (erf(x)-c)| <= 2**-59.06 81 * Remark: here we use the taylor series expansion at x=1. 82 * erf(1+s) = erf(1) + s*Poly(s) 83 * = 0.845.. + P1(s)/Q1(s) 84 * That is, we use rational approximation to approximate 85 * erf(1+s) - (c = (single)0.84506291151) 86 * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] 87 * where 88 * P1(s) = degree 7 poly in s 89 * 90 * 4. For x in [1.25,6], 91 * erf(x) = 1 - erfc(x) 92 * erfc(x) = exp(-x*x)*(1/x)*R1(1/x)/S1(1/x) 93 * where 94 * R1(y) = degree 7 poly in y, (y=1/x) 95 * S1(y) = degree 8 poly in y 96 * 97 * 5. For x in [6,28] 98 * erf(x) = 1.0 - tiny 99 * erfc(x) = (1/x)exp(-x*x-(.5*log(pi)+eps) + zP(z)) 100 * 101 * Where P is degree 9 polynomial in z = 1/(x*x) 102 * 103 * Notes: 104 * Here 4 and 5 make use of the asymptotic series 105 * exp(-x*x) 106 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ); 107 * x*sqrt(pi) 108 * 109 * where for z = 1/(x*x) 110 * P(z) ~ z/2*(-1 + z*3/2*(1 + z*5/2*(-1 + z*7/2*(1 +...)))) 111 * 112 * Thus we use rational approximation to approximate 113 * erfc*x*exp(x*x) ~ 1/sqrt(pi) 114 * 115 * The error bound for the target function, G(z) for 116 * case 5 is 117 * |eps + 1/(x*x)P(1/x*x) - G(x)| < 2**(-58.34) 118 * For case 4, 119 * |R2/S2 - erfc*x*exp(x*x)| < 2**(-61.52) 120 * 121 * 6. For inf > x >= 28 122 * erf(x) = 1 - tiny (raise inexact) 123 * erfc(x) = tiny*tiny (raise underflow) 124 * 125 * 7. Special cases: 126 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, 127 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, 128 * erfc/erf(NaN) is NaN 129 */ 130 131 #if defined(vax) || defined(tahoe) 132 #define _IEEE 0 133 #define TRUNC(x) (double) (float) (x) 134 #else 135 #define _IEEE 1 136 #define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 137 #define infnan(x) 0.0 138 #endif 139 140 #ifdef _IEEE_LIBM 141 /* 142 * redefining "___function" to "function" in _IEEE_LIBM mode 143 */ 144 #include "ieee_libm.h" 145 #endif 146 147 static double 148 tiny = 1e-300, 149 half = 0.5, 150 one = 1.0, 151 two = 2.0, 152 c = 8.45062911510467529297e-01, /* (float)0.84506291151 */ 153 /* 154 * Coefficients for approximation to erf on [0,0.84375] 155 */ 156 p0t8 = 1.02703333676410051049867154944018394163280, 157 p0 = 1.283791670955125638123339436800229927041e-0001, 158 p1 = -3.761263890318340796574473028946097022260e-0001, 159 p2 = 1.128379167093567004871858633779992337238e-0001, 160 p3 = -2.686617064084433642889526516177508374437e-0002, 161 p4 = 5.223977576966219409445780927846432273191e-0003, 162 p5 = -8.548323822001639515038738961618255438422e-0004, 163 p6 = 1.205520092530505090384383082516403772317e-0004, 164 p7 = -1.492214100762529635365672665955239554276e-0005, 165 p8 = 1.640186161764254363152286358441771740838e-0006, 166 p9 = -1.571599331700515057841960987689515895479e-0007, 167 p10= 1.073087585213621540635426191486561494058e-0008, 168 /* 169 * Coefficients for approximation to erf in [0.84375,1.25] 170 */ 171 pa0 = -2.362118560752659485957248365514511540287e-0003, 172 pa1 = 4.148561186837483359654781492060070469522e-0001, 173 pa2 = -3.722078760357013107593507594535478633044e-0001, 174 pa3 = 3.183466199011617316853636418691420262160e-0001, 175 pa4 = -1.108946942823966771253985510891237782544e-0001, 176 pa5 = 3.547830432561823343969797140537411825179e-0002, 177 pa6 = -2.166375594868790886906539848893221184820e-0003, 178 qa1 = 1.064208804008442270765369280952419863524e-0001, 179 qa2 = 5.403979177021710663441167681878575087235e-0001, 180 qa3 = 7.182865441419627066207655332170665812023e-0002, 181 qa4 = 1.261712198087616469108438860983447773726e-0001, 182 qa5 = 1.363708391202905087876983523620537833157e-0002, 183 qa6 = 1.198449984679910764099772682882189711364e-0002, 184 /* 185 * Coefficients for approximation to erfc in [1.25,6] 186 */ 187 ra0 = 5.641895806197543833169694096883621225329e-0001, 188 ra1 = 7.239004794325021293310782759791744583987e+0000, 189 ra2 = 4.615482605646378370356340497765510677914e+0001, 190 ra3 = 1.831130716384318567879039478746072928548e+0002, 191 ra4 = 4.827304689401256945023566678442020977744e+0002, 192 ra5 = 8.443683805001379929687313735294340282751e+0002, 193 ra6 = 9.151771804289399937165800774604677980269e+0002, 194 ra7 = 4.884236881266866025539987843147838061930e+0002, 195 sa1 = 1.283080158932067675016971332625972882793e+0001, 196 sa2 = 8.230730944985601552133528541648529041935e+0001, 197 sa3 = 3.309746710535947168967275132570416337810e+0002, 198 sa4 = 8.960238586988354676031385802384985611536e+0002, 199 sa5 = 1.652440076836585407285764071805622271834e+0003, 200 sa6 = 2.010492426273281289533320672757992672142e+0003, 201 sa7 = 1.466304171232599681829476641652969136592e+0003, 202 sa8 = 4.884237022526160104676542187698268809111e+0002; 203 /* 204 * Coefficients for approximation to erfc in [6,28] 205 */ 206 #define a0_hi -0.5723649429247001929610405568 /* ~-.5log(pi) */ 207 #define a0_lo -0.0000000000000000189783711362898601 208 209 #define P0 -4.99999999999749700219098258458e-0001 /* -1 /2 */ 210 #define P1 6.24999999807451578348604925850e-0001 /* 5/2 /4 */ 211 #define P2 -1.54166659013994022942029005208e+0001 /* -37/3 /8 */ 212 #define P3 5.51560710872094706047619183664e+0001 /* 353/4 /16 */ 213 #define P4 -2.55036053070125880992691236315e+0002 /* -4081/5 /32 */ 214 #define P5 1.43505282730286381820405949838e+0002 /* 55205/6 /64 */ 215 #define P6 -9.36421869861889035746571607888e+0002 /* ....etc.... */ 216 #define P7 6.51030087738772090233396738768e+0003 217 #define P8 -3.98835620275180117459967732430e+0004 218 #define P9 1.44460450428346201078966259956e+0005 219 220 221 double erf(x) 222 double x; 223 { 224 double R,S,P,Q,ax,s,y,z,odd,even,r,fabs(),exp(); 225 if(!finite(x)) { /* erf(nan)=nan */ 226 if (isnan(x)) 227 return(x); 228 return (x > 0 ? one : -one); /* erf(+/-inf)= +/-1 */ 229 } 230 if ((ax = x) < 0) 231 ax = - ax; 232 if (ax < .84375) { 233 if (ax < 3.7e-09) { 234 if (ax < 1.0e-308) 235 return 0.125*(8.0*x+p0t8*x); /*avoid underflow */ 236 return x + p0*x; 237 } 238 y = x*x; 239 z = y*y; 240 even = z*(p2+z*(p4+z*(p6+z*(p8+z*p10)))); 241 odd = p1+z*(p3+z*(p5+z*(p7+z*p9))); 242 r = y*odd+even; 243 return x + x*(p0+r); 244 } 245 if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ 246 s = fabs(x)-one; 247 P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); 248 Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); 249 if (x>=0) 250 return (c + P/Q); 251 else 252 return (-c - P/Q); 253 } 254 if (ax >= 6.0) { /* inf>|x|>=6 */ 255 if (x >= 0.0) 256 return (one-tiny); 257 else 258 return (tiny-one); 259 } 260 /* 1.25 <= |x| < 6 */ 261 s = one/fabs(x); 262 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7)))))); 263 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8))))))); 264 z = exp(-x*x)*(R/S)*s; 265 if (x >= 0) 266 return (one-z); 267 else 268 return (z-one); 269 } 270 271 double erfc(x) 272 double x; 273 { 274 double R,S,P,Q,s,ax,y,odd,even,z,r,fabs(),exp__D(); 275 if (!finite(x)) { 276 if (isnan(x)) /* erfc(NaN) = NaN */ 277 return(x); 278 else if (x > 0) /* erfc(+-inf)=0,2 */ 279 return 0.0; 280 else 281 return 2.0; 282 } 283 if ((ax = x) < 0) 284 ax = -ax; 285 if (ax < .84375) { /* |x|<0.84375 */ 286 if (ax < 1.38777878078144568e-17) /* |x|<2**-56 */ 287 return one-x; 288 y = x*x; 289 z = y*y; 290 even = z*(p2+z*(p4+z*(p6+z*(p8+z*p10)))); 291 odd = p1+z*(p3+z*(p5+z*(p7+z*p9))); 292 r = y*odd+even; 293 if (ax < .0625) { /* |x|<2**-4 */ 294 return (one-(x+x*(p0+r))); 295 } else { 296 r = x*(p0+r); 297 r += (x-half); 298 return (half - r); 299 } 300 } 301 if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ 302 s = ax-one; 303 P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); 304 Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); 305 if (x>=0) { 306 z = one-c; return z - P/Q; 307 } else { 308 z = c+P/Q; return one+z; 309 } 310 } 311 if (ax >= 28) /* Out of range */ 312 if (x>0) 313 return (tiny*tiny); 314 else 315 return (two-tiny); 316 z = ax; 317 TRUNC(z); 318 y = z - ax; y *= (ax+z); 319 z *= -z; /* Here z + y = -x^2 */ 320 if (ax >= 6) { /* 6 <= ax */ 321 s = one/(-z-y); /* 1/(x*x) */ 322 R = s*(P0+s*(P1+s*(P2+s*(P3+s*(P4+ 323 s*(P5+s*(P6+s*(P7+s*(P8+s*P9))))))))); 324 y += a0_lo; 325 /* return exp(-x^2 + a0_hi + R)/x; */ 326 s = ((R + y) + a0_hi) + z; 327 y = (((z-s) + a0_hi) + R) + y; 328 r = exp__D(s, y)/x; 329 } else { /* 1.25 <= ax <= 6 */ 330 s = one/(ax); 331 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+ 332 s*(ra5+s*(ra6+s*ra7)))))); 333 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+ 334 s*(sa5+s*(sa6+s*(sa7+s*sa8))))))); 335 r = (R/S)/x; 336 s = z + y; y = (z-s) + y; 337 r *= exp__D(s, y); 338 } 339 if (x>0) 340 return r; 341 else 342 return two-r; 343 } 344