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