nsz Git - libm/blob - src/math/erfl.c

   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 LD64
 104 long double erfl(long double x)
 105 {
 106         return erfl(x);
 107 }
 108 #elif LD80
 109 static const long double
 110 tiny = 1e-4931L,
 111 half = 0.5L,
 112 one = 1.0L,
 113 two = 2.0L,
 114 /* c = (float)0.84506291151 */
 115 erx = 0.845062911510467529296875L,
 116
 117 /*
 118  * Coefficients for approximation to  erf on [0,0.84375]
 119  */
 120 /* 2/sqrt(pi) - 1 */
 121 efx = 1.2837916709551257389615890312154517168810E-1L,
 122 /* 8 * (2/sqrt(pi) - 1) */
 123 efx8 = 1.0270333367641005911692712249723613735048E0L,
 124 pp[6] = {
 125         1.122751350964552113068262337278335028553E6L,
 126         -2.808533301997696164408397079650699163276E6L,
 127         -3.314325479115357458197119660818768924100E5L,
 128         -6.848684465326256109712135497895525446398E4L,
 129         -2.657817695110739185591505062971929859314E3L,
 130         -1.655310302737837556654146291646499062882E2L,
 131 },
 132 qq[6] = {
 133         8.745588372054466262548908189000448124232E6L,
 134         3.746038264792471129367533128637019611485E6L,
 135         7.066358783162407559861156173539693900031E5L,
 136         7.448928604824620999413120955705448117056E4L,
 137         4.511583986730994111992253980546131408924E3L,
 138         1.368902937933296323345610240009071254014E2L,
 139         /* 1.000000000000000000000000000000000000000E0 */
 140 },
 141
 142 /*
 143  * Coefficients for approximation to  erf  in [0.84375,1.25]
 144  */
 145 /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
 146    -0.15625 <= x <= +.25
 147    Peak relative error 8.5e-22  */
 148 pa[8] = {
 149         -1.076952146179812072156734957705102256059E0L,
 150          1.884814957770385593365179835059971587220E2L,
 151         -5.339153975012804282890066622962070115606E1L,
 152          4.435910679869176625928504532109635632618E1L,
 153          1.683219516032328828278557309642929135179E1L,
 154         -2.360236618396952560064259585299045804293E0L,
 155          1.852230047861891953244413872297940938041E0L,
 156          9.394994446747752308256773044667843200719E-2L,
 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 ra[] = {
 176         1.363566591833846324191000679620738857234E-1L,
 177         1.018203167219873573808450274314658434507E1L,
 178         1.862359362334248675526472871224778045594E2L,
 179         1.411622588180721285284945138667933330348E3L,
 180         5.088538459741511988784440103218342840478E3L,
 181         8.928251553922176506858267311750789273656E3L,
 182         7.264436000148052545243018622742770549982E3L,
 183         2.387492459664548651671894725748959751119E3L,
 184         2.220916652813908085449221282808458466556E2L,
 185 },
 186 sa[] = {
 187         -1.382234625202480685182526402169222331847E1L,
 188         -3.315638835627950255832519203687435946482E2L,
 189         -2.949124863912936259747237164260785326692E3L,
 190         -1.246622099070875940506391433635999693661E4L,
 191         -2.673079795851665428695842853070996219632E4L,
 192         -2.880269786660559337358397106518918220991E4L,
 193         -1.450600228493968044773354186390390823713E4L,
 194         -2.874539731125893533960680525192064277816E3L,
 195         -1.402241261419067750237395034116942296027E2L,
 196         /* 1.000000000000000000000000000000000000000E0 */
 197 },
 198
 199 /*
 200  * Coefficients for approximation to  erfc in [1/.35,107]
 201  */
 202 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
 203    1/6.6666259765625 < 1/x < 1/2.85711669921875
 204    Peak relative error 4.2e-22  */
 205 rb[] = {
 206         -4.869587348270494309550558460786501252369E-5L,
 207         -4.030199390527997378549161722412466959403E-3L,
 208         -9.434425866377037610206443566288917589122E-2L,
 209         -9.319032754357658601200655161585539404155E-1L,
 210         -4.273788174307459947350256581445442062291E0L,
 211         -8.842289940696150508373541814064198259278E0L,
 212         -7.069215249419887403187988144752613025255E0L,
 213         -1.401228723639514787920274427443330704764E0L,
 214 },
 215 sb[] = {
 216         4.936254964107175160157544545879293019085E-3L,
 217         1.583457624037795744377163924895349412015E-1L,
 218         1.850647991850328356622940552450636420484E0L,
 219         9.927611557279019463768050710008450625415E0L,
 220         2.531667257649436709617165336779212114570E1L,
 221         2.869752886406743386458304052862814690045E1L,
 222         1.182059497870819562441683560749192539345E1L,
 223         /* 1.000000000000000000000000000000000000000E0 */
 224 },
 225 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
 226    1/107 <= 1/x <= 1/6.6666259765625
 227    Peak relative error 1.1e-21  */
 228 rc[] = {
 229         -8.299617545269701963973537248996670806850E-5L,
 230         -6.243845685115818513578933902532056244108E-3L,
 231         -1.141667210620380223113693474478394397230E-1L,
 232         -7.521343797212024245375240432734425789409E-1L,
 233         -1.765321928311155824664963633786967602934E0L,
 234         -1.029403473103215800456761180695263439188E0L,
 235 },
 236 sc[] = {
 237         8.413244363014929493035952542677768808601E-3L,
 238         2.065114333816877479753334599639158060979E-1L,
 239         1.639064941530797583766364412782135680148E0L,
 240         4.936788463787115555582319302981666347450E0L,
 241         5.005177727208955487404729933261347679090E0L,
 242         /* 1.000000000000000000000000000000000000000E0 */
 243 };
 244
 245 long double erfl(long double x)
 246 {
 247         long double R, S, P, Q, s, y, z, r;
 248         int32_t ix, i;
 249         uint32_t se, i0, i1;
 250
 251         GET_LDOUBLE_WORDS (se, i0, i1, x);
 252         ix = se & 0x7fff;
 253
 254         if (ix >= 0x7fff) {  /* erf(nan)=nan */
 255                 i = ((se & 0xffff) >> 15) << 1;
 256                 return (long double)(1 - i) + one / x;  /* erf(+-inf)=+-1 */
 257         }
 258
 259         ix = (ix << 16) | (i0 >> 16);
 260         if (ix < 0x3ffed800) {  /* |x| < 0.84375 */
 261                 if (ix < 0x3fde8000) {  /* |x| < 2**-33 */
 262                         if (ix < 0x00080000)
 263                                 return 0.125 * (8.0 * x + efx8 * x);  /* avoid underflow */
 264                         return x + efx * x;
 265                 }
 266                 z = x * x;
 267                 r = pp[0] + z * (pp[1] +
 268                      z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
 269                 s = qq[0] + z * (qq[1] +
 270                      z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
 271                 y = r / s;
 272                 return x + x * y;
 273         }
 274         if (ix < 0x3fffa000) {  /* 0.84375 <= |x| < 1.25 */
 275                 s = fabsl (x) - one;
 276                 P = pa[0] + s * (pa[1] + s * (pa[2] +
 277                      s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
 278                 Q = qa[0] + s * (qa[1] + s * (qa[2] +
 279                      s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
 280                 if ((se & 0x8000) == 0)
 281                         return erx + P / Q;
 282                 return -erx - P / Q;
 283         }
 284         if (ix >= 0x4001d555) {  /* inf > |x| >= 6.6666259765625 */
 285                 if ((se & 0x8000) == 0)
 286                         return one - tiny;
 287                 return tiny - one;
 288         }
 289         x = fabsl (x);
 290         s = one / (x * x);
 291         if (ix < 0x4000b6db) {  /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
 292                 R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
 293                      s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
 294                 S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
 295                      s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
 296         } else { /* 2.857 <= |x| < 6.667 */
 297                 R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
 298                      s * (rb[5] + s * (rb[6] + s * rb[7]))))));
 299                 S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
 300                      s * (sb[5] + s * (sb[6] + s))))));
 301         }
 302         z = x;
 303         GET_LDOUBLE_WORDS(i, i0, i1, z);
 304         i1 = 0;
 305         SET_LDOUBLE_WORDS(z, i, i0, i1);
 306         r = expl(-z * z - 0.5625) * expl((z - x) * (z + x) + R / S);
 307         if ((se & 0x8000) == 0)
 308                 return one - r / x;
 309         return r / x - one;
 310 }
 311
 312 long double erfcl(long double x)
 313 {
 314         int32_t hx, ix;
 315         long double R, S, P, Q, s, y, z, r;
 316         uint32_t se, i0, i1;
 317
 318         GET_LDOUBLE_WORDS (se, i0, i1, x);
 319         ix = se & 0x7fff;
 320         if (ix >= 0x7fff) {  /* erfc(nan) = nan, erfc(+-inf) = 0,2 */
 321                 return (long double)(((se & 0xffff) >> 15) << 1) + one / x;
 322         }
 323
 324         ix = (ix << 16) | (i0 >> 16);
 325         if (ix < 0x3ffed800) {  /* |x| < 0.84375 */
 326                 if (ix < 0x3fbe0000)  /* |x| < 2**-65 */
 327                         return one - 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 one - (x + x * y);
 336                 r = x * y;
 337                 r += x - half;
 338                 return half - r;
 339         }
 340         if (ix < 0x3fffa000) {  /* 0.84375 <= |x| < 1.25 */
 341                 s = fabsl (x) - one;
 342                 P = pa[0] + s * (pa[1] + s * (pa[2] +
 343                      s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
 344                 Q = qa[0] + s * (qa[1] + s * (qa[2] +
 345                      s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
 346                 if ((se & 0x8000) == 0) {
 347                         z = one - erx;
 348                         return z - P / Q;
 349                 }
 350                 z = erx + P / Q;
 351                 return one + z;
 352         }
 353         if (ix < 0x4005d600) {  /* |x| < 107 */
 354                 x = fabsl (x);
 355                 s = one / (x * x);
 356                 if (ix < 0x4000b6db) {  /* 1.25 <= |x| < 2.85711669921875 ~ 1/.35 */
 357                         R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
 358                              s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
 359                         S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
 360                              s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
 361                 } else if (ix < 0x4001d555) {  /* 6.6666259765625 > |x| >= 1/.35 ~ 2.857143 */
 362                         R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
 363                              s * (rb[5] + s * (rb[6] + s * rb[7]))))));
 364                         S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
 365                              s * (sb[5] + s * (sb[6] + s))))));
 366                 } else { /* 107 > |x| >= 6.666 */
 367                         if (se & 0x8000)
 368                                 return two - tiny;/* x < -6.666 */
 369                         R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
 370                              s * (rc[4] + s * rc[5]))));
 371                         S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
 372                              s * (sc[4] + s))));
 373                 }
 374                 z = x;
 375                 GET_LDOUBLE_WORDS (hx, i0, i1, z);
 376                 i1 = 0;
 377                 i0 &= 0xffffff00;
 378                 SET_LDOUBLE_WORDS (z, hx, i0, i1);
 379                 r = expl (-z * z - 0.5625) *
 380                 expl ((z - x) * (z + x) + R / S);
 381                 if ((se & 0x8000) == 0)
 382                         return r / x;
 383                 return two - r / x;
 384         }
 385
 386         if ((se & 0x8000) == 0)
 387                 return tiny * tiny;
 388         return two - tiny;
 389 }
 390 #endif