sysdeps/ieee754/ldbl-128ibm/s_erfl.c

   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 /* Modifications and expansions for 128-bit long double are
  13    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14    and are incorporated herein by permission of the author.  The author
  15    reserves the right to distribute this material elsewhere under different
  16    copying permissions.  These modifications are distributed here under
  17    the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, see
  31     <http://www.gnu.org/licenses/>.  */
  32
  33 /* double erf(double x)
  34  * double erfc(double x)
  35  *                           x
  36  *                    2      |\
  37  *     erf(x)  =  ---------  | exp(-t*t)dt
  38  *                 sqrt(pi) \|
  39  *                           0
  40  *
  41  *     erfc(x) =  1-erf(x)
  42  *  Note that
  43  *              erf(-x) = -erf(x)
  44  *              erfc(-x) = 2 - erfc(x)
  45  *
  46  * Method:
  47  *      1.  erf(x)  = x + x*R(x^2) for |x| in [0, 7/8]
  48  *         Remark. The formula is derived by noting
  49  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
  50  *         and that
  51  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
  52  *         is close to one.
  53  *
  54  *      1a. erf(x)  = 1 - erfc(x), for |x| > 1.0
  55  *          erfc(x) = 1 - erf(x)  if |x| < 1/4
  56  *
  57  *      2. For |x| in [7/8, 1], let s = |x| - 1, and
  58  *         c = 0.84506291151 rounded to single (24 bits)
  59  *      erf(s + c)  = sign(x) * (c  + P1(s)/Q1(s))
  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/4, 5/4],
  66  *      erfc(s + const)  = erfc(const)  + s P1(s)/Q1(s)
  67  *              for const = 1/4, 3/8, ..., 9/8
  68  *              and 0 <= s <= 1/8 .
  69  *
  70  *      4. For x in [5/4, 107],
  71  *      erfc(x) = (1/x)*exp(-x*x-0.5625 + R(z))
  72  *              z=1/x^2
  73  *         The interval is partitioned into several segments
  74  *         of width 1/8 in 1/x.
  75  *      erf(x) = 1.0 - erfc(x) if x < 25.6283 else
  76  *      erf(x) = sign(x)*(1.0 - tiny)
  77  *
  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  *      Note3:
  91  *         For x higher than 25.6283, erf(x) underflows.
  92  *
  93  *      5. For inf > x >= 107
  94  *      erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
  95  *      erfc(x) = tiny*tiny (raise underflow) if x > 0
  96  *                      = 2 - tiny if x<0
  97  *
  98  *      7. Special case:
  99  *      erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
 100  *      erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
 101  *              erfc/erf(NaN) is NaN
 102  */
 103
 104 #include <errno.h>
 105 #include <float.h>
 106 #include <math.h>
 107 #include <math_private.h>
 108 #include <math_ldbl_opt.h>
 109
 110 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 111
 112 static long double
 113 neval (long double x, const long double *p, int n)
 114 {
 115   long double y;
 116
 117   p += n;
 118   y = *p--;
 119   do
 120     {
 121       y = y * x + *p--;
 122     }
 123   while (--n > 0);
 124   return y;
 125 }
 126
 127
 128 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 129
 130 static long double
 131 deval (long double x, const long double *p, int n)
 132 {
 133   long double y;
 134
 135   p += n;
 136   y = x + *p--;
 137   do
 138     {
 139       y = y * x + *p--;
 140     }
 141   while (--n > 0);
 142   return y;
 143 }
 144
 145
 146
 147 static const long double
 148 tiny = 1e-300L,
 149   half = 0.5L,
 150   one = 1.0L,
 151   two = 2.0L,
 152   /* 2/sqrt(pi) - 1 */
 153   efx = 1.2837916709551257389615890312154517168810E-1L;
 154
 155
 156 /* erf(x)  = x  + x R(x^2)
 157    0 <= x <= 7/8
 158    Peak relative error 1.8e-35  */
 159 #define NTN1 8
 160 static const long double TN1[NTN1 + 1] =
 161 {
 162  -3.858252324254637124543172907442106422373E10L,
 163   9.580319248590464682316366876952214879858E10L,
 164   1.302170519734879977595901236693040544854E10L,
 165   2.922956950426397417800321486727032845006E9L,
 166   1.764317520783319397868923218385468729799E8L,
 167   1.573436014601118630105796794840834145120E7L,
 168   4.028077380105721388745632295157816229289E5L,
 169   1.644056806467289066852135096352853491530E4L,
 170   3.390868480059991640235675479463287886081E1L
 171 };
 172 #define NTD1 8
 173 static const long double TD1[NTD1 + 1] =
 174 {
 175   -3.005357030696532927149885530689529032152E11L,
 176   -1.342602283126282827411658673839982164042E11L,
 177   -2.777153893355340961288511024443668743399E10L,
 178   -3.483826391033531996955620074072768276974E9L,
 179   -2.906321047071299585682722511260895227921E8L,
 180   -1.653347985722154162439387878512427542691E7L,
 181   -6.245520581562848778466500301865173123136E5L,
 182   -1.402124304177498828590239373389110545142E4L,
 183   -1.209368072473510674493129989468348633579E2L
 184 /* 1.0E0 */
 185 };
 186
 187
 188 /* erf(z+1)  = erf_const + P(z)/Q(z)
 189    -.125 <= z <= 0
 190    Peak relative error 7.3e-36  */
 191 static const long double erf_const = 0.845062911510467529296875L;
 192 #define NTN2 8
 193 static const long double TN2[NTN2 + 1] =
 194 {
 195  -4.088889697077485301010486931817357000235E1L,
 196   7.157046430681808553842307502826960051036E3L,
 197  -2.191561912574409865550015485451373731780E3L,
 198   2.180174916555316874988981177654057337219E3L,
 199   2.848578658049670668231333682379720943455E2L,
 200   1.630362490952512836762810462174798925274E2L,
 201   6.317712353961866974143739396865293596895E0L,
 202   2.450441034183492434655586496522857578066E1L,
 203   5.127662277706787664956025545897050896203E-1L
 204 };
 205 #define NTD2 8
 206 static const long double TD2[NTD2 + 1] =
 207 {
 208   1.731026445926834008273768924015161048885E4L,
 209   1.209682239007990370796112604286048173750E4L,
 210   1.160950290217993641320602282462976163857E4L,
 211   5.394294645127126577825507169061355698157E3L,
 212   2.791239340533632669442158497532521776093E3L,
 213   8.989365571337319032943005387378993827684E2L,
 214   2.974016493766349409725385710897298069677E2L,
 215   6.148192754590376378740261072533527271947E1L,
 216   1.178502892490738445655468927408440847480E1L
 217  /* 1.0E0 */
 218 };
 219
 220
 221 /* erfc(x + 0.25) = erfc(0.25) + x R(x)
 222    0 <= x < 0.125
 223    Peak relative error 1.4e-35  */
 224 #define NRNr13 8
 225 static const long double RNr13[NRNr13 + 1] =
 226 {
 227  -2.353707097641280550282633036456457014829E3L,
 228   3.871159656228743599994116143079870279866E2L,
 229  -3.888105134258266192210485617504098426679E2L,
 230  -2.129998539120061668038806696199343094971E1L,
 231  -8.125462263594034672468446317145384108734E1L,
 232   8.151549093983505810118308635926270319660E0L,
 233  -5.033362032729207310462422357772568553670E0L,
 234  -4.253956621135136090295893547735851168471E-2L,
 235  -8.098602878463854789780108161581050357814E-2L
 236 };
 237 #define NRDr13 7
 238 static const long double RDr13[NRDr13 + 1] =
 239 {
 240   2.220448796306693503549505450626652881752E3L,
 241   1.899133258779578688791041599040951431383E2L,
 242   1.061906712284961110196427571557149268454E3L,
 243   7.497086072306967965180978101974566760042E1L,
 244   2.146796115662672795876463568170441327274E2L,
 245   1.120156008362573736664338015952284925592E1L,
 246   2.211014952075052616409845051695042741074E1L,
 247   6.469655675326150785692908453094054988938E-1L
 248  /* 1.0E0 */
 249 };
 250 /* erfc(0.25) = C13a + C13b to extra precision.  */
 251 static const long double C13a = 0.723663330078125L;
 252 static const long double C13b = 1.0279753638067014931732235184287934646022E-5L;
 253
 254
 255 /* erfc(x + 0.375) = erfc(0.375) + x R(x)
 256    0 <= x < 0.125
 257    Peak relative error 1.2e-35  */
 258 #define NRNr14 8
 259 static const long double RNr14[NRNr14 + 1] =
 260 {
 261  -2.446164016404426277577283038988918202456E3L,
 262   6.718753324496563913392217011618096698140E2L,
 263  -4.581631138049836157425391886957389240794E2L,
 264  -2.382844088987092233033215402335026078208E1L,
 265  -7.119237852400600507927038680970936336458E1L,
 266   1.313609646108420136332418282286454287146E1L,
 267  -6.188608702082264389155862490056401365834E0L,
 268  -2.787116601106678287277373011101132659279E-2L,
 269  -2.230395570574153963203348263549700967918E-2L
 270 };
 271 #define NRDr14 7
 272 static const long double RDr14[NRDr14 + 1] =
 273 {
 274   2.495187439241869732696223349840963702875E3L,
 275   2.503549449872925580011284635695738412162E2L,
 276   1.159033560988895481698051531263861842461E3L,
 277   9.493751466542304491261487998684383688622E1L,
 278   2.276214929562354328261422263078480321204E2L,
 279   1.367697521219069280358984081407807931847E1L,
 280   2.276988395995528495055594829206582732682E1L,
 281   7.647745753648996559837591812375456641163E-1L
 282  /* 1.0E0 */
 283 };
 284 /* erfc(0.375) = C14a + C14b to extra precision.  */
 285 static const long double C14a = 0.5958709716796875L;
 286 static const long double C14b = 1.2118885490201676174914080878232469565953E-5L;
 287
 288 /* erfc(x + 0.5) = erfc(0.5) + x R(x)
 289    0 <= x < 0.125
 290    Peak relative error 4.7e-36  */
 291 #define NRNr15 8
 292 static const long double RNr15[NRNr15 + 1] =
 293 {
 294  -2.624212418011181487924855581955853461925E3L,
 295   8.473828904647825181073831556439301342756E2L,
 296  -5.286207458628380765099405359607331669027E2L,
 297  -3.895781234155315729088407259045269652318E1L,
 298  -6.200857908065163618041240848728398496256E1L,
 299   1.469324610346924001393137895116129204737E1L,
 300  -6.961356525370658572800674953305625578903E0L,
 301   5.145724386641163809595512876629030548495E-3L,
 302   1.990253655948179713415957791776180406812E-2L
 303 };
 304 #define NRDr15 7
 305 static const long double RDr15[NRDr15 + 1] =
 306 {
 307   2.986190760847974943034021764693341524962E3L,
 308   5.288262758961073066335410218650047725985E2L,
 309   1.363649178071006978355113026427856008978E3L,
 310   1.921707975649915894241864988942255320833E2L,
 311   2.588651100651029023069013885900085533226E2L,
 312   2.628752920321455606558942309396855629459E1L,
 313   2.455649035885114308978333741080991380610E1L,
 314   1.378826653595128464383127836412100939126E0L
 315   /* 1.0E0 */
 316 };
 317 /* erfc(0.5) = C15a + C15b to extra precision.  */
 318 static const long double C15a = 0.4794921875L;
 319 static const long double C15b = 7.9346869534623172533461080354712635484242E-6L;
 320
 321 /* erfc(x + 0.625) = erfc(0.625) + x R(x)
 322    0 <= x < 0.125
 323    Peak relative error 5.1e-36  */
 324 #define NRNr16 8
 325 static const long double RNr16[NRNr16 + 1] =
 326 {
 327  -2.347887943200680563784690094002722906820E3L,
 328   8.008590660692105004780722726421020136482E2L,
 329  -5.257363310384119728760181252132311447963E2L,
 330  -4.471737717857801230450290232600243795637E1L,
 331  -4.849540386452573306708795324759300320304E1L,
 332   1.140885264677134679275986782978655952843E1L,
 333  -6.731591085460269447926746876983786152300E0L,
 334   1.370831653033047440345050025876085121231E-1L,
 335   2.022958279982138755020825717073966576670E-2L,
 336 };
 337 #define NRDr16 7
 338 static const long double RDr16[NRDr16 + 1] =
 339 {
 340   3.075166170024837215399323264868308087281E3L,
 341   8.730468942160798031608053127270430036627E2L,
 342   1.458472799166340479742581949088453244767E3L,
 343   3.230423687568019709453130785873540386217E2L,
 344   2.804009872719893612081109617983169474655E2L,
 345   4.465334221323222943418085830026979293091E1L,
 346   2.612723259683205928103787842214809134746E1L,
 347   2.341526751185244109722204018543276124997E0L,
 348   /* 1.0E0 */
 349 };
 350 /* erfc(0.625) = C16a + C16b to extra precision.  */
 351 static const long double C16a = 0.3767547607421875L;
 352 static const long double C16b = 4.3570693945275513594941232097252997287766E-6L;
 353
 354 /* erfc(x + 0.75) = erfc(0.75) + x R(x)
 355    0 <= x < 0.125
 356    Peak relative error 1.7e-35  */
 357 #define NRNr17 8
 358 static const long double RNr17[NRNr17 + 1] =
 359 {
 360   -1.767068734220277728233364375724380366826E3L,
 361   6.693746645665242832426891888805363898707E2L,
 362   -4.746224241837275958126060307406616817753E2L,
 363   -2.274160637728782675145666064841883803196E1L,
 364   -3.541232266140939050094370552538987982637E1L,
 365   6.988950514747052676394491563585179503865E0L,
 366   -5.807687216836540830881352383529281215100E0L,
 367   3.631915988567346438830283503729569443642E-1L,
 368   -1.488945487149634820537348176770282391202E-2L
 369 };
 370 #define NRDr17 7
 371 static const long double RDr17[NRDr17 + 1] =
 372 {
 373   2.748457523498150741964464942246913394647E3L,
 374   1.020213390713477686776037331757871252652E3L,
 375   1.388857635935432621972601695296561952738E3L,
 376   3.903363681143817750895999579637315491087E2L,
 377   2.784568344378139499217928969529219886578E2L,
 378   5.555800830216764702779238020065345401144E1L,
 379   2.646215470959050279430447295801291168941E1L,
 380   2.984905282103517497081766758550112011265E0L,
 381   /* 1.0E0 */
 382 };
 383 /* erfc(0.75) = C17a + C17b to extra precision.  */
 384 static const long double C17a = 0.2888336181640625L;
 385 static const long double C17b = 1.0748182422368401062165408589222625794046E-5L;
 386
 387
 388 /* erfc(x + 0.875) = erfc(0.875) + x R(x)
 389    0 <= x < 0.125
 390    Peak relative error 2.2e-35  */
 391 #define NRNr18 8
 392 static const long double RNr18[NRNr18 + 1] =
 393 {
 394  -1.342044899087593397419622771847219619588E3L,
 395   6.127221294229172997509252330961641850598E2L,
 396  -4.519821356522291185621206350470820610727E2L,
 397   1.223275177825128732497510264197915160235E1L,
 398  -2.730789571382971355625020710543532867692E1L,
 399   4.045181204921538886880171727755445395862E0L,
 400  -4.925146477876592723401384464691452700539E0L,
 401   5.933878036611279244654299924101068088582E-1L,
 402  -5.557645435858916025452563379795159124753E-2L
 403 };
 404 #define NRDr18 7
 405 static const long double RDr18[NRDr18 + 1] =
 406 {
 407   2.557518000661700588758505116291983092951E3L,
 408   1.070171433382888994954602511991940418588E3L,
 409   1.344842834423493081054489613250688918709E3L,
 410   4.161144478449381901208660598266288188426E2L,
 411   2.763670252219855198052378138756906980422E2L,
 412   5.998153487868943708236273854747564557632E1L,
 413   2.657695108438628847733050476209037025318E1L,
 414   3.252140524394421868923289114410336976512E0L,
 415   /* 1.0E0 */
 416 };
 417 /* erfc(0.875) = C18a + C18b to extra precision.  */
 418 static const long double C18a = 0.215911865234375L;
 419 static const long double C18b = 1.3073705765341685464282101150637224028267E-5L;
 420
 421 /* erfc(x + 1.0) = erfc(1.0) + x R(x)
 422    0 <= x < 0.125
 423    Peak relative error 1.6e-35  */
 424 #define NRNr19 8
 425 static const long double RNr19[NRNr19 + 1] =
 426 {
 427  -1.139180936454157193495882956565663294826E3L,
 428   6.134903129086899737514712477207945973616E2L,
 429  -4.628909024715329562325555164720732868263E2L,
 430   4.165702387210732352564932347500364010833E1L,
 431  -2.286979913515229747204101330405771801610E1L,
 432   1.870695256449872743066783202326943667722E0L,
 433  -4.177486601273105752879868187237000032364E0L,
 434   7.533980372789646140112424811291782526263E-1L,
 435  -8.629945436917752003058064731308767664446E-2L
 436 };
 437 #define NRDr19 7
 438 static const long double RDr19[NRDr19 + 1] =
 439 {
 440   2.744303447981132701432716278363418643778E3L,
 441   1.266396359526187065222528050591302171471E3L,
 442   1.466739461422073351497972255511919814273E3L,
 443   4.868710570759693955597496520298058147162E2L,
 444   2.993694301559756046478189634131722579643E2L,
 445   6.868976819510254139741559102693828237440E1L,
 446   2.801505816247677193480190483913753613630E1L,
 447   3.604439909194350263552750347742663954481E0L,
 448   /* 1.0E0 */
 449 };
 450 /* erfc(1.0) = C19a + C19b to extra precision.  */
 451 static const long double C19a = 0.15728759765625L;
 452 static const long double C19b = 1.1609394035130658779364917390740703933002E-5L;
 453
 454 /* erfc(x + 1.125) = erfc(1.125) + x R(x)
 455    0 <= x < 0.125
 456    Peak relative error 3.6e-36  */
 457 #define NRNr20 8
 458 static const long double RNr20[NRNr20 + 1] =
 459 {
 460  -9.652706916457973956366721379612508047640E2L,
 461   5.577066396050932776683469951773643880634E2L,
 462  -4.406335508848496713572223098693575485978E2L,
 463   5.202893466490242733570232680736966655434E1L,
 464  -1.931311847665757913322495948705563937159E1L,
 465  -9.364318268748287664267341457164918090611E-2L,
 466  -3.306390351286352764891355375882586201069E0L,
 467   7.573806045289044647727613003096916516475E-1L,
 468  -9.611744011489092894027478899545635991213E-2L
 469 };
 470 #define NRDr20 7
 471 static const long double RDr20[NRDr20 + 1] =
 472 {
 473   3.032829629520142564106649167182428189014E3L,
 474   1.659648470721967719961167083684972196891E3L,
 475   1.703545128657284619402511356932569292535E3L,
 476   6.393465677731598872500200253155257708763E2L,
 477   3.489131397281030947405287112726059221934E2L,
 478   8.848641738570783406484348434387611713070E1L,
 479   3.132269062552392974833215844236160958502E1L,
 480   4.430131663290563523933419966185230513168E0L
 481  /* 1.0E0 */
 482 };
 483 /* erfc(1.125) = C20a + C20b to extra precision.  */
 484 static const long double C20a = 0.111602783203125L;
 485 static const long double C20b = 8.9850951672359304215530728365232161564636E-6L;
 486
 487 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
 488    7/8 <= 1/x < 1
 489    Peak relative error 1.4e-35  */
 490 #define NRNr8 9
 491 static const long double RNr8[NRNr8 + 1] =
 492 {
 493   3.587451489255356250759834295199296936784E1L,
 494   5.406249749087340431871378009874875889602E2L,
 495   2.931301290625250886238822286506381194157E3L,
 496   7.359254185241795584113047248898753470923E3L,
 497   9.201031849810636104112101947312492532314E3L,
 498   5.749697096193191467751650366613289284777E3L,
 499   1.710415234419860825710780802678697889231E3L,
 500   2.150753982543378580859546706243022719599E2L,
 501   8.740953582272147335100537849981160931197E0L,
 502   4.876422978828717219629814794707963640913E-2L
 503 };
 504 #define NRDr8 8
 505 static const long double RDr8[NRDr8 + 1] =
 506 {
 507   6.358593134096908350929496535931630140282E1L,
 508   9.900253816552450073757174323424051765523E2L,
 509   5.642928777856801020545245437089490805186E3L,
 510   1.524195375199570868195152698617273739609E4L,
 511   2.113829644500006749947332935305800887345E4L,
 512   1.526438562626465706267943737310282977138E4L,
 513   5.561370922149241457131421914140039411782E3L,
 514   9.394035530179705051609070428036834496942E2L,
 515   6.147019596150394577984175188032707343615E1L
 516   /* 1.0E0 */
 517 };
 518
 519 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
 520    0.75 <= 1/x <= 0.875
 521    Peak relative error 2.0e-36  */
 522 #define NRNr7 9
 523 static const long double RNr7[NRNr7 + 1] =
 524 {
 525  1.686222193385987690785945787708644476545E1L,
 526  1.178224543567604215602418571310612066594E3L,
 527  1.764550584290149466653899886088166091093E4L,
 528  1.073758321890334822002849369898232811561E5L,
 529  3.132840749205943137619839114451290324371E5L,
 530  4.607864939974100224615527007793867585915E5L,
 531  3.389781820105852303125270837910972384510E5L,
 532  1.174042187110565202875011358512564753399E5L,
 533  1.660013606011167144046604892622504338313E4L,
 534  6.700393957480661937695573729183733234400E2L
 535 };
 536 #define NRDr7 9
 537 static const long double RDr7[NRDr7 + 1] =
 538 {
 539 -1.709305024718358874701575813642933561169E3L,
 540 -3.280033887481333199580464617020514788369E4L,
 541 -2.345284228022521885093072363418750835214E5L,
 542 -8.086758123097763971926711729242327554917E5L,
 543 -1.456900414510108718402423999575992450138E6L,
 544 -1.391654264881255068392389037292702041855E6L,
 545 -6.842360801869939983674527468509852583855E5L,
 546 -1.597430214446573566179675395199807533371E5L,
 547 -1.488876130609876681421645314851760773480E4L,
 548 -3.511762950935060301403599443436465645703E2L
 549  /* 1.0E0 */
 550 };
 551
 552 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
 553    5/8 <= 1/x < 3/4
 554    Peak relative error 1.9e-35  */
 555 #define NRNr6 9
 556 static const long double RNr6[NRNr6 + 1] =
 557 {
 558  1.642076876176834390623842732352935761108E0L,
 559  1.207150003611117689000664385596211076662E2L,
 560  2.119260779316389904742873816462800103939E3L,
 561  1.562942227734663441801452930916044224174E4L,
 562  5.656779189549710079988084081145693580479E4L,
 563  1.052166241021481691922831746350942786299E5L,
 564  9.949798524786000595621602790068349165758E4L,
 565  4.491790734080265043407035220188849562856E4L,
 566  8.377074098301530326270432059434791287601E3L,
 567  4.506934806567986810091824791963991057083E2L
 568 };
 569 #define NRDr6 9
 570 static const long double RDr6[NRDr6 + 1] =
 571 {
 572 -1.664557643928263091879301304019826629067E2L,
 573 -3.800035902507656624590531122291160668452E3L,
 574 -3.277028191591734928360050685359277076056E4L,
 575 -1.381359471502885446400589109566587443987E5L,
 576 -3.082204287382581873532528989283748656546E5L,
 577 -3.691071488256738343008271448234631037095E5L,
 578 -2.300482443038349815750714219117566715043E5L,
 579 -6.873955300927636236692803579555752171530E4L,
 580 -8.262158817978334142081581542749986845399E3L,
 581 -2.517122254384430859629423488157361983661E2L
 582  /* 1.00 */
 583 };
 584
 585 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
 586    1/2 <= 1/x < 5/8
 587    Peak relative error 4.6e-36  */
 588 #define NRNr5 10
 589 static const long double RNr5[NRNr5 + 1] =
 590 {
 591 -3.332258927455285458355550878136506961608E-3L,
 592 -2.697100758900280402659586595884478660721E-1L,
 593 -6.083328551139621521416618424949137195536E0L,
 594 -6.119863528983308012970821226810162441263E1L,
 595 -3.176535282475593173248810678636522589861E2L,
 596 -8.933395175080560925809992467187963260693E2L,
 597 -1.360019508488475978060917477620199499560E3L,
 598 -1.075075579828188621541398761300910213280E3L,
 599 -4.017346561586014822824459436695197089916E2L,
 600 -5.857581368145266249509589726077645791341E1L,
 601 -2.077715925587834606379119585995758954399E0L
 602 };
 603 #define NRDr5 9
 604 static const long double RDr5[NRDr5 + 1] =
 605 {
 606  3.377879570417399341550710467744693125385E-1L,
 607  1.021963322742390735430008860602594456187E1L,
 608  1.200847646592942095192766255154827011939E2L,
 609  7.118915528142927104078182863387116942836E2L,
 610  2.318159380062066469386544552429625026238E3L,
 611  4.238729853534009221025582008928765281620E3L,
 612  4.279114907284825886266493994833515580782E3L,
 613  2.257277186663261531053293222591851737504E3L,
 614  5.570475501285054293371908382916063822957E2L,
 615  5.142189243856288981145786492585432443560E1L
 616  /* 1.0E0 */
 617 };
 618
 619 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
 620    3/8 <= 1/x < 1/2
 621    Peak relative error 2.0e-36  */
 622 #define NRNr4 10
 623 static const long double RNr4[NRNr4 + 1] =
 624 {
 625  3.258530712024527835089319075288494524465E-3L,
 626  2.987056016877277929720231688689431056567E-1L,
 627  8.738729089340199750734409156830371528862E0L,
 628  1.207211160148647782396337792426311125923E2L,
 629  8.997558632489032902250523945248208224445E2L,
 630  3.798025197699757225978410230530640879762E3L,
 631  9.113203668683080975637043118209210146846E3L,
 632  1.203285891339933238608683715194034900149E4L,
 633  8.100647057919140328536743641735339740855E3L,
 634  2.383888249907144945837976899822927411769E3L,
 635  2.127493573166454249221983582495245662319E2L
 636 };
 637 #define NRDr4 10
 638 static const long double RDr4[NRDr4 + 1] =
 639 {
 640 -3.303141981514540274165450687270180479586E-1L,
 641 -1.353768629363605300707949368917687066724E1L,
 642 -2.206127630303621521950193783894598987033E2L,
 643 -1.861800338758066696514480386180875607204E3L,
 644 -8.889048775872605708249140016201753255599E3L,
 645 -2.465888106627948210478692168261494857089E4L,
 646 -3.934642211710774494879042116768390014289E4L,
 647 -3.455077258242252974937480623730228841003E4L,
 648 -1.524083977439690284820586063729912653196E4L,
 649 -2.810541887397984804237552337349093953857E3L,
 650 -1.343929553541159933824901621702567066156E2L
 651  /* 1.0E0 */
 652 };
 653
 654 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
 655    1/4 <= 1/x < 3/8
 656    Peak relative error 8.4e-37  */
 657 #define NRNr3 11
 658 static const long double RNr3[NRNr3 + 1] =
 659 {
 660 -1.952401126551202208698629992497306292987E-6L,
 661 -2.130881743066372952515162564941682716125E-4L,
 662 -8.376493958090190943737529486107282224387E-3L,
 663 -1.650592646560987700661598877522831234791E-1L,
 664 -1.839290818933317338111364667708678163199E0L,
 665 -1.216278715570882422410442318517814388470E1L,
 666 -4.818759344462360427612133632533779091386E1L,
 667 -1.120994661297476876804405329172164436784E2L,
 668 -1.452850765662319264191141091859300126931E2L,
 669 -9.485207851128957108648038238656777241333E1L,
 670 -2.563663855025796641216191848818620020073E1L,
 671 -1.787995944187565676837847610706317833247E0L
 672 };
 673 #define NRDr3 10
 674 static const long double RDr3[NRDr3 + 1] =
 675 {
 676  1.979130686770349481460559711878399476903E-4L,
 677  1.156941716128488266238105813374635099057E-2L,
 678  2.752657634309886336431266395637285974292E-1L,
 679  3.482245457248318787349778336603569327521E0L,
 680  2.569347069372696358578399521203959253162E1L,
 681  1.142279000180457419740314694631879921561E2L,
 682  3.056503977190564294341422623108332700840E2L,
 683  4.780844020923794821656358157128719184422E2L,
 684  4.105972727212554277496256802312730410518E2L,
 685  1.724072188063746970865027817017067646246E2L,
 686  2.815939183464818198705278118326590370435E1L
 687  /* 1.0E0 */
 688 };
 689
 690 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
 691    1/8 <= 1/x < 1/4
 692    Peak relative error 1.5e-36  */
 693 #define NRNr2 11
 694 static const long double RNr2[NRNr2 + 1] =
 695 {
 696 -2.638914383420287212401687401284326363787E-8L,
 697 -3.479198370260633977258201271399116766619E-6L,
 698 -1.783985295335697686382487087502222519983E-4L,
 699 -4.777876933122576014266349277217559356276E-3L,
 700 -7.450634738987325004070761301045014986520E-2L,
 701 -7.068318854874733315971973707247467326619E-1L,
 702 -4.113919921935944795764071670806867038732E0L,
 703 -1.440447573226906222417767283691888875082E1L,
 704 -2.883484031530718428417168042141288943905E1L,
 705 -2.990886974328476387277797361464279931446E1L,
 706 -1.325283914915104866248279787536128997331E1L,
 707 -1.572436106228070195510230310658206154374E0L
 708 };
 709 #define NRDr2 10
 710 static const long double RDr2[NRDr2 + 1] =
 711 {
 712  2.675042728136731923554119302571867799673E-6L,
 713  2.170997868451812708585443282998329996268E-4L,
 714  7.249969752687540289422684951196241427445E-3L,
 715  1.302040375859768674620410563307838448508E-1L,
 716  1.380202483082910888897654537144485285549E0L,
 717  8.926594113174165352623847870299170069350E0L,
 718  3.521089584782616472372909095331572607185E1L,
 719  8.233547427533181375185259050330809105570E1L,
 720  1.072971579885803033079469639073292840135E2L,
 721  6.943803113337964469736022094105143158033E1L,
 722  1.775695341031607738233608307835017282662E1L
 723  /* 1.0E0 */
 724 };
 725
 726 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
 727    1/128 <= 1/x < 1/8
 728    Peak relative error 2.2e-36  */
 729 #define NRNr1 9
 730 static const long double RNr1[NRNr1 + 1] =
 731 {
 732 -4.250780883202361946697751475473042685782E-8L,
 733 -5.375777053288612282487696975623206383019E-6L,
 734 -2.573645949220896816208565944117382460452E-4L,
 735 -6.199032928113542080263152610799113086319E-3L,
 736 -8.262721198693404060380104048479916247786E-2L,
 737 -6.242615227257324746371284637695778043982E-1L,
 738 -2.609874739199595400225113299437099626386E0L,
 739 -5.581967563336676737146358534602770006970E0L,
 740 -5.124398923356022609707490956634280573882E0L,
 741 -1.290865243944292370661544030414667556649E0L
 742 };
 743 #define NRDr1 8
 744 static const long double RDr1[NRDr1 + 1] =
 745 {
 746  4.308976661749509034845251315983612976224E-6L,
 747  3.265390126432780184125233455960049294580E-4L,
 748  9.811328839187040701901866531796570418691E-3L,
 749  1.511222515036021033410078631914783519649E-1L,
 750  1.289264341917429958858379585970225092274E0L,
 751  6.147640356182230769548007536914983522270E0L,
 752  1.573966871337739784518246317003956180750E1L,
 753  1.955534123435095067199574045529218238263E1L,
 754  9.472613121363135472247929109615785855865E0L
 755   /* 1.0E0 */
 756 };
 757
 758
 759 long double
 760 __erfl (long double x)
 761 {
 762   long double a, y, z;
 763   int32_t i, ix, hx;
 764   double xhi;
 765
 766   xhi = ldbl_high (x);
 767   GET_HIGH_WORD (hx, xhi);
 768   ix = hx & 0x7fffffff;
 769
 770   if (ix >= 0x7ff00000)
 771     {                           /* erf(nan)=nan */
 772       i = ((uint32_t) hx >> 31) << 1;
 773       return (long double) (1 - i) + one / x;   /* erf(+-inf)=+-1 */
 774     }
 775
 776   if (ix >= 0x3ff00000) /* |x| >= 1.0 */
 777     {
 778       if (ix >= 0x4039A0DE)
 779         {
 780         /* __erfcl (x) underflows if x > 25.6283 */
 781           if ((hx & 0x80000000) == 0)
 782             return one-tiny;
 783           else
 784             return tiny-one;
 785         }
 786       else
 787         {
 788           y = __erfcl (x);
 789           return (one - y);
 790         }
 791     }
 792   a = x;
 793   if ((hx & 0x80000000) != 0)
 794     a = -a;
 795   z = x * x;
 796   if (ix < 0x3fec0000)  /* a < 0.875 */
 797     {
 798       if (ix < 0x3c600000) /* |x|<2**-57 */
 799         {
 800           if (ix < 0x00800000)
 801             {
 802               /* erf (-0) = -0.  Unfortunately, for IBM extended double
 803                  0.0625 * (16.0 * x + (16.0 * efx) * x) for x = -0
 804                  evaluates to 0.  */
 805               if (x == 0)
 806                 return x;
 807               long double ret = 0.0625 * (16.0 * x + (16.0 * efx) * x);
 808               if (fabsl (ret) < LDBL_MIN)
 809                 {
 810                   long double force_underflow = ret * ret;
 811                   math_force_eval (force_underflow);
 812                 }
 813               return ret;
 814             }
 815           return x + efx * x;
 816         }
 817       y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
 818     }
 819   else
 820     {
 821       a = a - one;
 822       y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2);
 823     }
 824
 825   if (hx & 0x80000000) /* x < 0 */
 826     y = -y;
 827   return( y );
 828 }
 829
 830 long_double_symbol (libm, __erfl, erfl);
 831 long double
 832 __erfcl (long double x)
 833 {
 834   long double y, z, p, r;
 835   int32_t i, ix;
 836   uint32_t hx;
 837   double xhi;
 838
 839   xhi = ldbl_high (x);
 840   GET_HIGH_WORD (hx, xhi);
 841   ix = hx & 0x7fffffff;
 842
 843   if (ix >= 0x7ff00000)
 844     {                           /* erfc(nan)=nan */
 845       /* erfc(+-inf)=0,2 */
 846       return (long double) ((hx >> 31) << 1) + one / x;
 847     }
 848
 849   if (ix < 0x3fd00000) /* |x| <1/4 */
 850     {
 851       if (ix < 0x38d00000) /* |x|<2**-114 */
 852         return one - x;
 853       return one - __erfl (x);
 854     }
 855   if (ix < 0x3ff40000) /* 1.25 */
 856     {
 857       if ((hx & 0x80000000) != 0)
 858         x = -x;
 859       i = 8.0 * x;
 860       switch (i)
 861         {
 862         case 2:
 863           z = x - 0.25L;
 864           y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13);
 865           y += C13a;
 866           break;
 867         case 3:
 868           z = x - 0.375L;
 869           y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14);
 870           y += C14a;
 871           break;
 872         case 4:
 873           z = x - 0.5L;
 874           y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15);
 875           y += C15a;
 876           break;
 877         case 5:
 878           z = x - 0.625L;
 879           y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16);
 880           y += C16a;
 881           break;
 882         case 6:
 883           z = x - 0.75L;
 884           y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17);
 885           y += C17a;
 886           break;
 887         case 7:
 888           z = x - 0.875L;
 889           y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
 890           y += C18a;
 891           break;
 892         case 8:
 893           z = x - 1.0L;
 894           y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
 895           y += C19a;
 896           break;
 897         default: /* i == 9.  */
 898           z = x - 1.125L;
 899           y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
 900           y += C20a;
 901           break;
 902         }
 903       if (hx & 0x80000000)
 904         y = 2.0L - y;
 905       return y;
 906     }
 907   /* 1.25 < |x| < 107 */
 908   if (ix < 0x405ac000)
 909     {
 910       /* x < -9 */
 911       if (hx >= 0xc0220000)
 912         return two - tiny;
 913
 914       if ((hx & 0x80000000) != 0)
 915         x = -x;
 916       z = one / (x * x);
 917       i = 8.0 / x;
 918       switch (i)
 919         {
 920         default:
 921         case 0:
 922           p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1);
 923           break;
 924         case 1:
 925           p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2);
 926           break;
 927         case 2:
 928           p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3);
 929           break;
 930         case 3:
 931           p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4);
 932           break;
 933         case 4:
 934           p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5);
 935           break;
 936         case 5:
 937           p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6);
 938           break;
 939         case 6:
 940           p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7);
 941           break;
 942         case 7:
 943           p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8);
 944           break;
 945         }
 946       z = (float) x;
 947       r = __ieee754_expl (-z * z - 0.5625) *
 948         __ieee754_expl ((z - x) * (z + x) + p);
 949       if ((hx & 0x80000000) == 0)
 950         {
 951           long double ret = r / x;
 952           if (ret == 0)
 953             __set_errno (ERANGE);
 954           return ret;
 955         }
 956       else
 957         return two - r / x;
 958     }
 959   else
 960     {
 961       if ((hx & 0x80000000) == 0)
 962         {
 963           __set_errno (ERANGE);
 964           return tiny * tiny;
 965         }
 966       else
 967         return two - tiny;
 968     }
 969 }
 970
 971 long_double_symbol (libm, __erfcl, erfcl);