1 /* 2 * ==================================================== 3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * 5 * Developed at SunPro, a Sun Microsystems, Inc. business. 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 /* 13 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 14 * 15 * Permission to use, copy, modify, and distribute this software for any 16 * purpose with or without fee is hereby granted, provided that the above 17 * copyright notice and this permission notice appear in all copies. 18 * 19 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 20 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 21 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 22 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 23 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 24 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 25 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 26 */ 27 28 /* double erf(double x) 29 * double erfc(double x) 30 * x 31 * 2 |\ 32 * erf(x) = --------- | exp(-t*t)dt 33 * sqrt(pi) \| 34 * 0 35 * 36 * erfc(x) = 1-erf(x) 37 * Note that 38 * erf(-x) = -erf(x) 39 * erfc(-x) = 2 - erfc(x) 40 * 41 * Method: 42 * 1. For |x| in [0, 0.84375] 43 * erf(x) = x + x*R(x^2) 44 * erfc(x) = 1 - erf(x) if x in [-.84375,0.25] 45 * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375] 46 * Remark. The formula is derived by noting 47 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) 48 * and that 49 * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 50 * is close to one. The interval is chosen because the fix 51 * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is 52 * near 0.6174), and by some experiment, 0.84375 is chosen to 53 * guarantee the error is less than one ulp for erf. 54 * 55 * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and 56 * c = 0.84506291151 rounded to single (24 bits) 57 * erf(x) = sign(x) * (c + P1(s)/Q1(s)) 58 * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0 59 * 1+(c+P1(s)/Q1(s)) if x < 0 60 * Remark: here we use the taylor series expansion at x=1. 61 * erf(1+s) = erf(1) + s*Poly(s) 62 * = 0.845.. + P1(s)/Q1(s) 63 * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] 64 * 65 * 3. For x in [1.25,1/0.35(~2.857143)], 66 * erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z)) 67 * z=1/x^2 68 * erf(x) = 1 - erfc(x) 69 * 70 * 4. For x in [1/0.35,107] 71 * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0 72 * = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z)) 73 * if -6.666<x<0 74 * = 2.0 - tiny (if x <= -6.666) 75 * z=1/x^2 76 * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else 77 * erf(x) = sign(x)*(1.0 - tiny) 78 * Note1: 79 * To compute exp(-x*x-0.5625+R/S), let s be a single 80 * precision number and s := x; then 81 * -x*x = -s*s + (s-x)*(s+x) 82 * exp(-x*x-0.5626+R/S) = 83 * exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S); 84 * Note2: 85 * Here 4 and 5 make use of the asymptotic series 86 * exp(-x*x) 87 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ) 88 * x*sqrt(pi) 89 * 90 * 5. For inf > x >= 107 91 * erf(x) = sign(x) *(1 - tiny) (raise inexact) 92 * erfc(x) = tiny*tiny (raise underflow) if x > 0 93 * = 2 - tiny if x<0 94 * 95 * 7. Special case: 96 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, 97 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, 98 * erfc/erf(NaN) is NaN 99 */ 100 101 102 #include <math.h> 103 104 #include "math_private.h" 105 106 static const long double 107 tiny = 1e-4931L, 108 half = 0.5L, 109 one = 1.0L, 110 two = 2.0L, 111 /* c = (float)0.84506291151 */ 112 erx = 0.845062911510467529296875L, 113 /* 114 * Coefficients for approximation to erf on [0,0.84375] 115 */ 116 /* 2/sqrt(pi) - 1 */ 117 efx = 1.2837916709551257389615890312154517168810E-1L, 118 /* 8 * (2/sqrt(pi) - 1) */ 119 efx8 = 1.0270333367641005911692712249723613735048E0L, 120 121 pp[6] = { 122 1.122751350964552113068262337278335028553E6L, 123 -2.808533301997696164408397079650699163276E6L, 124 -3.314325479115357458197119660818768924100E5L, 125 -6.848684465326256109712135497895525446398E4L, 126 -2.657817695110739185591505062971929859314E3L, 127 -1.655310302737837556654146291646499062882E2L, 128 }, 129 130 qq[6] = { 131 8.745588372054466262548908189000448124232E6L, 132 3.746038264792471129367533128637019611485E6L, 133 7.066358783162407559861156173539693900031E5L, 134 7.448928604824620999413120955705448117056E4L, 135 4.511583986730994111992253980546131408924E3L, 136 1.368902937933296323345610240009071254014E2L, 137 /* 1.000000000000000000000000000000000000000E0 */ 138 }, 139 140 /* 141 * Coefficients for approximation to erf in [0.84375,1.25] 142 */ 143 /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x) 144 -0.15625 <= x <= +.25 145 Peak relative error 8.5e-22 */ 146 147 pa[8] = { 148 -1.076952146179812072156734957705102256059E0L, 149 1.884814957770385593365179835059971587220E2L, 150 -5.339153975012804282890066622962070115606E1L, 151 4.435910679869176625928504532109635632618E1L, 152 1.683219516032328828278557309642929135179E1L, 153 -2.360236618396952560064259585299045804293E0L, 154 1.852230047861891953244413872297940938041E0L, 155 9.394994446747752308256773044667843200719E-2L, 156 }, 157 158 qa[7] = { 159 4.559263722294508998149925774781887811255E2L, 160 3.289248982200800575749795055149780689738E2L, 161 2.846070965875643009598627918383314457912E2L, 162 1.398715859064535039433275722017479994465E2L, 163 6.060190733759793706299079050985358190726E1L, 164 2.078695677795422351040502569964299664233E1L, 165 4.641271134150895940966798357442234498546E0L, 166 /* 1.000000000000000000000000000000000000000E0 */ 167 }, 168 169 /* 170 * Coefficients for approximation to erfc in [1.25,1/0.35] 171 */ 172 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2)) 173 1/2.85711669921875 < 1/x < 1/1.25 174 Peak relative error 3.1e-21 */ 175 176 ra[] = { 177 1.363566591833846324191000679620738857234E-1L, 178 1.018203167219873573808450274314658434507E1L, 179 1.862359362334248675526472871224778045594E2L, 180 1.411622588180721285284945138667933330348E3L, 181 5.088538459741511988784440103218342840478E3L, 182 8.928251553922176506858267311750789273656E3L, 183 7.264436000148052545243018622742770549982E3L, 184 2.387492459664548651671894725748959751119E3L, 185 2.220916652813908085449221282808458466556E2L, 186 }, 187 188 sa[] = { 189 -1.382234625202480685182526402169222331847E1L, 190 -3.315638835627950255832519203687435946482E2L, 191 -2.949124863912936259747237164260785326692E3L, 192 -1.246622099070875940506391433635999693661E4L, 193 -2.673079795851665428695842853070996219632E4L, 194 -2.880269786660559337358397106518918220991E4L, 195 -1.450600228493968044773354186390390823713E4L, 196 -2.874539731125893533960680525192064277816E3L, 197 -1.402241261419067750237395034116942296027E2L, 198 /* 1.000000000000000000000000000000000000000E0 */ 199 }, 200 /* 201 * Coefficients for approximation to erfc in [1/.35,107] 202 */ 203 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2)) 204 1/6.6666259765625 < 1/x < 1/2.85711669921875 205 Peak relative error 4.2e-22 */ 206 rb[] = { 207 -4.869587348270494309550558460786501252369E-5L, 208 -4.030199390527997378549161722412466959403E-3L, 209 -9.434425866377037610206443566288917589122E-2L, 210 -9.319032754357658601200655161585539404155E-1L, 211 -4.273788174307459947350256581445442062291E0L, 212 -8.842289940696150508373541814064198259278E0L, 213 -7.069215249419887403187988144752613025255E0L, 214 -1.401228723639514787920274427443330704764E0L, 215 }, 216 217 sb[] = { 218 4.936254964107175160157544545879293019085E-3L, 219 1.583457624037795744377163924895349412015E-1L, 220 1.850647991850328356622940552450636420484E0L, 221 9.927611557279019463768050710008450625415E0L, 222 2.531667257649436709617165336779212114570E1L, 223 2.869752886406743386458304052862814690045E1L, 224 1.182059497870819562441683560749192539345E1L, 225 /* 1.000000000000000000000000000000000000000E0 */ 226 }, 227 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2)) 228 1/107 <= 1/x <= 1/6.6666259765625 229 Peak relative error 1.1e-21 */ 230 rc[] = { 231 -8.299617545269701963973537248996670806850E-5L, 232 -6.243845685115818513578933902532056244108E-3L, 233 -1.141667210620380223113693474478394397230E-1L, 234 -7.521343797212024245375240432734425789409E-1L, 235 -1.765321928311155824664963633786967602934E0L, 236 -1.029403473103215800456761180695263439188E0L, 237 }, 238 239 sc[] = { 240 8.413244363014929493035952542677768808601E-3L, 241 2.065114333816877479753334599639158060979E-1L, 242 1.639064941530797583766364412782135680148E0L, 243 4.936788463787115555582319302981666347450E0L, 244 5.005177727208955487404729933261347679090E0L, 245 /* 1.000000000000000000000000000000000000000E0 */ 246 }; 247 248 long double 249 erfl(long double x) 250 { 251 long double R, S, P, Q, s, y, z, r; 252 int32_t ix, i; 253 u_int32_t se, i0, i1; 254 255 GET_LDOUBLE_WORDS (se, i0, i1, x); 256 ix = se & 0x7fff; 257 258 if (ix >= 0x7fff) 259 { /* erf(nan)=nan */ 260 i = ((se & 0xffff) >> 15) << 1; 261 return (long double) (1 - i) + one / x; /* erf(+-inf)=+-1 */ 262 } 263 264 ix = (ix << 16) | (i0 >> 16); 265 if (ix < 0x3ffed800) /* |x|<0.84375 */ 266 { 267 if (ix < 0x3fde8000) /* |x|<2**-33 */ 268 { 269 if (ix < 0x00080000) 270 return 0.125 * (8.0 * x + efx8 * x); /*avoid underflow */ 271 return x + efx * x; 272 } 273 z = x * x; 274 r = pp[0] + z * (pp[1] 275 + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5])))); 276 s = qq[0] + z * (qq[1] 277 + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z))))); 278 y = r / s; 279 return x + x * y; 280 } 281 if (ix < 0x3fffa000) /* 1.25 */ 282 { /* 0.84375 <= |x| < 1.25 */ 283 s = fabsl (x) - one; 284 P = pa[0] + s * (pa[1] + s * (pa[2] 285 + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); 286 Q = qa[0] + s * (qa[1] + s * (qa[2] 287 + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s)))))); 288 if ((se & 0x8000) == 0) 289 return erx + P / Q; 290 else 291 return -erx - P / Q; 292 } 293 if (ix >= 0x4001d555) /* 6.6666259765625 */ 294 { /* inf>|x|>=6.666 */ 295 if ((se & 0x8000) == 0) 296 return one - tiny; 297 else 298 return tiny - one; 299 } 300 x = fabsl (x); 301 s = one / (x * x); 302 if (ix < 0x4000b6db) /* 2.85711669921875 */ 303 { 304 R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] + 305 s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8]))))))); 306 S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] + 307 s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s)))))))); 308 } 309 else 310 { /* |x| >= 1/0.35 */ 311 R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] + 312 s * (rb[5] + s * (rb[6] + s * rb[7])))))); 313 S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] + 314 s * (sb[5] + s * (sb[6] + s)))))); 315 } 316 z = x; 317 GET_LDOUBLE_WORDS (i, i0, i1, z); 318 i1 = 0; 319 SET_LDOUBLE_WORDS (z, i, i0, i1); 320 r = 321 expl (-z * z - 0.5625) * expl ((z - x) * (z + x) + R / S); 322 if ((se & 0x8000) == 0) 323 return one - r / x; 324 else 325 return r / x - one; 326 } 327 328 long double 329 erfcl(long double x) 330 { 331 int32_t hx, ix; 332 long double R, S, P, Q, s, y, z, r; 333 u_int32_t se, i0, i1; 334 335 GET_LDOUBLE_WORDS (se, i0, i1, x); 336 ix = se & 0x7fff; 337 if (ix >= 0x7fff) 338 { /* erfc(nan)=nan */ 339 /* erfc(+-inf)=0,2 */ 340 return (long double) (((se & 0xffff) >> 15) << 1) + one / x; 341 } 342 343 ix = (ix << 16) | (i0 >> 16); 344 if (ix < 0x3ffed800) /* |x|<0.84375 */ 345 { 346 if (ix < 0x3fbe0000) /* |x|<2**-65 */ 347 return one - x; 348 z = x * x; 349 r = pp[0] + z * (pp[1] 350 + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5])))); 351 s = qq[0] + z * (qq[1] 352 + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z))))); 353 y = r / s; 354 if (ix < 0x3ffd8000) /* x<1/4 */ 355 { 356 return one - (x + x * y); 357 } 358 else 359 { 360 r = x * y; 361 r += (x - half); 362 return half - r; 363 } 364 } 365 if (ix < 0x3fffa000) /* 1.25 */ 366 { /* 0.84375 <= |x| < 1.25 */ 367 s = fabsl (x) - one; 368 P = pa[0] + s * (pa[1] + s * (pa[2] 369 + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7])))))); 370 Q = qa[0] + s * (qa[1] + s * (qa[2] 371 + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s)))))); 372 if ((se & 0x8000) == 0) 373 { 374 z = one - erx; 375 return z - P / Q; 376 } 377 else 378 { 379 z = erx + P / Q; 380 return one + z; 381 } 382 } 383 if (ix < 0x4005d600) /* 107 */ 384 { /* |x|<107 */ 385 x = fabsl (x); 386 s = one / (x * x); 387 if (ix < 0x4000b6db) /* 2.85711669921875 */ 388 { /* |x| < 1/.35 ~ 2.857143 */ 389 R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] + 390 s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8]))))))); 391 S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] + 392 s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s)))))))); 393 } 394 else if (ix < 0x4001d555) /* 6.6666259765625 */ 395 { /* 6.666 > |x| >= 1/.35 ~ 2.857143 */ 396 R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] + 397 s * (rb[5] + s * (rb[6] + s * rb[7])))))); 398 S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] + 399 s * (sb[5] + s * (sb[6] + s)))))); 400 } 401 else 402 { /* |x| >= 6.666 */ 403 if (se & 0x8000) 404 return two - tiny; /* x < -6.666 */ 405 406 R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] + 407 s * (rc[4] + s * rc[5])))); 408 S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] + 409 s * (sc[4] + s)))); 410 } 411 z = x; 412 GET_LDOUBLE_WORDS (hx, i0, i1, z); 413 i1 = 0; 414 i0 &= 0xffffff00; 415 SET_LDOUBLE_WORDS (z, hx, i0, i1); 416 r = expl (-z * z - 0.5625) * 417 expl ((z - x) * (z + x) + R / S); 418 if ((se & 0x8000) == 0) 419 return r / x; 420 else 421 return two - r / x; 422 } 423 else 424 { 425 if ((se & 0x8000) == 0) 426 return tiny * tiny; 427 else 428 return two - tiny; 429 } 430 } 431