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
neval(__float128 x,const __float128 * p,int n)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
deval(__float128 x,const __float128 * p,int n)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
erfq(__float128 x)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
erfcq(__float128 x)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 }
949