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