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