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
100 #include "math_private.h"
101 #include <math_ldbl_opt.h>
103 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
106 neval (long double x
, const long double *p
, int n
)
121 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
124 deval (long double x
, const long double *p
, int n
)
140 static const long double
146 efx
= 1.2837916709551257389615890312154517168810E-1L,
147 /* 8 * (2/sqrt(pi) - 1) */
148 efx8
= 1.0270333367641005911692712249723613735048E0L
;
151 /* erf(x) = x + x R(x^2)
153 Peak relative error 1.8e-35 */
155 static const long double TN1
[NTN1
+ 1] =
157 -3.858252324254637124543172907442106422373E10L
,
158 9.580319248590464682316366876952214879858E10L
,
159 1.302170519734879977595901236693040544854E10L
,
160 2.922956950426397417800321486727032845006E9L
,
161 1.764317520783319397868923218385468729799E8L
,
162 1.573436014601118630105796794840834145120E7L
,
163 4.028077380105721388745632295157816229289E5L
,
164 1.644056806467289066852135096352853491530E4L
,
165 3.390868480059991640235675479463287886081E1L
168 static const long double TD1
[NTD1
+ 1] =
170 -3.005357030696532927149885530689529032152E11L
,
171 -1.342602283126282827411658673839982164042E11L
,
172 -2.777153893355340961288511024443668743399E10L
,
173 -3.483826391033531996955620074072768276974E9L
,
174 -2.906321047071299585682722511260895227921E8L
,
175 -1.653347985722154162439387878512427542691E7L
,
176 -6.245520581562848778466500301865173123136E5L
,
177 -1.402124304177498828590239373389110545142E4L
,
178 -1.209368072473510674493129989468348633579E2L
183 /* erf(z+1) = erf_const + P(z)/Q(z)
185 Peak relative error 7.3e-36 */
186 static const long double erf_const
= 0.845062911510467529296875L;
188 static const long double TN2
[NTN2
+ 1] =
190 -4.088889697077485301010486931817357000235E1L
,
191 7.157046430681808553842307502826960051036E3L
,
192 -2.191561912574409865550015485451373731780E3L
,
193 2.180174916555316874988981177654057337219E3L
,
194 2.848578658049670668231333682379720943455E2L
,
195 1.630362490952512836762810462174798925274E2L
,
196 6.317712353961866974143739396865293596895E0L
,
197 2.450441034183492434655586496522857578066E1L
,
198 5.127662277706787664956025545897050896203E-1L
201 static const long double TD2
[NTD2
+ 1] =
203 1.731026445926834008273768924015161048885E4L
,
204 1.209682239007990370796112604286048173750E4L
,
205 1.160950290217993641320602282462976163857E4L
,
206 5.394294645127126577825507169061355698157E3L
,
207 2.791239340533632669442158497532521776093E3L
,
208 8.989365571337319032943005387378993827684E2L
,
209 2.974016493766349409725385710897298069677E2L
,
210 6.148192754590376378740261072533527271947E1L
,
211 1.178502892490738445655468927408440847480E1L
216 /* erfc(x + 0.25) = erfc(0.25) + x R(x)
218 Peak relative error 1.4e-35 */
220 static const long double RNr13
[NRNr13
+ 1] =
222 -2.353707097641280550282633036456457014829E3L
,
223 3.871159656228743599994116143079870279866E2L
,
224 -3.888105134258266192210485617504098426679E2L
,
225 -2.129998539120061668038806696199343094971E1L
,
226 -8.125462263594034672468446317145384108734E1L
,
227 8.151549093983505810118308635926270319660E0L
,
228 -5.033362032729207310462422357772568553670E0L
,
229 -4.253956621135136090295893547735851168471E-2L,
230 -8.098602878463854789780108161581050357814E-2L
233 static const long double RDr13
[NRDr13
+ 1] =
235 2.220448796306693503549505450626652881752E3L
,
236 1.899133258779578688791041599040951431383E2L
,
237 1.061906712284961110196427571557149268454E3L
,
238 7.497086072306967965180978101974566760042E1L
,
239 2.146796115662672795876463568170441327274E2L
,
240 1.120156008362573736664338015952284925592E1L
,
241 2.211014952075052616409845051695042741074E1L
,
242 6.469655675326150785692908453094054988938E-1L
245 /* erfc(0.25) = C13a + C13b to extra precision. */
246 static const long double C13a
= 0.723663330078125L;
247 static const long double C13b
= 1.0279753638067014931732235184287934646022E-5L;
250 /* erfc(x + 0.375) = erfc(0.375) + x R(x)
252 Peak relative error 1.2e-35 */
254 static const long double RNr14
[NRNr14
+ 1] =
256 -2.446164016404426277577283038988918202456E3L
,
257 6.718753324496563913392217011618096698140E2L
,
258 -4.581631138049836157425391886957389240794E2L
,
259 -2.382844088987092233033215402335026078208E1L
,
260 -7.119237852400600507927038680970936336458E1L
,
261 1.313609646108420136332418282286454287146E1L
,
262 -6.188608702082264389155862490056401365834E0L
,
263 -2.787116601106678287277373011101132659279E-2L,
264 -2.230395570574153963203348263549700967918E-2L
267 static const long double RDr14
[NRDr14
+ 1] =
269 2.495187439241869732696223349840963702875E3L
,
270 2.503549449872925580011284635695738412162E2L
,
271 1.159033560988895481698051531263861842461E3L
,
272 9.493751466542304491261487998684383688622E1L
,
273 2.276214929562354328261422263078480321204E2L
,
274 1.367697521219069280358984081407807931847E1L
,
275 2.276988395995528495055594829206582732682E1L
,
276 7.647745753648996559837591812375456641163E-1L
279 /* erfc(0.375) = C14a + C14b to extra precision. */
280 static const long double C14a
= 0.5958709716796875L;
281 static const long double C14b
= 1.2118885490201676174914080878232469565953E-5L;
283 /* erfc(x + 0.5) = erfc(0.5) + x R(x)
285 Peak relative error 4.7e-36 */
287 static const long double RNr15
[NRNr15
+ 1] =
289 -2.624212418011181487924855581955853461925E3L
,
290 8.473828904647825181073831556439301342756E2L
,
291 -5.286207458628380765099405359607331669027E2L
,
292 -3.895781234155315729088407259045269652318E1L
,
293 -6.200857908065163618041240848728398496256E1L
,
294 1.469324610346924001393137895116129204737E1L
,
295 -6.961356525370658572800674953305625578903E0L
,
296 5.145724386641163809595512876629030548495E-3L,
297 1.990253655948179713415957791776180406812E-2L
300 static const long double RDr15
[NRDr15
+ 1] =
302 2.986190760847974943034021764693341524962E3L
,
303 5.288262758961073066335410218650047725985E2L
,
304 1.363649178071006978355113026427856008978E3L
,
305 1.921707975649915894241864988942255320833E2L
,
306 2.588651100651029023069013885900085533226E2L
,
307 2.628752920321455606558942309396855629459E1L
,
308 2.455649035885114308978333741080991380610E1L
,
309 1.378826653595128464383127836412100939126E0L
312 /* erfc(0.5) = C15a + C15b to extra precision. */
313 static const long double C15a
= 0.4794921875L;
314 static const long double C15b
= 7.9346869534623172533461080354712635484242E-6L;
316 /* erfc(x + 0.625) = erfc(0.625) + x R(x)
318 Peak relative error 5.1e-36 */
320 static const long double RNr16
[NRNr16
+ 1] =
322 -2.347887943200680563784690094002722906820E3L
,
323 8.008590660692105004780722726421020136482E2L
,
324 -5.257363310384119728760181252132311447963E2L
,
325 -4.471737717857801230450290232600243795637E1L
,
326 -4.849540386452573306708795324759300320304E1L
,
327 1.140885264677134679275986782978655952843E1L
,
328 -6.731591085460269447926746876983786152300E0L
,
329 1.370831653033047440345050025876085121231E-1L,
330 2.022958279982138755020825717073966576670E-2L,
333 static const long double RDr16
[NRDr16
+ 1] =
335 3.075166170024837215399323264868308087281E3L
,
336 8.730468942160798031608053127270430036627E2L
,
337 1.458472799166340479742581949088453244767E3L
,
338 3.230423687568019709453130785873540386217E2L
,
339 2.804009872719893612081109617983169474655E2L
,
340 4.465334221323222943418085830026979293091E1L
,
341 2.612723259683205928103787842214809134746E1L
,
342 2.341526751185244109722204018543276124997E0L
,
345 /* erfc(0.625) = C16a + C16b to extra precision. */
346 static const long double C16a
= 0.3767547607421875L;
347 static const long double C16b
= 4.3570693945275513594941232097252997287766E-6L;
349 /* erfc(x + 0.75) = erfc(0.75) + x R(x)
351 Peak relative error 1.7e-35 */
353 static const long double RNr17
[NRNr17
+ 1] =
355 -1.767068734220277728233364375724380366826E3L
,
356 6.693746645665242832426891888805363898707E2L
,
357 -4.746224241837275958126060307406616817753E2L
,
358 -2.274160637728782675145666064841883803196E1L
,
359 -3.541232266140939050094370552538987982637E1L
,
360 6.988950514747052676394491563585179503865E0L
,
361 -5.807687216836540830881352383529281215100E0L
,
362 3.631915988567346438830283503729569443642E-1L,
363 -1.488945487149634820537348176770282391202E-2L
366 static const long double RDr17
[NRDr17
+ 1] =
368 2.748457523498150741964464942246913394647E3L
,
369 1.020213390713477686776037331757871252652E3L
,
370 1.388857635935432621972601695296561952738E3L
,
371 3.903363681143817750895999579637315491087E2L
,
372 2.784568344378139499217928969529219886578E2L
,
373 5.555800830216764702779238020065345401144E1L
,
374 2.646215470959050279430447295801291168941E1L
,
375 2.984905282103517497081766758550112011265E0L
,
378 /* erfc(0.75) = C17a + C17b to extra precision. */
379 static const long double C17a
= 0.2888336181640625L;
380 static const long double C17b
= 1.0748182422368401062165408589222625794046E-5L;
383 /* erfc(x + 0.875) = erfc(0.875) + x R(x)
385 Peak relative error 2.2e-35 */
387 static const long double RNr18
[NRNr18
+ 1] =
389 -1.342044899087593397419622771847219619588E3L
,
390 6.127221294229172997509252330961641850598E2L
,
391 -4.519821356522291185621206350470820610727E2L
,
392 1.223275177825128732497510264197915160235E1L
,
393 -2.730789571382971355625020710543532867692E1L
,
394 4.045181204921538886880171727755445395862E0L
,
395 -4.925146477876592723401384464691452700539E0L
,
396 5.933878036611279244654299924101068088582E-1L,
397 -5.557645435858916025452563379795159124753E-2L
400 static const long double RDr18
[NRDr18
+ 1] =
402 2.557518000661700588758505116291983092951E3L
,
403 1.070171433382888994954602511991940418588E3L
,
404 1.344842834423493081054489613250688918709E3L
,
405 4.161144478449381901208660598266288188426E2L
,
406 2.763670252219855198052378138756906980422E2L
,
407 5.998153487868943708236273854747564557632E1L
,
408 2.657695108438628847733050476209037025318E1L
,
409 3.252140524394421868923289114410336976512E0L
,
412 /* erfc(0.875) = C18a + C18b to extra precision. */
413 static const long double C18a
= 0.215911865234375L;
414 static const long double C18b
= 1.3073705765341685464282101150637224028267E-5L;
416 /* erfc(x + 1.0) = erfc(1.0) + x R(x)
418 Peak relative error 1.6e-35 */
420 static const long double RNr19
[NRNr19
+ 1] =
422 -1.139180936454157193495882956565663294826E3L
,
423 6.134903129086899737514712477207945973616E2L
,
424 -4.628909024715329562325555164720732868263E2L
,
425 4.165702387210732352564932347500364010833E1L
,
426 -2.286979913515229747204101330405771801610E1L
,
427 1.870695256449872743066783202326943667722E0L
,
428 -4.177486601273105752879868187237000032364E0L
,
429 7.533980372789646140112424811291782526263E-1L,
430 -8.629945436917752003058064731308767664446E-2L
433 static const long double RDr19
[NRDr19
+ 1] =
435 2.744303447981132701432716278363418643778E3L
,
436 1.266396359526187065222528050591302171471E3L
,
437 1.466739461422073351497972255511919814273E3L
,
438 4.868710570759693955597496520298058147162E2L
,
439 2.993694301559756046478189634131722579643E2L
,
440 6.868976819510254139741559102693828237440E1L
,
441 2.801505816247677193480190483913753613630E1L
,
442 3.604439909194350263552750347742663954481E0L
,
445 /* erfc(1.0) = C19a + C19b to extra precision. */
446 static const long double C19a
= 0.15728759765625L;
447 static const long double C19b
= 1.1609394035130658779364917390740703933002E-5L;
449 /* erfc(x + 1.125) = erfc(1.125) + x R(x)
451 Peak relative error 3.6e-36 */
453 static const long double RNr20
[NRNr20
+ 1] =
455 -9.652706916457973956366721379612508047640E2L
,
456 5.577066396050932776683469951773643880634E2L
,
457 -4.406335508848496713572223098693575485978E2L
,
458 5.202893466490242733570232680736966655434E1L
,
459 -1.931311847665757913322495948705563937159E1L
,
460 -9.364318268748287664267341457164918090611E-2L,
461 -3.306390351286352764891355375882586201069E0L
,
462 7.573806045289044647727613003096916516475E-1L,
463 -9.611744011489092894027478899545635991213E-2L
466 static const long double RDr20
[NRDr20
+ 1] =
468 3.032829629520142564106649167182428189014E3L
,
469 1.659648470721967719961167083684972196891E3L
,
470 1.703545128657284619402511356932569292535E3L
,
471 6.393465677731598872500200253155257708763E2L
,
472 3.489131397281030947405287112726059221934E2L
,
473 8.848641738570783406484348434387611713070E1L
,
474 3.132269062552392974833215844236160958502E1L
,
475 4.430131663290563523933419966185230513168E0L
478 /* erfc(1.125) = C20a + C20b to extra precision. */
479 static const long double C20a
= 0.111602783203125L;
480 static const long double C20b
= 8.9850951672359304215530728365232161564636E-6L;
482 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
484 Peak relative error 1.4e-35 */
486 static const long double RNr8
[NRNr8
+ 1] =
488 3.587451489255356250759834295199296936784E1L
,
489 5.406249749087340431871378009874875889602E2L
,
490 2.931301290625250886238822286506381194157E3L
,
491 7.359254185241795584113047248898753470923E3L
,
492 9.201031849810636104112101947312492532314E3L
,
493 5.749697096193191467751650366613289284777E3L
,
494 1.710415234419860825710780802678697889231E3L
,
495 2.150753982543378580859546706243022719599E2L
,
496 8.740953582272147335100537849981160931197E0L
,
497 4.876422978828717219629814794707963640913E-2L
500 static const long double RDr8
[NRDr8
+ 1] =
502 6.358593134096908350929496535931630140282E1L
,
503 9.900253816552450073757174323424051765523E2L
,
504 5.642928777856801020545245437089490805186E3L
,
505 1.524195375199570868195152698617273739609E4L
,
506 2.113829644500006749947332935305800887345E4L
,
507 1.526438562626465706267943737310282977138E4L
,
508 5.561370922149241457131421914140039411782E3L
,
509 9.394035530179705051609070428036834496942E2L
,
510 6.147019596150394577984175188032707343615E1L
514 /* erfc(1/x) = 1/x exp (-1/x^2 - 0.5625 + R(1/x^2))
516 Peak relative error 2.0e-36 */
518 static const long double RNr7
[NRNr7
+ 1] =
520 1.686222193385987690785945787708644476545E1L
,
521 1.178224543567604215602418571310612066594E3L
,
522 1.764550584290149466653899886088166091093E4L
,
523 1.073758321890334822002849369898232811561E5L
,
524 3.132840749205943137619839114451290324371E5L
,
525 4.607864939974100224615527007793867585915E5L
,
526 3.389781820105852303125270837910972384510E5L
,
527 1.174042187110565202875011358512564753399E5L
,
528 1.660013606011167144046604892622504338313E4L
,
529 6.700393957480661937695573729183733234400E2L
532 static const long double RDr7
[NRDr7
+ 1] =
534 -1.709305024718358874701575813642933561169E3L
,
535 -3.280033887481333199580464617020514788369E4L
,
536 -2.345284228022521885093072363418750835214E5L
,
537 -8.086758123097763971926711729242327554917E5L
,
538 -1.456900414510108718402423999575992450138E6L
,
539 -1.391654264881255068392389037292702041855E6L
,
540 -6.842360801869939983674527468509852583855E5L
,
541 -1.597430214446573566179675395199807533371E5L
,
542 -1.488876130609876681421645314851760773480E4L
,
543 -3.511762950935060301403599443436465645703E2L
547 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
549 Peak relative error 1.9e-35 */
551 static const long double RNr6
[NRNr6
+ 1] =
553 1.642076876176834390623842732352935761108E0L
,
554 1.207150003611117689000664385596211076662E2L
,
555 2.119260779316389904742873816462800103939E3L
,
556 1.562942227734663441801452930916044224174E4L
,
557 5.656779189549710079988084081145693580479E4L
,
558 1.052166241021481691922831746350942786299E5L
,
559 9.949798524786000595621602790068349165758E4L
,
560 4.491790734080265043407035220188849562856E4L
,
561 8.377074098301530326270432059434791287601E3L
,
562 4.506934806567986810091824791963991057083E2L
565 static const long double RDr6
[NRDr6
+ 1] =
567 -1.664557643928263091879301304019826629067E2L
,
568 -3.800035902507656624590531122291160668452E3L
,
569 -3.277028191591734928360050685359277076056E4L
,
570 -1.381359471502885446400589109566587443987E5L
,
571 -3.082204287382581873532528989283748656546E5L
,
572 -3.691071488256738343008271448234631037095E5L
,
573 -2.300482443038349815750714219117566715043E5L
,
574 -6.873955300927636236692803579555752171530E4L
,
575 -8.262158817978334142081581542749986845399E3L
,
576 -2.517122254384430859629423488157361983661E2L
580 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
582 Peak relative error 4.6e-36 */
584 static const long double RNr5
[NRNr5
+ 1] =
586 -3.332258927455285458355550878136506961608E-3L,
587 -2.697100758900280402659586595884478660721E-1L,
588 -6.083328551139621521416618424949137195536E0L
,
589 -6.119863528983308012970821226810162441263E1L
,
590 -3.176535282475593173248810678636522589861E2L
,
591 -8.933395175080560925809992467187963260693E2L
,
592 -1.360019508488475978060917477620199499560E3L
,
593 -1.075075579828188621541398761300910213280E3L
,
594 -4.017346561586014822824459436695197089916E2L
,
595 -5.857581368145266249509589726077645791341E1L
,
596 -2.077715925587834606379119585995758954399E0L
599 static const long double RDr5
[NRDr5
+ 1] =
601 3.377879570417399341550710467744693125385E-1L,
602 1.021963322742390735430008860602594456187E1L
,
603 1.200847646592942095192766255154827011939E2L
,
604 7.118915528142927104078182863387116942836E2L
,
605 2.318159380062066469386544552429625026238E3L
,
606 4.238729853534009221025582008928765281620E3L
,
607 4.279114907284825886266493994833515580782E3L
,
608 2.257277186663261531053293222591851737504E3L
,
609 5.570475501285054293371908382916063822957E2L
,
610 5.142189243856288981145786492585432443560E1L
614 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
616 Peak relative error 2.0e-36 */
618 static const long double RNr4
[NRNr4
+ 1] =
620 3.258530712024527835089319075288494524465E-3L,
621 2.987056016877277929720231688689431056567E-1L,
622 8.738729089340199750734409156830371528862E0L
,
623 1.207211160148647782396337792426311125923E2L
,
624 8.997558632489032902250523945248208224445E2L
,
625 3.798025197699757225978410230530640879762E3L
,
626 9.113203668683080975637043118209210146846E3L
,
627 1.203285891339933238608683715194034900149E4L
,
628 8.100647057919140328536743641735339740855E3L
,
629 2.383888249907144945837976899822927411769E3L
,
630 2.127493573166454249221983582495245662319E2L
633 static const long double RDr4
[NRDr4
+ 1] =
635 -3.303141981514540274165450687270180479586E-1L,
636 -1.353768629363605300707949368917687066724E1L
,
637 -2.206127630303621521950193783894598987033E2L
,
638 -1.861800338758066696514480386180875607204E3L
,
639 -8.889048775872605708249140016201753255599E3L
,
640 -2.465888106627948210478692168261494857089E4L
,
641 -3.934642211710774494879042116768390014289E4L
,
642 -3.455077258242252974937480623730228841003E4L
,
643 -1.524083977439690284820586063729912653196E4L
,
644 -2.810541887397984804237552337349093953857E3L
,
645 -1.343929553541159933824901621702567066156E2L
649 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
651 Peak relative error 8.4e-37 */
653 static const long double RNr3
[NRNr3
+ 1] =
655 -1.952401126551202208698629992497306292987E-6L,
656 -2.130881743066372952515162564941682716125E-4L,
657 -8.376493958090190943737529486107282224387E-3L,
658 -1.650592646560987700661598877522831234791E-1L,
659 -1.839290818933317338111364667708678163199E0L
,
660 -1.216278715570882422410442318517814388470E1L
,
661 -4.818759344462360427612133632533779091386E1L
,
662 -1.120994661297476876804405329172164436784E2L
,
663 -1.452850765662319264191141091859300126931E2L
,
664 -9.485207851128957108648038238656777241333E1L
,
665 -2.563663855025796641216191848818620020073E1L
,
666 -1.787995944187565676837847610706317833247E0L
669 static const long double RDr3
[NRDr3
+ 1] =
671 1.979130686770349481460559711878399476903E-4L,
672 1.156941716128488266238105813374635099057E-2L,
673 2.752657634309886336431266395637285974292E-1L,
674 3.482245457248318787349778336603569327521E0L
,
675 2.569347069372696358578399521203959253162E1L
,
676 1.142279000180457419740314694631879921561E2L
,
677 3.056503977190564294341422623108332700840E2L
,
678 4.780844020923794821656358157128719184422E2L
,
679 4.105972727212554277496256802312730410518E2L
,
680 1.724072188063746970865027817017067646246E2L
,
681 2.815939183464818198705278118326590370435E1L
685 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
687 Peak relative error 1.5e-36 */
689 static const long double RNr2
[NRNr2
+ 1] =
691 -2.638914383420287212401687401284326363787E-8L,
692 -3.479198370260633977258201271399116766619E-6L,
693 -1.783985295335697686382487087502222519983E-4L,
694 -4.777876933122576014266349277217559356276E-3L,
695 -7.450634738987325004070761301045014986520E-2L,
696 -7.068318854874733315971973707247467326619E-1L,
697 -4.113919921935944795764071670806867038732E0L
,
698 -1.440447573226906222417767283691888875082E1L
,
699 -2.883484031530718428417168042141288943905E1L
,
700 -2.990886974328476387277797361464279931446E1L
,
701 -1.325283914915104866248279787536128997331E1L
,
702 -1.572436106228070195510230310658206154374E0L
705 static const long double RDr2
[NRDr2
+ 1] =
707 2.675042728136731923554119302571867799673E-6L,
708 2.170997868451812708585443282998329996268E-4L,
709 7.249969752687540289422684951196241427445E-3L,
710 1.302040375859768674620410563307838448508E-1L,
711 1.380202483082910888897654537144485285549E0L
,
712 8.926594113174165352623847870299170069350E0L
,
713 3.521089584782616472372909095331572607185E1L
,
714 8.233547427533181375185259050330809105570E1L
,
715 1.072971579885803033079469639073292840135E2L
,
716 6.943803113337964469736022094105143158033E1L
,
717 1.775695341031607738233608307835017282662E1L
721 /* erfc(1/x) = 1/x exp(-1/x^2 - 0.5625 + R(1/x^2))
723 Peak relative error 2.2e-36 */
725 static const long double RNr1
[NRNr1
+ 1] =
727 -4.250780883202361946697751475473042685782E-8L,
728 -5.375777053288612282487696975623206383019E-6L,
729 -2.573645949220896816208565944117382460452E-4L,
730 -6.199032928113542080263152610799113086319E-3L,
731 -8.262721198693404060380104048479916247786E-2L,
732 -6.242615227257324746371284637695778043982E-1L,
733 -2.609874739199595400225113299437099626386E0L
,
734 -5.581967563336676737146358534602770006970E0L
,
735 -5.124398923356022609707490956634280573882E0L
,
736 -1.290865243944292370661544030414667556649E0L
739 static const long double RDr1
[NRDr1
+ 1] =
741 4.308976661749509034845251315983612976224E-6L,
742 3.265390126432780184125233455960049294580E-4L,
743 9.811328839187040701901866531796570418691E-3L,
744 1.511222515036021033410078631914783519649E-1L,
745 1.289264341917429958858379585970225092274E0L
,
746 6.147640356182230769548007536914983522270E0L
,
747 1.573966871337739784518246317003956180750E1L
,
748 1.955534123435095067199574045529218238263E1L
,
749 9.472613121363135472247929109615785855865E0L
755 __erfl (long double x
)
759 ieee854_long_double_shape_type u
;
763 ix
= sign
& 0x7fffffff;
765 if (ix
>= 0x7ff00000)
767 i
= ((sign
& 0xfff00000) >> 31) << 1;
768 return (long double) (1 - i
) + one
/ x
; /* erf(+-inf)=+-1 */
771 if (ix
>= 0x3ff00000) /* |x| >= 1.0 */
775 /* return (one - __erfcl (x)); */
780 if (ix
< 0x3fec0000) /* a < 0.875 */
782 if (ix
< 0x3c600000) /* |x|<2**-57 */
786 /* erf (-0) = -0. Unfortunately, for IBM extended double
787 0.125 * (8.0 * x + efx8 * x) for x = -0 evaluates to 0. */
790 return 0.125 * (8.0 * x
+ efx8
* x
); /*avoid underflow */
794 y
= a
+ a
* neval (z
, TN1
, NTN1
) / deval (z
, TD1
, NTD1
);
799 y
= erf_const
+ neval (a
, TN2
, NTN2
) / deval (a
, TD2
, NTD2
);
802 if (sign
& 0x80000000) /* x < 0 */
807 long_double_symbol (libm
, __erfl
, erfl
);
809 __erfcl (long double x
)
811 long double y
, z
, p
, r
;
813 ieee854_long_double_shape_type u
;
817 ix
= sign
& 0x7fffffff;
820 if (ix
>= 0x7ff00000)
821 { /* erfc(nan)=nan */
822 /* erfc(+-inf)=0,2 */
823 return (long double) (((u_int32_t
) sign
>> 31) << 1) + one
/ x
;
826 if (ix
< 0x3fd00000) /* |x| <1/4 */
828 if (ix
< 0x38d00000) /* |x|<2**-114 */
830 return one
- __erfl (x
);
832 if (ix
< 0x3ff40000) /* 1.25 */
840 y
= C13b
+ z
* neval (z
, RNr13
, NRNr13
) / deval (z
, RDr13
, NRDr13
);
845 y
= C14b
+ z
* neval (z
, RNr14
, NRNr14
) / deval (z
, RDr14
, NRDr14
);
850 y
= C15b
+ z
* neval (z
, RNr15
, NRNr15
) / deval (z
, RDr15
, NRDr15
);
855 y
= C16b
+ z
* neval (z
, RNr16
, NRNr16
) / deval (z
, RDr16
, NRDr16
);
860 y
= C17b
+ z
* neval (z
, RNr17
, NRNr17
) / deval (z
, RDr17
, NRDr17
);
865 y
= C18b
+ z
* neval (z
, RNr18
, NRNr18
) / deval (z
, RDr18
, NRDr18
);
870 y
= C19b
+ z
* neval (z
, RNr19
, NRNr19
) / deval (z
, RDr19
, NRDr19
);
875 y
= C20b
+ z
* neval (z
, RNr20
, NRNr20
) / deval (z
, RDr20
, NRDr20
);
879 if (sign
& 0x80000000)
883 /* 1.25 < |x| < 107 */
887 if ((ix
>= 0x40220000) && (sign
& 0x80000000))
897 p
= neval (z
, RNr1
, NRNr1
) / deval (z
, RDr1
, NRDr1
);
900 p
= neval (z
, RNr2
, NRNr2
) / deval (z
, RDr2
, NRDr2
);
903 p
= neval (z
, RNr3
, NRNr3
) / deval (z
, RDr3
, NRDr3
);
906 p
= neval (z
, RNr4
, NRNr4
) / deval (z
, RDr4
, NRDr4
);
909 p
= neval (z
, RNr5
, NRNr5
) / deval (z
, RDr5
, NRDr5
);
912 p
= neval (z
, RNr6
, NRNr6
) / deval (z
, RDr6
, NRDr6
);
915 p
= neval (z
, RNr7
, NRNr7
) / deval (z
, RDr7
, NRDr7
);
918 p
= neval (z
, RNr8
, NRNr8
) / deval (z
, RDr8
, NRDr8
);
923 u
.parts32
.w2
&= 0xffffe000;
925 r
= __ieee754_expl (-z
* z
- 0.5625) *
926 __ieee754_expl ((z
- x
) * (z
+ x
) + p
);
927 if ((sign
& 0x80000000) == 0)
934 if ((sign
& 0x80000000) == 0)
941 long_double_symbol (libm
, __erfcl
, erfcl
);