2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
9 * ====================================================
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
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.
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.
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/>. */
33 /* double erf(double x)
34 * double erfc(double x)
37 * erf(x) = --------- | exp(-t*t)dt
44 * erfc(-x) = 2 - erfc(x)
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 + ....)
51 * 2/sqrt(pi) = 1.128379167095512573896158903121545171688
54 * 1a. erf(x) = 1 - erfc(x), for |x| > 1.0
55 * erfc(x) = 1 - erf(x) if |x| < 1/4
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]
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
70 * 4. For x in [5/4, 107],
71 * erfc(x) = (1/x)*exp(-x*x-0.5625 + R(z))
73 * The interval is partitioned into several segments
74 * of width 1/8 in 1/x.
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);
83 * Here 4 and 5 make use of the asymptotic series
85 * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
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
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
102 #include <math_private.h>
104 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
107 neval (_Float128 x
, const _Float128
*p
, int n
)
122 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
125 deval (_Float128 x
, const _Float128
*p
, int n
)
141 static const _Float128
146 efx
= L(1.2837916709551257389615890312154517168810E-1);
149 /* erf(x) = x + x R(x^2)
151 Peak relative error 1.8e-35 */
153 static const _Float128 TN1
[NTN1
+ 1] =
155 L(-3.858252324254637124543172907442106422373E10
),
156 L(9.580319248590464682316366876952214879858E10
),
157 L(1.302170519734879977595901236693040544854E10
),
158 L(2.922956950426397417800321486727032845006E9
),
159 L(1.764317520783319397868923218385468729799E8
),
160 L(1.573436014601118630105796794840834145120E7
),
161 L(4.028077380105721388745632295157816229289E5
),
162 L(1.644056806467289066852135096352853491530E4
),
163 L(3.390868480059991640235675479463287886081E1
)
166 static const _Float128 TD1
[NTD1
+ 1] =
168 L(-3.005357030696532927149885530689529032152E11
),
169 L(-1.342602283126282827411658673839982164042E11
),
170 L(-2.777153893355340961288511024443668743399E10
),
171 L(-3.483826391033531996955620074072768276974E9
),
172 L(-2.906321047071299585682722511260895227921E8
),
173 L(-1.653347985722154162439387878512427542691E7
),
174 L(-6.245520581562848778466500301865173123136E5
),
175 L(-1.402124304177498828590239373389110545142E4
),
176 L(-1.209368072473510674493129989468348633579E2
)
181 /* erf(z+1) = erf_const + P(z)/Q(z)
183 Peak relative error 7.3e-36 */
184 static const _Float128 erf_const
= L(0.845062911510467529296875);
186 static const _Float128 TN2
[NTN2
+ 1] =
188 L(-4.088889697077485301010486931817357000235E1
),
189 L(7.157046430681808553842307502826960051036E3
),
190 L(-2.191561912574409865550015485451373731780E3
),
191 L(2.180174916555316874988981177654057337219E3
),
192 L(2.848578658049670668231333682379720943455E2
),
193 L(1.630362490952512836762810462174798925274E2
),
194 L(6.317712353961866974143739396865293596895E0
),
195 L(2.450441034183492434655586496522857578066E1
),
196 L(5.127662277706787664956025545897050896203E-1)
199 static const _Float128 TD2
[NTD2
+ 1] =
201 L(1.731026445926834008273768924015161048885E4
),
202 L(1.209682239007990370796112604286048173750E4
),
203 L(1.160950290217993641320602282462976163857E4
),
204 L(5.394294645127126577825507169061355698157E3
),
205 L(2.791239340533632669442158497532521776093E3
),
206 L(8.989365571337319032943005387378993827684E2
),
207 L(2.974016493766349409725385710897298069677E2
),
208 L(6.148192754590376378740261072533527271947E1
),
209 L(1.178502892490738445655468927408440847480E1
)
214 /* erfc(x + 0.25) = erfc(0.25) + x R(x)
216 Peak relative error 1.4e-35 */
218 static const _Float128 RNr13
[NRNr13
+ 1] =
220 L(-2.353707097641280550282633036456457014829E3
),
221 L(3.871159656228743599994116143079870279866E2
),
222 L(-3.888105134258266192210485617504098426679E2
),
223 L(-2.129998539120061668038806696199343094971E1
),
224 L(-8.125462263594034672468446317145384108734E1
),
225 L(8.151549093983505810118308635926270319660E0
),
226 L(-5.033362032729207310462422357772568553670E0
),
227 L(-4.253956621135136090295893547735851168471E-2),
228 L(-8.098602878463854789780108161581050357814E-2)
231 static const _Float128 RDr13
[NRDr13
+ 1] =
233 L(2.220448796306693503549505450626652881752E3
),
234 L(1.899133258779578688791041599040951431383E2
),
235 L(1.061906712284961110196427571557149268454E3
),
236 L(7.497086072306967965180978101974566760042E1
),
237 L(2.146796115662672795876463568170441327274E2
),
238 L(1.120156008362573736664338015952284925592E1
),
239 L(2.211014952075052616409845051695042741074E1
),
240 L(6.469655675326150785692908453094054988938E-1)
243 /* erfc(0.25) = C13a + C13b to extra precision. */
244 static const _Float128 C13a
= L(0.723663330078125);
245 static const _Float128 C13b
= L(1.0279753638067014931732235184287934646022E-5);
248 /* erfc(x + 0.375) = erfc(0.375) + x R(x)
250 Peak relative error 1.2e-35 */
252 static const _Float128 RNr14
[NRNr14
+ 1] =
254 L(-2.446164016404426277577283038988918202456E3
),
255 L(6.718753324496563913392217011618096698140E2
),
256 L(-4.581631138049836157425391886957389240794E2
),
257 L(-2.382844088987092233033215402335026078208E1
),
258 L(-7.119237852400600507927038680970936336458E1
),
259 L(1.313609646108420136332418282286454287146E1
),
260 L(-6.188608702082264389155862490056401365834E0
),
261 L(-2.787116601106678287277373011101132659279E-2),
262 L(-2.230395570574153963203348263549700967918E-2)
265 static const _Float128 RDr14
[NRDr14
+ 1] =
267 L(2.495187439241869732696223349840963702875E3
),
268 L(2.503549449872925580011284635695738412162E2
),
269 L(1.159033560988895481698051531263861842461E3
),
270 L(9.493751466542304491261487998684383688622E1
),
271 L(2.276214929562354328261422263078480321204E2
),
272 L(1.367697521219069280358984081407807931847E1
),
273 L(2.276988395995528495055594829206582732682E1
),
274 L(7.647745753648996559837591812375456641163E-1)
277 /* erfc(0.375) = C14a + C14b to extra precision. */
278 static const _Float128 C14a
= L(0.5958709716796875);
279 static const _Float128 C14b
= L(1.2118885490201676174914080878232469565953E-5);
281 /* erfc(x + 0.5) = erfc(0.5) + x R(x)
283 Peak relative error 4.7e-36 */
285 static const _Float128 RNr15
[NRNr15
+ 1] =
287 L(-2.624212418011181487924855581955853461925E3
),
288 L(8.473828904647825181073831556439301342756E2
),
289 L(-5.286207458628380765099405359607331669027E2
),
290 L(-3.895781234155315729088407259045269652318E1
),
291 L(-6.200857908065163618041240848728398496256E1
),
292 L(1.469324610346924001393137895116129204737E1
),
293 L(-6.961356525370658572800674953305625578903E0
),
294 L(5.145724386641163809595512876629030548495E-3),
295 L(1.990253655948179713415957791776180406812E-2)
298 static const _Float128 RDr15
[NRDr15
+ 1] =
300 L(2.986190760847974943034021764693341524962E3
),
301 L(5.288262758961073066335410218650047725985E2
),
302 L(1.363649178071006978355113026427856008978E3
),
303 L(1.921707975649915894241864988942255320833E2
),
304 L(2.588651100651029023069013885900085533226E2
),
305 L(2.628752920321455606558942309396855629459E1
),
306 L(2.455649035885114308978333741080991380610E1
),
307 L(1.378826653595128464383127836412100939126E0
)
310 /* erfc(0.5) = C15a + C15b to extra precision. */
311 static const _Float128 C15a
= L(0.4794921875);
312 static const _Float128 C15b
= L(7.9346869534623172533461080354712635484242E-6);
314 /* erfc(x + 0.625) = erfc(0.625) + x R(x)
316 Peak relative error 5.1e-36 */
318 static const _Float128 RNr16
[NRNr16
+ 1] =
320 L(-2.347887943200680563784690094002722906820E3
),
321 L(8.008590660692105004780722726421020136482E2
),
322 L(-5.257363310384119728760181252132311447963E2
),
323 L(-4.471737717857801230450290232600243795637E1
),
324 L(-4.849540386452573306708795324759300320304E1
),
325 L(1.140885264677134679275986782978655952843E1
),
326 L(-6.731591085460269447926746876983786152300E0
),
327 L(1.370831653033047440345050025876085121231E-1),
328 L(2.022958279982138755020825717073966576670E-2),
331 static const _Float128 RDr16
[NRDr16
+ 1] =
333 L(3.075166170024837215399323264868308087281E3
),
334 L(8.730468942160798031608053127270430036627E2
),
335 L(1.458472799166340479742581949088453244767E3
),
336 L(3.230423687568019709453130785873540386217E2
),
337 L(2.804009872719893612081109617983169474655E2
),
338 L(4.465334221323222943418085830026979293091E1
),
339 L(2.612723259683205928103787842214809134746E1
),
340 L(2.341526751185244109722204018543276124997E0
),
343 /* erfc(0.625) = C16a + C16b to extra precision. */
344 static const _Float128 C16a
= L(0.3767547607421875);
345 static const _Float128 C16b
= L(4.3570693945275513594941232097252997287766E-6);
347 /* erfc(x + 0.75) = erfc(0.75) + x R(x)
349 Peak relative error 1.7e-35 */
351 static const _Float128 RNr17
[NRNr17
+ 1] =
353 L(-1.767068734220277728233364375724380366826E3
),
354 L(6.693746645665242832426891888805363898707E2
),
355 L(-4.746224241837275958126060307406616817753E2
),
356 L(-2.274160637728782675145666064841883803196E1
),
357 L(-3.541232266140939050094370552538987982637E1
),
358 L(6.988950514747052676394491563585179503865E0
),
359 L(-5.807687216836540830881352383529281215100E0
),
360 L(3.631915988567346438830283503729569443642E-1),
361 L(-1.488945487149634820537348176770282391202E-2)
364 static const _Float128 RDr17
[NRDr17
+ 1] =
366 L(2.748457523498150741964464942246913394647E3
),
367 L(1.020213390713477686776037331757871252652E3
),
368 L(1.388857635935432621972601695296561952738E3
),
369 L(3.903363681143817750895999579637315491087E2
),
370 L(2.784568344378139499217928969529219886578E2
),
371 L(5.555800830216764702779238020065345401144E1
),
372 L(2.646215470959050279430447295801291168941E1
),
373 L(2.984905282103517497081766758550112011265E0
),
376 /* erfc(0.75) = C17a + C17b to extra precision. */
377 static const _Float128 C17a
= L(0.2888336181640625);
378 static const _Float128 C17b
= L(1.0748182422368401062165408589222625794046E-5);
381 /* erfc(x + 0.875) = erfc(0.875) + x R(x)
383 Peak relative error 2.2e-35 */
385 static const _Float128 RNr18
[NRNr18
+ 1] =
387 L(-1.342044899087593397419622771847219619588E3
),
388 L(6.127221294229172997509252330961641850598E2
),
389 L(-4.519821356522291185621206350470820610727E2
),
390 L(1.223275177825128732497510264197915160235E1
),
391 L(-2.730789571382971355625020710543532867692E1
),
392 L(4.045181204921538886880171727755445395862E0
),
393 L(-4.925146477876592723401384464691452700539E0
),
394 L(5.933878036611279244654299924101068088582E-1),
395 L(-5.557645435858916025452563379795159124753E-2)
398 static const _Float128 RDr18
[NRDr18
+ 1] =
400 L(2.557518000661700588758505116291983092951E3
),
401 L(1.070171433382888994954602511991940418588E3
),
402 L(1.344842834423493081054489613250688918709E3
),
403 L(4.161144478449381901208660598266288188426E2
),
404 L(2.763670252219855198052378138756906980422E2
),
405 L(5.998153487868943708236273854747564557632E1
),
406 L(2.657695108438628847733050476209037025318E1
),
407 L(3.252140524394421868923289114410336976512E0
),
410 /* erfc(0.875) = C18a + C18b to extra precision. */
411 static const _Float128 C18a
= L(0.215911865234375);
412 static const _Float128 C18b
= L(1.3073705765341685464282101150637224028267E-5);
414 /* erfc(x + 1.0) = erfc(1.0) + x R(x)
416 Peak relative error 1.6e-35 */
418 static const _Float128 RNr19
[NRNr19
+ 1] =
420 L(-1.139180936454157193495882956565663294826E3
),
421 L(6.134903129086899737514712477207945973616E2
),
422 L(-4.628909024715329562325555164720732868263E2
),
423 L(4.165702387210732352564932347500364010833E1
),
424 L(-2.286979913515229747204101330405771801610E1
),
425 L(1.870695256449872743066783202326943667722E0
),
426 L(-4.177486601273105752879868187237000032364E0
),
427 L(7.533980372789646140112424811291782526263E-1),
428 L(-8.629945436917752003058064731308767664446E-2)
431 static const _Float128 RDr19
[NRDr19
+ 1] =
433 L(2.744303447981132701432716278363418643778E3
),
434 L(1.266396359526187065222528050591302171471E3
),
435 L(1.466739461422073351497972255511919814273E3
),
436 L(4.868710570759693955597496520298058147162E2
),
437 L(2.993694301559756046478189634131722579643E2
),
438 L(6.868976819510254139741559102693828237440E1
),
439 L(2.801505816247677193480190483913753613630E1
),
440 L(3.604439909194350263552750347742663954481E0
),
443 /* erfc(1.0) = C19a + C19b to extra precision. */
444 static const _Float128 C19a
= L(0.15728759765625);
445 static const _Float128 C19b
= L(1.1609394035130658779364917390740703933002E-5);
447 /* erfc(x + 1.125) = erfc(1.125) + x R(x)
449 Peak relative error 3.6e-36 */
451 static const _Float128 RNr20
[NRNr20
+ 1] =
453 L(-9.652706916457973956366721379612508047640E2
),
454 L(5.577066396050932776683469951773643880634E2
),
455 L(-4.406335508848496713572223098693575485978E2
),
456 L(5.202893466490242733570232680736966655434E1
),
457 L(-1.931311847665757913322495948705563937159E1
),
458 L(-9.364318268748287664267341457164918090611E-2),
459 L(-3.306390351286352764891355375882586201069E0
),
460 L(7.573806045289044647727613003096916516475E-1),
461 L(-9.611744011489092894027478899545635991213E-2)
464 static const _Float128 RDr20
[NRDr20
+ 1] =
466 L(3.032829629520142564106649167182428189014E3
),
467 L(1.659648470721967719961167083684972196891E3
),
468 L(1.703545128657284619402511356932569292535E3
),
469 L(6.393465677731598872500200253155257708763E2
),
470 L(3.489131397281030947405287112726059221934E2
),
471 L(8.848641738570783406484348434387611713070E1
),
472 L(3.132269062552392974833215844236160958502E1
),
473 L(4.430131663290563523933419966185230513168E0
)
476 /* erfc(1.125) = C20a + C20b to extra precision. */
477 static const _Float128 C20a
= L(0.111602783203125);
478 static const _Float128 C20b
= L(8.9850951672359304215530728365232161564636E-6);
480 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
482 Peak relative error 1.4e-35 */
484 static const _Float128 RNr8
[NRNr8
+ 1] =
486 L(3.587451489255356250759834295199296936784E1
),
487 L(5.406249749087340431871378009874875889602E2
),
488 L(2.931301290625250886238822286506381194157E3
),
489 L(7.359254185241795584113047248898753470923E3
),
490 L(9.201031849810636104112101947312492532314E3
),
491 L(5.749697096193191467751650366613289284777E3
),
492 L(1.710415234419860825710780802678697889231E3
),
493 L(2.150753982543378580859546706243022719599E2
),
494 L(8.740953582272147335100537849981160931197E0
),
495 L(4.876422978828717219629814794707963640913E-2)
498 static const _Float128 RDr8
[NRDr8
+ 1] =
500 L(6.358593134096908350929496535931630140282E1
),
501 L(9.900253816552450073757174323424051765523E2
),
502 L(5.642928777856801020545245437089490805186E3
),
503 L(1.524195375199570868195152698617273739609E4
),
504 L(2.113829644500006749947332935305800887345E4
),
505 L(1.526438562626465706267943737310282977138E4
),
506 L(5.561370922149241457131421914140039411782E3
),
507 L(9.394035530179705051609070428036834496942E2
),
508 L(6.147019596150394577984175188032707343615E1
)
512 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
514 Peak relative error 2.0e-36 */
516 static const _Float128 RNr7
[NRNr7
+ 1] =
518 L(1.686222193385987690785945787708644476545E1
),
519 L(1.178224543567604215602418571310612066594E3
),
520 L(1.764550584290149466653899886088166091093E4
),
521 L(1.073758321890334822002849369898232811561E5
),
522 L(3.132840749205943137619839114451290324371E5
),
523 L(4.607864939974100224615527007793867585915E5
),
524 L(3.389781820105852303125270837910972384510E5
),
525 L(1.174042187110565202875011358512564753399E5
),
526 L(1.660013606011167144046604892622504338313E4
),
527 L(6.700393957480661937695573729183733234400E2
)
530 static const _Float128 RDr7
[NRDr7
+ 1] =
532 L(-1.709305024718358874701575813642933561169E3
),
533 L(-3.280033887481333199580464617020514788369E4
),
534 L(-2.345284228022521885093072363418750835214E5
),
535 L(-8.086758123097763971926711729242327554917E5
),
536 L(-1.456900414510108718402423999575992450138E6
),
537 L(-1.391654264881255068392389037292702041855E6
),
538 L(-6.842360801869939983674527468509852583855E5
),
539 L(-1.597430214446573566179675395199807533371E5
),
540 L(-1.488876130609876681421645314851760773480E4
),
541 L(-3.511762950935060301403599443436465645703E2
)
545 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
547 Peak relative error 1.9e-35 */
549 static const _Float128 RNr6
[NRNr6
+ 1] =
551 L(1.642076876176834390623842732352935761108E0
),
552 L(1.207150003611117689000664385596211076662E2
),
553 L(2.119260779316389904742873816462800103939E3
),
554 L(1.562942227734663441801452930916044224174E4
),
555 L(5.656779189549710079988084081145693580479E4
),
556 L(1.052166241021481691922831746350942786299E5
),
557 L(9.949798524786000595621602790068349165758E4
),
558 L(4.491790734080265043407035220188849562856E4
),
559 L(8.377074098301530326270432059434791287601E3
),
560 L(4.506934806567986810091824791963991057083E2
)
563 static const _Float128 RDr6
[NRDr6
+ 1] =
565 L(-1.664557643928263091879301304019826629067E2
),
566 L(-3.800035902507656624590531122291160668452E3
),
567 L(-3.277028191591734928360050685359277076056E4
),
568 L(-1.381359471502885446400589109566587443987E5
),
569 L(-3.082204287382581873532528989283748656546E5
),
570 L(-3.691071488256738343008271448234631037095E5
),
571 L(-2.300482443038349815750714219117566715043E5
),
572 L(-6.873955300927636236692803579555752171530E4
),
573 L(-8.262158817978334142081581542749986845399E3
),
574 L(-2.517122254384430859629423488157361983661E2
)
578 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
580 Peak relative error 4.6e-36 */
582 static const _Float128 RNr5
[NRNr5
+ 1] =
584 L(-3.332258927455285458355550878136506961608E-3),
585 L(-2.697100758900280402659586595884478660721E-1),
586 L(-6.083328551139621521416618424949137195536E0
),
587 L(-6.119863528983308012970821226810162441263E1
),
588 L(-3.176535282475593173248810678636522589861E2
),
589 L(-8.933395175080560925809992467187963260693E2
),
590 L(-1.360019508488475978060917477620199499560E3
),
591 L(-1.075075579828188621541398761300910213280E3
),
592 L(-4.017346561586014822824459436695197089916E2
),
593 L(-5.857581368145266249509589726077645791341E1
),
594 L(-2.077715925587834606379119585995758954399E0
)
597 static const _Float128 RDr5
[NRDr5
+ 1] =
599 L(3.377879570417399341550710467744693125385E-1),
600 L(1.021963322742390735430008860602594456187E1
),
601 L(1.200847646592942095192766255154827011939E2
),
602 L(7.118915528142927104078182863387116942836E2
),
603 L(2.318159380062066469386544552429625026238E3
),
604 L(4.238729853534009221025582008928765281620E3
),
605 L(4.279114907284825886266493994833515580782E3
),
606 L(2.257277186663261531053293222591851737504E3
),
607 L(5.570475501285054293371908382916063822957E2
),
608 L(5.142189243856288981145786492585432443560E1
)
612 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
614 Peak relative error 2.0e-36 */
616 static const _Float128 RNr4
[NRNr4
+ 1] =
618 L(3.258530712024527835089319075288494524465E-3),
619 L(2.987056016877277929720231688689431056567E-1),
620 L(8.738729089340199750734409156830371528862E0
),
621 L(1.207211160148647782396337792426311125923E2
),
622 L(8.997558632489032902250523945248208224445E2
),
623 L(3.798025197699757225978410230530640879762E3
),
624 L(9.113203668683080975637043118209210146846E3
),
625 L(1.203285891339933238608683715194034900149E4
),
626 L(8.100647057919140328536743641735339740855E3
),
627 L(2.383888249907144945837976899822927411769E3
),
628 L(2.127493573166454249221983582495245662319E2
)
631 static const _Float128 RDr4
[NRDr4
+ 1] =
633 L(-3.303141981514540274165450687270180479586E-1),
634 L(-1.353768629363605300707949368917687066724E1
),
635 L(-2.206127630303621521950193783894598987033E2
),
636 L(-1.861800338758066696514480386180875607204E3
),
637 L(-8.889048775872605708249140016201753255599E3
),
638 L(-2.465888106627948210478692168261494857089E4
),
639 L(-3.934642211710774494879042116768390014289E4
),
640 L(-3.455077258242252974937480623730228841003E4
),
641 L(-1.524083977439690284820586063729912653196E4
),
642 L(-2.810541887397984804237552337349093953857E3
),
643 L(-1.343929553541159933824901621702567066156E2
)
647 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
649 Peak relative error 8.4e-37 */
651 static const _Float128 RNr3
[NRNr3
+ 1] =
653 L(-1.952401126551202208698629992497306292987E-6),
654 L(-2.130881743066372952515162564941682716125E-4),
655 L(-8.376493958090190943737529486107282224387E-3),
656 L(-1.650592646560987700661598877522831234791E-1),
657 L(-1.839290818933317338111364667708678163199E0
),
658 L(-1.216278715570882422410442318517814388470E1
),
659 L(-4.818759344462360427612133632533779091386E1
),
660 L(-1.120994661297476876804405329172164436784E2
),
661 L(-1.452850765662319264191141091859300126931E2
),
662 L(-9.485207851128957108648038238656777241333E1
),
663 L(-2.563663855025796641216191848818620020073E1
),
664 L(-1.787995944187565676837847610706317833247E0
)
667 static const _Float128 RDr3
[NRDr3
+ 1] =
669 L(1.979130686770349481460559711878399476903E-4),
670 L(1.156941716128488266238105813374635099057E-2),
671 L(2.752657634309886336431266395637285974292E-1),
672 L(3.482245457248318787349778336603569327521E0
),
673 L(2.569347069372696358578399521203959253162E1
),
674 L(1.142279000180457419740314694631879921561E2
),
675 L(3.056503977190564294341422623108332700840E2
),
676 L(4.780844020923794821656358157128719184422E2
),
677 L(4.105972727212554277496256802312730410518E2
),
678 L(1.724072188063746970865027817017067646246E2
),
679 L(2.815939183464818198705278118326590370435E1
)
683 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
685 Peak relative error 1.5e-36 */
687 static const _Float128 RNr2
[NRNr2
+ 1] =
689 L(-2.638914383420287212401687401284326363787E-8),
690 L(-3.479198370260633977258201271399116766619E-6),
691 L(-1.783985295335697686382487087502222519983E-4),
692 L(-4.777876933122576014266349277217559356276E-3),
693 L(-7.450634738987325004070761301045014986520E-2),
694 L(-7.068318854874733315971973707247467326619E-1),
695 L(-4.113919921935944795764071670806867038732E0
),
696 L(-1.440447573226906222417767283691888875082E1
),
697 L(-2.883484031530718428417168042141288943905E1
),
698 L(-2.990886974328476387277797361464279931446E1
),
699 L(-1.325283914915104866248279787536128997331E1
),
700 L(-1.572436106228070195510230310658206154374E0
)
703 static const _Float128 RDr2
[NRDr2
+ 1] =
705 L(2.675042728136731923554119302571867799673E-6),
706 L(2.170997868451812708585443282998329996268E-4),
707 L(7.249969752687540289422684951196241427445E-3),
708 L(1.302040375859768674620410563307838448508E-1),
709 L(1.380202483082910888897654537144485285549E0
),
710 L(8.926594113174165352623847870299170069350E0
),
711 L(3.521089584782616472372909095331572607185E1
),
712 L(8.233547427533181375185259050330809105570E1
),
713 L(1.072971579885803033079469639073292840135E2
),
714 L(6.943803113337964469736022094105143158033E1
),
715 L(1.775695341031607738233608307835017282662E1
)
719 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
721 Peak relative error 2.2e-36 */
723 static const _Float128 RNr1
[NRNr1
+ 1] =
725 L(-4.250780883202361946697751475473042685782E-8),
726 L(-5.375777053288612282487696975623206383019E-6),
727 L(-2.573645949220896816208565944117382460452E-4),
728 L(-6.199032928113542080263152610799113086319E-3),
729 L(-8.262721198693404060380104048479916247786E-2),
730 L(-6.242615227257324746371284637695778043982E-1),
731 L(-2.609874739199595400225113299437099626386E0
),
732 L(-5.581967563336676737146358534602770006970E0
),
733 L(-5.124398923356022609707490956634280573882E0
),
734 L(-1.290865243944292370661544030414667556649E0
)
737 static const _Float128 RDr1
[NRDr1
+ 1] =
739 L(4.308976661749509034845251315983612976224E-6),
740 L(3.265390126432780184125233455960049294580E-4),
741 L(9.811328839187040701901866531796570418691E-3),
742 L(1.511222515036021033410078631914783519649E-1),
743 L(1.289264341917429958858379585970225092274E0
),
744 L(6.147640356182230769548007536914983522270E0
),
745 L(1.573966871337739784518246317003956180750E1
),
746 L(1.955534123435095067199574045529218238263E1
),
747 L(9.472613121363135472247929109615785855865E0
)
757 ieee854_long_double_shape_type u
;
761 ix
= sign
& 0x7fffffff;
763 if (ix
>= 0x7fff0000)
765 i
= ((sign
& 0xffff0000) >> 31) << 1;
766 return (_Float128
) (1 - i
) + one
/ x
; /* erf(+-inf)=+-1 */
769 if (ix
>= 0x3fff0000) /* |x| >= 1.0 */
771 if (ix
>= 0x40030000 && sign
> 0)
772 return one
; /* x >= 16, avoid spurious underflow from erfc. */
775 /* return (one - __erfcl (x)); */
780 if (ix
< 0x3ffec000) /* a < 0.875 */
782 if (ix
< 0x3fc60000) /* |x|<2**-57 */
786 /* Avoid spurious underflow. */
787 _Float128 ret
= 0.0625 * (16.0 * x
+ (16.0 * efx
) * x
);
788 math_check_force_underflow (ret
);
793 y
= a
+ a
* neval (z
, TN1
, NTN1
) / deval (z
, TD1
, NTD1
);
798 y
= erf_const
+ neval (a
, TN2
, NTN2
) / deval (a
, TD2
, NTD2
);
801 if (sign
& 0x80000000) /* x < 0 */
806 weak_alias (__erfl
, erfl
)
808 __erfcl (_Float128 x
)
810 _Float128 y
, z
, p
, r
;
812 ieee854_long_double_shape_type u
;
816 ix
= sign
& 0x7fffffff;
819 if (ix
>= 0x7fff0000)
820 { /* erfc(nan)=nan */
821 /* erfc(+-inf)=0,2 */
822 return (_Float128
) (((u_int32_t
) sign
>> 31) << 1) + one
/ x
;
825 if (ix
< 0x3ffd0000) /* |x| <1/4 */
827 if (ix
< 0x3f8d0000) /* |x|<2**-114 */
829 return one
- __erfl (x
);
831 if (ix
< 0x3fff4000) /* 1.25 */
839 y
= C13b
+ z
* neval (z
, RNr13
, NRNr13
) / deval (z
, RDr13
, NRDr13
);
844 y
= C14b
+ z
* neval (z
, RNr14
, NRNr14
) / deval (z
, RDr14
, NRDr14
);
849 y
= C15b
+ z
* neval (z
, RNr15
, NRNr15
) / deval (z
, RDr15
, NRDr15
);
854 y
= C16b
+ z
* neval (z
, RNr16
, NRNr16
) / deval (z
, RDr16
, NRDr16
);
859 y
= C17b
+ z
* neval (z
, RNr17
, NRNr17
) / deval (z
, RDr17
, NRDr17
);
864 y
= C18b
+ z
* neval (z
, RNr18
, NRNr18
) / deval (z
, RDr18
, NRDr18
);
869 y
= C19b
+ z
* neval (z
, RNr19
, NRNr19
) / deval (z
, RDr19
, NRDr19
);
872 default: /* i == 9. */
874 y
= C20b
+ z
* neval (z
, RNr20
, NRNr20
) / deval (z
, RDr20
, NRDr20
);
878 if (sign
& 0x80000000)
882 /* 1.25 < |x| < 107 */
886 if ((ix
>= 0x40022000) && (sign
& 0x80000000))
896 p
= neval (z
, RNr1
, NRNr1
) / deval (z
, RDr1
, NRDr1
);
899 p
= neval (z
, RNr2
, NRNr2
) / deval (z
, RDr2
, NRDr2
);
902 p
= neval (z
, RNr3
, NRNr3
) / deval (z
, RDr3
, NRDr3
);
905 p
= neval (z
, RNr4
, NRNr4
) / deval (z
, RDr4
, NRDr4
);
908 p
= neval (z
, RNr5
, NRNr5
) / deval (z
, RDr5
, NRDr5
);
911 p
= neval (z
, RNr6
, NRNr6
) / deval (z
, RDr6
, NRDr6
);
914 p
= neval (z
, RNr7
, NRNr7
) / deval (z
, RDr7
, NRDr7
);
917 p
= neval (z
, RNr8
, NRNr8
) / deval (z
, RDr8
, NRDr8
);
922 u
.parts32
.w2
&= 0xfe000000;
924 r
= __ieee754_expl (-z
* z
- 0.5625) *
925 __ieee754_expl ((z
- x
) * (z
+ x
) + p
);
926 if ((sign
& 0x80000000) == 0)
928 _Float128 ret
= r
/ x
;
930 __set_errno (ERANGE
);
938 if ((sign
& 0x80000000) == 0)
940 __set_errno (ERANGE
);
948 weak_alias (__erfcl
, erfcl
)