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