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