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 #include <fix-int-fp-convert-zero.h>
 110
 111 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 112
 113 static long double
 114 neval (long double x, const long double *p, int n)
 115 {
 116   long double y;
 117
 118   p += n;
 119   y = *p--;
 120   do
 121     {
 122       y = y * x + *p--;
 123     }
 124   while (--n > 0);
 125   return y;
 126 }
 127
 128
 129 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
 130
 131 static long double
 132 deval (long double x, const long double *p, int n)
 133 {
 134   long double y;
 135
 136   p += n;
 137   y = x + *p--;
 138   do
 139     {
 140       y = y * x + *p--;
 141     }
 142   while (--n > 0);
 143   return y;
 144 }
 145
 146
 147
 148 static const long double
 149 tiny = 1e-300L,
 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               math_check_force_underflow (ret);
 809               return ret;
 810             }
 811           return x + efx * x;
 812         }
 813       y = a + a * neval (z, TN1, NTN1) / deval (z, TD1, NTD1);
 814     }
 815   else
 816     {
 817       a = a - one;
 818       y = erf_const + neval (a, TN2, NTN2) / deval (a, TD2, NTD2);
 819     }
 820
 821   if (hx & 0x80000000) /* x < 0 */
 822     y = -y;
 823   return( y );
 824 }
 825
 826 long_double_symbol (libm, __erfl, erfl);
 827 long double
 828 __erfcl (long double x)
 829 {
 830   long double y, z, p, r;
 831   int32_t i, ix;
 832   uint32_t hx;
 833   double xhi;
 834
 835   xhi = ldbl_high (x);
 836   GET_HIGH_WORD (hx, xhi);
 837   ix = hx & 0x7fffffff;
 838
 839   if (ix >= 0x7ff00000)
 840     {                           /* erfc(nan)=nan */
 841       /* erfc(+-inf)=0,2 */
 842       long double ret = (long double) ((hx >> 31) << 1) + one / x;
 843       if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0L)
 844         return 0.0L;
 845       return ret;
 846     }
 847
 848   if (ix < 0x3fd00000) /* |x| <1/4 */
 849     {
 850       if (ix < 0x38d00000) /* |x|<2**-114 */
 851         return one - x;
 852       return one - __erfl (x);
 853     }
 854   if (ix < 0x3ff40000) /* 1.25 */
 855     {
 856       if ((hx & 0x80000000) != 0)
 857         x = -x;
 858       i = 8.0 * x;
 859       switch (i)
 860         {
 861         case 2:
 862           z = x - 0.25L;
 863           y = C13b + z * neval (z, RNr13, NRNr13) / deval (z, RDr13, NRDr13);
 864           y += C13a;
 865           break;
 866         case 3:
 867           z = x - 0.375L;
 868           y = C14b + z * neval (z, RNr14, NRNr14) / deval (z, RDr14, NRDr14);
 869           y += C14a;
 870           break;
 871         case 4:
 872           z = x - 0.5L;
 873           y = C15b + z * neval (z, RNr15, NRNr15) / deval (z, RDr15, NRDr15);
 874           y += C15a;
 875           break;
 876         case 5:
 877           z = x - 0.625L;
 878           y = C16b + z * neval (z, RNr16, NRNr16) / deval (z, RDr16, NRDr16);
 879           y += C16a;
 880           break;
 881         case 6:
 882           z = x - 0.75L;
 883           y = C17b + z * neval (z, RNr17, NRNr17) / deval (z, RDr17, NRDr17);
 884           y += C17a;
 885           break;
 886         case 7:
 887           z = x - 0.875L;
 888           y = C18b + z * neval (z, RNr18, NRNr18) / deval (z, RDr18, NRDr18);
 889           y += C18a;
 890           break;
 891         case 8:
 892           z = x - 1.0L;
 893           y = C19b + z * neval (z, RNr19, NRNr19) / deval (z, RDr19, NRDr19);
 894           y += C19a;
 895           break;
 896         default: /* i == 9.  */
 897           z = x - 1.125L;
 898           y = C20b + z * neval (z, RNr20, NRNr20) / deval (z, RDr20, NRDr20);
 899           y += C20a;
 900           break;
 901         }
 902       if (hx & 0x80000000)
 903         y = 2.0L - y;
 904       return y;
 905     }
 906   /* 1.25 < |x| < 107 */
 907   if (ix < 0x405ac000)
 908     {
 909       /* x < -9 */
 910       if (hx >= 0xc0220000)
 911         return two - tiny;
 912
 913       if ((hx & 0x80000000) != 0)
 914         x = -x;
 915       z = one / (x * x);
 916       i = 8.0 / x;
 917       switch (i)
 918         {
 919         default:
 920         case 0:
 921           p = neval (z, RNr1, NRNr1) / deval (z, RDr1, NRDr1);
 922           break;
 923         case 1:
 924           p = neval (z, RNr2, NRNr2) / deval (z, RDr2, NRDr2);
 925           break;
 926         case 2:
 927           p = neval (z, RNr3, NRNr3) / deval (z, RDr3, NRDr3);
 928           break;
 929         case 3:
 930           p = neval (z, RNr4, NRNr4) / deval (z, RDr4, NRDr4);
 931           break;
 932         case 4:
 933           p = neval (z, RNr5, NRNr5) / deval (z, RDr5, NRDr5);
 934           break;
 935         case 5:
 936           p = neval (z, RNr6, NRNr6) / deval (z, RDr6, NRDr6);
 937           break;
 938         case 6:
 939           p = neval (z, RNr7, NRNr7) / deval (z, RDr7, NRDr7);
 940           break;
 941         case 7:
 942           p = neval (z, RNr8, NRNr8) / deval (z, RDr8, NRDr8);
 943           break;
 944         }
 945       z = (float) x;
 946       r = __ieee754_expl (-z * z - 0.5625) *
 947         __ieee754_expl ((z - x) * (z + x) + p);
 948       if ((hx & 0x80000000) == 0)
 949         {
 950           long double ret = r / x;
 951           if (ret == 0)
 952             __set_errno (ERANGE);
 953           return ret;
 954         }
 955       else
 956         return two - r / x;
 957     }
 958   else
 959     {
 960       if ((hx & 0x80000000) == 0)
 961         {
 962           __set_errno (ERANGE);
 963           return tiny * tiny;
 964         }
 965       else
 966         return two - tiny;
 967     }
 968 }
 969
 970 long_double_symbol (libm, __erfcl, erfcl);