| OLD | NEW |
|---|
| (Empty) | |
| 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_erfl.c */ |
| 2 /* |
| 3 * ==================================================== |
| 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 5 * |
| 6 * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 7 * Permission to use, copy, modify, and distribute this |
| 8 * software is freely granted, provided that this notice |
| 9 * is preserved. |
| 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 /* double erf(double x) |
| 28 * double erfc(double x) |
| 29 * x |
| 30 * 2 |\ |
| 31 * erf(x) = --------- | exp(-t*t)dt |
| 32 * sqrt(pi) \| |
| 33 * 0 |
| 34 * |
| 35 * erfc(x) = 1-erf(x) |
| 36 * Note that |
| 37 * erf(-x) = -erf(x) |
| 38 * erfc(-x) = 2 - erfc(x) |
| 39 * |
| 40 * Method: |
| 41 * 1. For |x| in [0, 0.84375] |
| 42 * erf(x) = x + x*R(x^2) |
| 43 * erfc(x) = 1 - erf(x) if x in [-.84375,0.25] |
| 44 * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375] |
| 45 * Remark. The formula is derived by noting |
| 46 * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) |
| 47 * and that |
| 48 * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 |
| 49 * is close to one. The interval is chosen because the fix |
| 50 * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is |
| 51 * near 0.6174), and by some experiment, 0.84375 is chosen to |
| 52 * guarantee the error is less than one ulp for erf. |
| 53 * |
| 54 * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and |
| 55 * c = 0.84506291151 rounded to single (24 bits) |
| 56 * erf(x) = sign(x) * (c + P1(s)/Q1(s)) |
| 57 * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0 |
| 58 * 1+(c+P1(s)/Q1(s)) if x < 0 |
| 59 * Remark: here we use the taylor series expansion at x=1. |
| 60 * erf(1+s) = erf(1) + s*Poly(s) |
| 61 * = 0.845.. + P1(s)/Q1(s) |
| 62 * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] |
| 63 * |
| 64 * 3. For x in [1.25,1/0.35(~2.857143)], |
| 65 * erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z)) |
| 66 * z=1/x^2 |
| 67 * erf(x) = 1 - erfc(x) |
| 68 * |
| 69 * 4. For x in [1/0.35,107] |
| 70 * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0 |
| 71 * = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z)) |
| 72 * if -6.666<x<0 |
| 73 * = 2.0 - tiny (if x <= -6.666) |
| 74 * z=1/x^2 |
| 75 * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6.666, else |
| 76 * erf(x) = sign(x)*(1.0 - tiny) |
| 77 * Note1: |
| 78 * To compute exp(-x*x-0.5625+R/S), let s be a single |
| 79 * precision number and s := x; then |
| 80 * -x*x = -s*s + (s-x)*(s+x) |
| 81 * exp(-x*x-0.5626+R/S) = |
| 82 * exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S); |
| 83 * Note2: |
| 84 * Here 4 and 5 make use of the asymptotic series |
| 85 * exp(-x*x) |
| 86 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ) |
| 87 * x*sqrt(pi) |
| 88 * |
| 89 * 5. For inf > x >= 107 |
| 90 * erf(x) = sign(x) *(1 - tiny) (raise inexact) |
| 91 * erfc(x) = tiny*tiny (raise underflow) if x > 0 |
| 92 * = 2 - tiny if x<0 |
| 93 * |
| 94 * 7. Special case: |
| 95 * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, |
| 96 * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, |
| 97 * erfc/erf(NaN) is NaN |
| 98 */ |
| 99 |
| 100 |
| 101 #include "libm.h" |
| 102 |
| 103 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 104 long double erfl(long double x) |
| 105 { |
| 106 return erf(x); |
| 107 } |
| 108 long double erfcl(long double x) |
| 109 { |
| 110 return erfc(x); |
| 111 } |
| 112 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
| 113 static const long double |
| 114 erx = 0.845062911510467529296875L, |
| 115 |
| 116 /* |
| 117 * Coefficients for approximation to erf on [0,0.84375] |
| 118 */ |
| 119 /* 8 * (2/sqrt(pi) - 1) */ |
| 120 efx8 = 1.0270333367641005911692712249723613735048E0L, |
| 121 pp[6] = { |
| 122 1.122751350964552113068262337278335028553E6L, |
| 123 -2.808533301997696164408397079650699163276E6L, |
| 124 -3.314325479115357458197119660818768924100E5L, |
| 125 -6.848684465326256109712135497895525446398E4L, |
| 126 -2.657817695110739185591505062971929859314E3L, |
| 127 -1.655310302737837556654146291646499062882E2L, |
| 128 }, |
| 129 qq[6] = { |
| 130 8.745588372054466262548908189000448124232E6L, |
| 131 3.746038264792471129367533128637019611485E6L, |
| 132 7.066358783162407559861156173539693900031E5L, |
| 133 7.448928604824620999413120955705448117056E4L, |
| 134 4.511583986730994111992253980546131408924E3L, |
| 135 1.368902937933296323345610240009071254014E2L, |
| 136 /* 1.000000000000000000000000000000000000000E0 */ |
| 137 }, |
| 138 |
| 139 /* |
| 140 * Coefficients for approximation to erf in [0.84375,1.25] |
| 141 */ |
| 142 /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x) |
| 143 -0.15625 <= x <= +.25 |
| 144 Peak relative error 8.5e-22 */ |
| 145 pa[8] = { |
| 146 -1.076952146179812072156734957705102256059E0L, |
| 147 1.884814957770385593365179835059971587220E2L, |
| 148 -5.339153975012804282890066622962070115606E1L, |
| 149 4.435910679869176625928504532109635632618E1L, |
| 150 1.683219516032328828278557309642929135179E1L, |
| 151 -2.360236618396952560064259585299045804293E0L, |
| 152 1.852230047861891953244413872297940938041E0L, |
| 153 9.394994446747752308256773044667843200719E-2L, |
| 154 }, |
| 155 qa[7] = { |
| 156 4.559263722294508998149925774781887811255E2L, |
| 157 3.289248982200800575749795055149780689738E2L, |
| 158 2.846070965875643009598627918383314457912E2L, |
| 159 1.398715859064535039433275722017479994465E2L, |
| 160 6.060190733759793706299079050985358190726E1L, |
| 161 2.078695677795422351040502569964299664233E1L, |
| 162 4.641271134150895940966798357442234498546E0L, |
| 163 /* 1.000000000000000000000000000000000000000E0 */ |
| 164 }, |
| 165 |
| 166 /* |
| 167 * Coefficients for approximation to erfc in [1.25,1/0.35] |
| 168 */ |
| 169 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2)) |
| 170 1/2.85711669921875 < 1/x < 1/1.25 |
| 171 Peak relative error 3.1e-21 */ |
| 172 ra[] = { |
| 173 1.363566591833846324191000679620738857234E-1L, |
| 174 1.018203167219873573808450274314658434507E1L, |
| 175 1.862359362334248675526472871224778045594E2L, |
| 176 1.411622588180721285284945138667933330348E3L, |
| 177 5.088538459741511988784440103218342840478E3L, |
| 178 8.928251553922176506858267311750789273656E3L, |
| 179 7.264436000148052545243018622742770549982E3L, |
| 180 2.387492459664548651671894725748959751119E3L, |
| 181 2.220916652813908085449221282808458466556E2L, |
| 182 }, |
| 183 sa[] = { |
| 184 -1.382234625202480685182526402169222331847E1L, |
| 185 -3.315638835627950255832519203687435946482E2L, |
| 186 -2.949124863912936259747237164260785326692E3L, |
| 187 -1.246622099070875940506391433635999693661E4L, |
| 188 -2.673079795851665428695842853070996219632E4L, |
| 189 -2.880269786660559337358397106518918220991E4L, |
| 190 -1.450600228493968044773354186390390823713E4L, |
| 191 -2.874539731125893533960680525192064277816E3L, |
| 192 -1.402241261419067750237395034116942296027E2L, |
| 193 /* 1.000000000000000000000000000000000000000E0 */ |
| 194 }, |
| 195 |
| 196 /* |
| 197 * Coefficients for approximation to erfc in [1/.35,107] |
| 198 */ |
| 199 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2)) |
| 200 1/6.6666259765625 < 1/x < 1/2.85711669921875 |
| 201 Peak relative error 4.2e-22 */ |
| 202 rb[] = { |
| 203 -4.869587348270494309550558460786501252369E-5L, |
| 204 -4.030199390527997378549161722412466959403E-3L, |
| 205 -9.434425866377037610206443566288917589122E-2L, |
| 206 -9.319032754357658601200655161585539404155E-1L, |
| 207 -4.273788174307459947350256581445442062291E0L, |
| 208 -8.842289940696150508373541814064198259278E0L, |
| 209 -7.069215249419887403187988144752613025255E0L, |
| 210 -1.401228723639514787920274427443330704764E0L, |
| 211 }, |
| 212 sb[] = { |
| 213 4.936254964107175160157544545879293019085E-3L, |
| 214 1.583457624037795744377163924895349412015E-1L, |
| 215 1.850647991850328356622940552450636420484E0L, |
| 216 9.927611557279019463768050710008450625415E0L, |
| 217 2.531667257649436709617165336779212114570E1L, |
| 218 2.869752886406743386458304052862814690045E1L, |
| 219 1.182059497870819562441683560749192539345E1L, |
| 220 /* 1.000000000000000000000000000000000000000E0 */ |
| 221 }, |
| 222 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2)) |
| 223 1/107 <= 1/x <= 1/6.6666259765625 |
| 224 Peak relative error 1.1e-21 */ |
| 225 rc[] = { |
| 226 -8.299617545269701963973537248996670806850E-5L, |
| 227 -6.243845685115818513578933902532056244108E-3L, |
| 228 -1.141667210620380223113693474478394397230E-1L, |
| 229 -7.521343797212024245375240432734425789409E-1L, |
| 230 -1.765321928311155824664963633786967602934E0L, |
| 231 -1.029403473103215800456761180695263439188E0L, |
| 232 }, |
| 233 sc[] = { |
| 234 8.413244363014929493035952542677768808601E-3L, |
| 235 2.065114333816877479753334599639158060979E-1L, |
| 236 1.639064941530797583766364412782135680148E0L, |
| 237 4.936788463787115555582319302981666347450E0L, |
| 238 5.005177727208955487404729933261347679090E0L, |
| 239 /* 1.000000000000000000000000000000000000000E0 */ |
| 240 }; |
| 241 |
| 242 static long double erfc1(long double x) |
| 243 { |
| 244 long double s,P,Q; |
| 245 |
| 246 s = fabsl(x) - 1; |
| 247 P = pa[0] + s * (pa[1] + s * (pa[2] + |
| 248 s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))))
; |
| 249 Q = qa[0] + s * (qa[1] + s * (qa[2] + |
| 250 s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s)))))); |
| 251 return 1 - erx - P / Q; |
| 252 } |
| 253 |
| 254 static long double erfc2(uint32_t ix, long double x) |
| 255 { |
| 256 union ldshape u; |
| 257 long double s,z,R,S; |
| 258 |
| 259 if (ix < 0x3fffa000) /* 0.84375 <= |x| < 1.25 */ |
| 260 return erfc1(x); |
| 261 |
| 262 x = fabsl(x); |
| 263 s = 1 / (x * x); |
| 264 if (ix < 0x4000b6db) { /* 1.25 <= |x| < 2.857 ~ 1/.35 */ |
| 265 R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] + |
| 266 s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8]))))))); |
| 267 S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] + |
| 268 s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))
)); |
| 269 } else if (ix < 0x4001d555) { /* 2.857 <= |x| < 6.6666259765625 */ |
| 270 R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] + |
| 271 s * (rb[5] + s * (rb[6] + s * rb[7])))))); |
| 272 S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] + |
| 273 s * (sb[5] + s * (sb[6] + s)))))); |
| 274 } else { /* 6.666 <= |x| < 107 (erfc only) */ |
| 275 R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] + |
| 276 s * (rc[4] + s * rc[5])))); |
| 277 S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] + |
| 278 s * (sc[4] + s)))); |
| 279 } |
| 280 u.f = x; |
| 281 u.i.m &= -1ULL << 40; |
| 282 z = u.f; |
| 283 return expl(-z*z - 0.5625) * expl((z - x) * (z + x) + R / S) / x; |
| 284 } |
| 285 |
| 286 long double erfl(long double x) |
| 287 { |
| 288 long double r, s, z, y; |
| 289 union ldshape u = {x}; |
| 290 uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48; |
| 291 int sign = u.i.se >> 15; |
| 292 |
| 293 if (ix >= 0x7fff0000) |
| 294 /* erf(nan)=nan, erf(+-inf)=+-1 */ |
| 295 return 1 - 2*sign + 1/x; |
| 296 if (ix < 0x3ffed800) { /* |x| < 0.84375 */ |
| 297 if (ix < 0x3fde8000) { /* |x| < 2**-33 */ |
| 298 return 0.125 * (8 * x + efx8 * x); /* avoid underflow *
/ |
| 299 } |
| 300 z = x * x; |
| 301 r = pp[0] + z * (pp[1] + |
| 302 z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5])))); |
| 303 s = qq[0] + z * (qq[1] + |
| 304 z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z))))); |
| 305 y = r / s; |
| 306 return x + x * y; |
| 307 } |
| 308 if (ix < 0x4001d555) /* |x| < 6.6666259765625 */ |
| 309 y = 1 - erfc2(ix,x); |
| 310 else |
| 311 y = 1 - 0x1p-16382L; |
| 312 return sign ? -y : y; |
| 313 } |
| 314 |
| 315 long double erfcl(long double x) |
| 316 { |
| 317 long double r, s, z, y; |
| 318 union ldshape u = {x}; |
| 319 uint32_t ix = (u.i.se & 0x7fffU)<<16 | u.i.m>>48; |
| 320 int sign = u.i.se >> 15; |
| 321 |
| 322 if (ix >= 0x7fff0000) |
| 323 /* erfc(nan) = nan, erfc(+-inf) = 0,2 */ |
| 324 return 2*sign + 1/x; |
| 325 if (ix < 0x3ffed800) { /* |x| < 0.84375 */ |
| 326 if (ix < 0x3fbe0000) /* |x| < 2**-65 */ |
| 327 return 1.0 - x; |
| 328 z = x * x; |
| 329 r = pp[0] + z * (pp[1] + |
| 330 z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5])))); |
| 331 s = qq[0] + z * (qq[1] + |
| 332 z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z))))); |
| 333 y = r / s; |
| 334 if (ix < 0x3ffd8000) /* x < 1/4 */ |
| 335 return 1.0 - (x + x * y); |
| 336 return 0.5 - (x - 0.5 + x * y); |
| 337 } |
| 338 if (ix < 0x4005d600) /* |x| < 107 */ |
| 339 return sign ? 2 - erfc2(ix,x) : erfc2(ix,x); |
| 340 y = 0x1p-16382L; |
| 341 return sign ? 2 - y : y*y; |
| 342 } |
| 343 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
| 344 // TODO: broken implementation to make things compile |
| 345 long double erfl(long double x) |
| 346 { |
| 347 return erf(x); |
| 348 } |
| 349 long double erfcl(long double x) |
| 350 { |
| 351 return erfc(x); |
| 352 } |
| 353 #endif |
| OLD | NEW |
|---|
Chromium Code Reviews
Issue 1573973002:
Add a "fork" of musl as //fusl. (Closed)
Base URL: https://github.com/domokit/mojo.git@master