Back to index

glibc  2.9
s_erfl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /* Long double expansions are
00013   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00014   and are incorporated herein by permission of the author.  The author 
00015   reserves the right to distribute this material elsewhere under different
00016   copying permissions.  These modifications are distributed here under 
00017   the following terms:
00018 
00019     This library is free software; you can redistribute it and/or
00020     modify it under the terms of the GNU Lesser General Public
00021     License as published by the Free Software Foundation; either
00022     version 2.1 of the License, or (at your option) any later version.
00023 
00024     This library is distributed in the hope that it will be useful,
00025     but WITHOUT ANY WARRANTY; without even the implied warranty of
00026     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027     Lesser General Public License for more details.
00028 
00029     You should have received a copy of the GNU Lesser General Public
00030     License along with this library; if not, write to the Free Software
00031     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00032 
00033 /* double erf(double x)
00034  * double erfc(double x)
00035  *                        x
00036  *                  2      |\
00037  *     erf(x)  =  ---------  | exp(-t*t)dt
00038  *               sqrt(pi) \|
00039  *                        0
00040  *
00041  *     erfc(x) =  1-erf(x)
00042  *  Note that
00043  *            erf(-x) = -erf(x)
00044  *            erfc(-x) = 2 - erfc(x)
00045  *
00046  * Method:
00047  *     1. For |x| in [0, 0.84375]
00048  *         erf(x)  = x + x*R(x^2)
00049  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
00050  *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
00051  *        Remark. The formula is derived by noting
00052  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
00053  *        and that
00054  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
00055  *        is close to one. The interval is chosen because the fix
00056  *        point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
00057  *        near 0.6174), and by some experiment, 0.84375 is chosen to
00058  *        guarantee the error is less than one ulp for erf.
00059  *
00060  *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
00061  *         c = 0.84506291151 rounded to single (24 bits)
00062  *     erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
00063  *     erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
00064  *                     1+(c+P1(s)/Q1(s))    if x < 0
00065  *        Remark: here we use the taylor series expansion at x=1.
00066  *            erf(1+s) = erf(1) + s*Poly(s)
00067  *                    = 0.845.. + P1(s)/Q1(s)
00068  *        Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
00069  *
00070  *      3. For x in [1.25,1/0.35(~2.857143)],
00071  *     erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
00072  *              z=1/x^2
00073  *     erf(x)  = 1 - erfc(x)
00074  *
00075  *      4. For x in [1/0.35,107]
00076  *     erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
00077  *                   = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
00078  *                             if -6.666<x<0
00079  *                   = 2.0 - tiny         (if x <= -6.666)
00080  *              z=1/x^2
00081  *     erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
00082  *     erf(x)  = sign(x)*(1.0 - tiny)
00083  *      Note1:
00084  *        To compute exp(-x*x-0.5625+R/S), let s be a single
00085  *        precision number and s := x; then
00086  *            -x*x = -s*s + (s-x)*(s+x)
00087  *             exp(-x*x-0.5626+R/S) =
00088  *                   exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
00089  *      Note2:
00090  *        Here 4 and 5 make use of the asymptotic series
00091  *                     exp(-x*x)
00092  *            erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
00093  *                     x*sqrt(pi)
00094  *
00095  *      5. For inf > x >= 107
00096  *     erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
00097  *     erfc(x) = tiny*tiny (raise underflow) if x > 0
00098  *                   = 2 - tiny if x<0
00099  *
00100  *      7. Special case:
00101  *     erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
00102  *     erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
00103  *            erfc/erf(NaN) is NaN
00104  */
00105 
00106 
00107 #include "math.h"
00108 #include "math_private.h"
00109 
00110 #ifdef __STDC__
00111 static const long double
00112 #else
00113 static long double
00114 #endif
00115 tiny = 1e-4931L,
00116   half = 0.5L,
00117   one = 1.0L,
00118   two = 2.0L,
00119        /* c = (float)0.84506291151 */
00120   erx = 0.845062911510467529296875L,
00121 /*
00122  * Coefficients for approximation to  erf on [0,0.84375]
00123  */
00124   /* 2/sqrt(pi) - 1 */
00125   efx = 1.2837916709551257389615890312154517168810E-1L,
00126   /* 8 * (2/sqrt(pi) - 1) */
00127   efx8 = 1.0270333367641005911692712249723613735048E0L,
00128 
00129   pp[6] = {
00130     1.122751350964552113068262337278335028553E6L,
00131     -2.808533301997696164408397079650699163276E6L,
00132     -3.314325479115357458197119660818768924100E5L,
00133     -6.848684465326256109712135497895525446398E4L,
00134     -2.657817695110739185591505062971929859314E3L,
00135     -1.655310302737837556654146291646499062882E2L,
00136   },
00137 
00138   qq[6] = {
00139     8.745588372054466262548908189000448124232E6L,
00140     3.746038264792471129367533128637019611485E6L,
00141     7.066358783162407559861156173539693900031E5L,
00142     7.448928604824620999413120955705448117056E4L,
00143     4.511583986730994111992253980546131408924E3L,
00144     1.368902937933296323345610240009071254014E2L,
00145     /* 1.000000000000000000000000000000000000000E0 */
00146   },
00147 
00148 /*
00149  * Coefficients for approximation to  erf  in [0.84375,1.25]
00150  */
00151 /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
00152    -0.15625 <= x <= +.25
00153    Peak relative error 8.5e-22  */
00154 
00155   pa[8] = {
00156     -1.076952146179812072156734957705102256059E0L,
00157      1.884814957770385593365179835059971587220E2L,
00158     -5.339153975012804282890066622962070115606E1L,
00159      4.435910679869176625928504532109635632618E1L,
00160      1.683219516032328828278557309642929135179E1L,
00161     -2.360236618396952560064259585299045804293E0L,
00162      1.852230047861891953244413872297940938041E0L,
00163      9.394994446747752308256773044667843200719E-2L,
00164   },
00165 
00166   qa[7] =  {
00167     4.559263722294508998149925774781887811255E2L,
00168     3.289248982200800575749795055149780689738E2L,
00169     2.846070965875643009598627918383314457912E2L,
00170     1.398715859064535039433275722017479994465E2L,
00171     6.060190733759793706299079050985358190726E1L,
00172     2.078695677795422351040502569964299664233E1L,
00173     4.641271134150895940966798357442234498546E0L,
00174     /* 1.000000000000000000000000000000000000000E0 */
00175   },
00176 
00177 /*
00178  * Coefficients for approximation to  erfc in [1.25,1/0.35]
00179  */
00180 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
00181    1/2.85711669921875 < 1/x < 1/1.25
00182    Peak relative error 3.1e-21  */
00183 
00184     ra[] = {
00185       1.363566591833846324191000679620738857234E-1L,
00186       1.018203167219873573808450274314658434507E1L,
00187       1.862359362334248675526472871224778045594E2L,
00188       1.411622588180721285284945138667933330348E3L,
00189       5.088538459741511988784440103218342840478E3L,
00190       8.928251553922176506858267311750789273656E3L,
00191       7.264436000148052545243018622742770549982E3L,
00192       2.387492459664548651671894725748959751119E3L,
00193       2.220916652813908085449221282808458466556E2L,
00194     },
00195 
00196     sa[] = {
00197       -1.382234625202480685182526402169222331847E1L,
00198       -3.315638835627950255832519203687435946482E2L,
00199       -2.949124863912936259747237164260785326692E3L,
00200       -1.246622099070875940506391433635999693661E4L,
00201       -2.673079795851665428695842853070996219632E4L,
00202       -2.880269786660559337358397106518918220991E4L,
00203       -1.450600228493968044773354186390390823713E4L,
00204       -2.874539731125893533960680525192064277816E3L,
00205       -1.402241261419067750237395034116942296027E2L,
00206       /* 1.000000000000000000000000000000000000000E0 */
00207     },
00208 /*
00209  * Coefficients for approximation to  erfc in [1/.35,107]
00210  */
00211 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
00212    1/6.6666259765625 < 1/x < 1/2.85711669921875
00213    Peak relative error 4.2e-22  */
00214     rb[] = {
00215       -4.869587348270494309550558460786501252369E-5L,
00216       -4.030199390527997378549161722412466959403E-3L,
00217       -9.434425866377037610206443566288917589122E-2L,
00218       -9.319032754357658601200655161585539404155E-1L,
00219       -4.273788174307459947350256581445442062291E0L,
00220       -8.842289940696150508373541814064198259278E0L,
00221       -7.069215249419887403187988144752613025255E0L,
00222       -1.401228723639514787920274427443330704764E0L,
00223     },
00224 
00225     sb[] = {
00226       4.936254964107175160157544545879293019085E-3L,
00227       1.583457624037795744377163924895349412015E-1L,
00228       1.850647991850328356622940552450636420484E0L,
00229       9.927611557279019463768050710008450625415E0L,
00230       2.531667257649436709617165336779212114570E1L,
00231       2.869752886406743386458304052862814690045E1L,
00232       1.182059497870819562441683560749192539345E1L,
00233       /* 1.000000000000000000000000000000000000000E0 */
00234     },
00235 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
00236    1/107 <= 1/x <= 1/6.6666259765625
00237    Peak relative error 1.1e-21  */
00238     rc[] = {
00239       -8.299617545269701963973537248996670806850E-5L,
00240       -6.243845685115818513578933902532056244108E-3L,
00241       -1.141667210620380223113693474478394397230E-1L,
00242       -7.521343797212024245375240432734425789409E-1L,
00243       -1.765321928311155824664963633786967602934E0L,
00244       -1.029403473103215800456761180695263439188E0L,
00245     },
00246 
00247     sc[] = {
00248       8.413244363014929493035952542677768808601E-3L,
00249       2.065114333816877479753334599639158060979E-1L,
00250       1.639064941530797583766364412782135680148E0L,
00251       4.936788463787115555582319302981666347450E0L,
00252       5.005177727208955487404729933261347679090E0L,
00253       /* 1.000000000000000000000000000000000000000E0 */
00254     };
00255 
00256 #ifdef __STDC__
00257 long double
00258 __erfl (long double x)
00259 #else
00260 long double
00261 __erfl (x)
00262      long double x;
00263 #endif
00264 {
00265   long double R, S, P, Q, s, y, z, r;
00266   int32_t ix, i;
00267   u_int32_t se, i0, i1;
00268 
00269   GET_LDOUBLE_WORDS (se, i0, i1, x);
00270   ix = se & 0x7fff;
00271 
00272   if (ix >= 0x7fff)
00273     {                       /* erf(nan)=nan */
00274       i = ((se & 0xffff) >> 15) << 1;
00275       return (long double) (1 - i) + one / x;    /* erf(+-inf)=+-1 */
00276     }
00277 
00278   ix = (ix << 16) | (i0 >> 16);
00279   if (ix < 0x3ffed800) /* |x|<0.84375 */
00280     {
00281       if (ix < 0x3fde8000) /* |x|<2**-33 */
00282        {
00283          if (ix < 0x00080000)
00284            return 0.125 * (8.0 * x + efx8 * x);  /*avoid underflow */
00285          return x + efx * x;
00286        }
00287       z = x * x;
00288       r = pp[0] + z * (pp[1]
00289           + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
00290       s = qq[0] + z * (qq[1]
00291          + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
00292       y = r / s;
00293       return x + x * y;
00294     }
00295   if (ix < 0x3fffa000) /* 1.25 */
00296     {                       /* 0.84375 <= |x| < 1.25 */
00297       s = fabsl (x) - one;
00298       P = pa[0] + s * (pa[1] + s * (pa[2]
00299         + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
00300       Q = qa[0] + s * (qa[1] + s * (qa[2]
00301         + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
00302       if ((se & 0x8000) == 0)
00303        return erx + P / Q;
00304       else
00305        return -erx - P / Q;
00306     }
00307   if (ix >= 0x4001d555) /* 6.6666259765625 */
00308     {                       /* inf>|x|>=6.666 */
00309       if ((se & 0x8000) == 0)
00310        return one - tiny;
00311       else
00312        return tiny - one;
00313     }
00314   x = fabsl (x);
00315   s = one / (x * x);
00316   if (ix < 0x4000b6db) /* 2.85711669921875 */
00317     {
00318       R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
00319           s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
00320       S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
00321           s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
00322     }
00323   else
00324     {                       /* |x| >= 1/0.35 */
00325       R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
00326          s * (rb[5] + s * (rb[6] + s * rb[7]))))));
00327       S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
00328          s * (sb[5] + s * (sb[6] + s))))));
00329     }
00330   z = x;
00331   GET_LDOUBLE_WORDS (i, i0, i1, z);
00332   i1 = 0;
00333   SET_LDOUBLE_WORDS (z, i, i0, i1);
00334   r =
00335     __ieee754_expl (-z * z - 0.5625) * __ieee754_expl ((z - x) * (z + x) +
00336                                                R / S);
00337   if ((se & 0x8000) == 0)
00338     return one - r / x;
00339   else
00340     return r / x - one;
00341 }
00342 
00343 weak_alias (__erfl, erfl)
00344 #ifdef __STDC__
00345      long double
00346      __erfcl (long double x)
00347 #else
00348      long double
00349      __erfcl (x)
00350      long double x;
00351 #endif
00352 {
00353   int32_t hx, ix;
00354   long double R, S, P, Q, s, y, z, r;
00355   u_int32_t se, i0, i1;
00356 
00357   GET_LDOUBLE_WORDS (se, i0, i1, x);
00358   ix = se & 0x7fff;
00359   if (ix >= 0x7fff)
00360     {                       /* erfc(nan)=nan */
00361       /* erfc(+-inf)=0,2 */
00362       return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
00363     }
00364 
00365   ix = (ix << 16) | (i0 >> 16);
00366   if (ix < 0x3ffed800) /* |x|<0.84375 */
00367     {
00368       if (ix < 0x3fbe0000) /* |x|<2**-65 */
00369        return one - x;
00370       z = x * x;
00371       r = pp[0] + z * (pp[1]
00372           + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
00373       s = qq[0] + z * (qq[1]
00374          + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
00375       y = r / s;
00376       if (ix < 0x3ffd8000) /* x<1/4 */
00377        {
00378          return one - (x + x * y);
00379        }
00380       else
00381        {
00382          r = x * y;
00383          r += (x - half);
00384          return half - r;
00385        }
00386     }
00387   if (ix < 0x3fffa000) /* 1.25 */
00388     {                       /* 0.84375 <= |x| < 1.25 */
00389       s = fabsl (x) - one;
00390       P = pa[0] + s * (pa[1] + s * (pa[2]
00391         + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
00392       Q = qa[0] + s * (qa[1] + s * (qa[2]
00393         + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
00394       if ((se & 0x8000) == 0)
00395        {
00396          z = one - erx;
00397          return z - P / Q;
00398        }
00399       else
00400        {
00401          z = erx + P / Q;
00402          return one + z;
00403        }
00404     }
00405   if (ix < 0x4005d600) /* 107 */
00406     {                       /* |x|<107 */
00407       x = fabsl (x);
00408       s = one / (x * x);
00409       if (ix < 0x4000b6db) /* 2.85711669921875 */
00410        {                    /* |x| < 1/.35 ~ 2.857143 */
00411          R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
00412               s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
00413          S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
00414               s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
00415        }
00416       else if (ix < 0x4001d555) /* 6.6666259765625 */
00417        {                    /* 6.666 > |x| >= 1/.35 ~ 2.857143 */
00418          R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
00419              s * (rb[5] + s * (rb[6] + s * rb[7]))))));
00420          S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
00421               s * (sb[5] + s * (sb[6] + s))))));
00422        }
00423       else
00424        {                    /* |x| >= 6.666 */
00425          if (se & 0x8000)
00426            return two - tiny;      /* x < -6.666 */
00427 
00428          R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
00429                                               s * (rc[4] + s * rc[5]))));
00430          S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
00431                                               s * (sc[4] + s))));
00432        }
00433       z = x;
00434       GET_LDOUBLE_WORDS (hx, i0, i1, z);
00435       i1 = 0;
00436       i0 &= 0xffffff00;
00437       SET_LDOUBLE_WORDS (z, hx, i0, i1);
00438       r = __ieee754_expl (-z * z - 0.5625) *
00439        __ieee754_expl ((z - x) * (z + x) + R / S);
00440       if ((se & 0x8000) == 0)
00441        return r / x;
00442       else
00443        return two - r / x;
00444     }
00445   else
00446     {
00447       if ((se & 0x8000) == 0)
00448        return tiny * tiny;
00449       else
00450        return two - tiny;
00451     }
00452 }
00453 
00454 weak_alias (__erfcl, erfcl)