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