[352] | 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
|
---|