1/**
2 * This file has no copyright assigned and is placed in the Public Domain.
3 * This file is part of the mingw-w64 runtime package.
4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
5 */
6/* erfl.c
7 *
8 * Error function
9 *
10 *
11 *
12 * SYNOPSIS:
13 *
14 * long double x, y, erfl();
15 *
16 * y = erfl( x );
17 *
18 *
19 *
20 * DESCRIPTION:
21 *
22 * The integral is
23 *
24 * x
25 * -
26 * 2 | | 2
27 * erf(x) = -------- | exp( - t ) dt.
28 * sqrt(pi) | |
29 * -
30 * 0
31 *
32 * The magnitude of x is limited to about 106.56 for IEEE
33 * arithmetic; 1 or -1 is returned outside this range.
34 *
35 * For 0 <= |x| < 1, erf(x) = x * P6(x^2)/Q6(x^2);
36 * Otherwise: erf(x) = 1 - erfc(x).
37 *
38 *
39 *
40 * ACCURACY:
41 *
42 * Relative error:
43 * arithmetic domain # trials peak rms
44 * IEEE 0,1 50000 2.0e-19 5.7e-20
45 *
46 */
47
48/* erfcl.c
49 *
50 * Complementary error function
51 *
52 *
53 *
54 * SYNOPSIS:
55 *
56 * long double x, y, erfcl();
57 *
58 * y = erfcl( x );
59 *
60 *
61 *
62 * DESCRIPTION:
63 *
64 *
65 * 1 - erf(x) =
66 *
67 * inf.
68 * -
69 * 2 | | 2
70 * erfc(x) = -------- | exp( - t ) dt
71 * sqrt(pi) | |
72 * -
73 * x
74 *
75 *
76 * For small x, erfc(x) = 1 - erf(x); otherwise rational
77 * approximations are computed.
78 *
79 * A special function expx2l.c is used to suppress error amplification
80 * in computing exp(-x^2).
81 *
82 *
83 * ACCURACY:
84 *
85 * Relative error:
86 * arithmetic domain # trials peak rms
87 * IEEE 0,13 50000 8.4e-19 9.7e-20
88 * IEEE 6,106.56 20000 2.9e-19 7.1e-20
89 *
90 *
91 * ERROR MESSAGES:
92 *
93 * message condition value returned
94 * erfcl underflow x^2 > MAXLOGL 0.0
95 *
96 *
97 */
98
99
100/*
101Modified from file ndtrl.c
102Cephes Math Library Release 2.3: January, 1995
103Copyright 1984, 1995 by Stephen L. Moshier
104*/
105
106#include <math.h>
107#include "cephes_mconf.h"
108
109long double erfl(long double x);
110
111#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__
112long double erfcl(long double x)
113{
114 return erfc(x);
115}
116
117long double erfl(long double x)
118{
119 return erf(x);
120}
121#else
122/* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
123 1/8 <= 1/x <= 1
124 Peak relative error 5.8e-21 */
125
126static const uLD P[10] = {
127 { { 0x4bf0,0x9ad8,0x7a03,0x86c7,0x401d, 0, 0, 0 } },
128 { { 0xdf23,0xd843,0x4032,0x8881,0x401e, 0, 0, 0 } },
129 { { 0xd025,0xcfd5,0x8494,0x88d3,0x401e, 0, 0, 0 } },
130 { { 0xb6d0,0xc92b,0x5417,0xacb1,0x401d, 0, 0, 0 } },
131 { { 0xada8,0x356a,0x4982,0x94a6,0x401c, 0, 0, 0 } },
132 { { 0x4e13,0xcaee,0x9e31,0xb258,0x401a, 0, 0, 0 } },
133 { { 0x5840,0x554d,0x37a3,0x9239,0x4018, 0, 0, 0 } },
134 { { 0x3b58,0x3da2,0xaf02,0x9780,0x4015, 0, 0, 0 } },
135 { { 0x0144,0x489e,0xbe68,0x9c31,0x4011, 0, 0, 0 } },
136 { { 0x333b,0xd9e6,0xd404,0x986f,0xbfee, 0, 0, 0 } }
137};
138static const uLD Q[] = {
139 { { 0x0e43,0x302d,0x79ed,0x86c7,0x401d, 0, 0, 0 } },
140 { { 0xf817,0x9128,0xc0f8,0xd48b,0x401e, 0, 0, 0 } },
141 { { 0x8eae,0x8dad,0x6eb4,0x9aa2,0x401f, 0, 0, 0 } },
142 { { 0x00e7,0x7595,0xcd06,0x88bb,0x401f, 0, 0, 0 } },
143 { { 0x4991,0xcfda,0x52f1,0xa2a9,0x401e, 0, 0, 0 } },
144 { { 0xc39d,0xe415,0xc43d,0x87c0,0x401d, 0, 0, 0 } },
145 { { 0xa75d,0x436f,0x30dd,0xa027,0x401b, 0, 0, 0 } },
146 { { 0xc4cb,0x305a,0xbf78,0x8220,0x4019, 0, 0, 0 } },
147 { { 0x3708,0x33b1,0x07fa,0x8644,0x4016, 0, 0, 0 } },
148 { { 0x24fa,0x96f6,0x7153,0x8a6c,0x4012, 0, 0, 0 } }
149};
150
151/* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
152 1/128 <= 1/x < 1/8
153 Peak relative error 1.9e-21 */
154
155static const uLD R[] = {
156 { { 0x260a,0xab95,0x2fc7,0xe7c4,0x4000, 0, 0, 0 } },
157 { { 0x4761,0x613e,0xdf6d,0xe58e,0x4001, 0, 0, 0 } },
158 { { 0x0615,0x4b00,0x575f,0xdc7b,0x4000, 0, 0, 0 } },
159 { { 0x521d,0x8527,0x3435,0x8dc2,0x3ffe, 0, 0, 0 } },
160 { { 0x22cf,0xc711,0x6c5b,0xdcfb,0x3ff9, 0, 0, 0 } }
161};
162static const uLD S[] = {
163 { { 0x5de6,0x17d7,0x54d6,0xaba9,0x4002, 0, 0, 0 } },
164 { { 0x55d5,0xd300,0xe71e,0xf564,0x4002, 0, 0, 0 } },
165 { { 0xb611,0x8f76,0xf020,0xd255,0x4001, 0, 0, 0 } },
166 { { 0x3684,0x3798,0xb793,0x80b0,0x3fff, 0, 0, 0 } },
167 { { 0xf5af,0x2fb2,0x1e57,0xc3d7,0x3ffa, 0, 0, 0 } }
168};
169
170/* erf(x) = x T(x^2)/U(x^2)
171 0 <= x <= 1
172 Peak relative error 7.6e-23 */
173
174static const uLD T[] = {
175 { { 0xfd7a,0x3a1a,0x705b,0xe0c4,0x3ffb, 0, 0, 0 } },
176 { { 0x3128,0xc337,0x3716,0xace5,0x4001, 0, 0, 0 } },
177 { { 0x9517,0x4e93,0x540e,0x8f97,0x4007, 0, 0, 0 } },
178 { { 0x6118,0x6059,0x9093,0xa757,0x400a, 0, 0, 0 } },
179 { { 0xb954,0xa987,0xc60c,0xbc83,0x400e, 0, 0, 0 } },
180 { { 0x7a56,0xe45a,0xa4bd,0x975b,0x4010, 0, 0, 0 } },
181 { { 0xc446,0x6bab,0x0b2a,0x86d0,0x4013, 0, 0, 0 } }
182};
183
184static const uLD U[] = {
185 { { 0x3453,0x1f8e,0xf688,0xb507,0x4004, 0, 0, 0 } },
186 { { 0x71ac,0xb12f,0x21ca,0xf2e2,0x4008, 0, 0, 0 } },
187 { { 0xffe8,0x9cac,0x3b84,0xc2ac,0x400c, 0, 0, 0 } },
188 { { 0x481d,0x445b,0xc807,0xc232,0x400f, 0, 0, 0 } },
189 { { 0x9ad5,0x1aef,0x45b1,0xe25e,0x4011, 0, 0, 0 } },
190 { { 0x71a7,0x1cad,0x012e,0xeef3,0x4012, 0, 0, 0 } }
191};
192
193/* expx2l.c
194 *
195 * Exponential of squared argument
196 *
197 *
198 *
199 * SYNOPSIS:
200 *
201 * long double x, y, expmx2l();
202 * int sign;
203 *
204 * y = expx2l( x );
205 *
206 *
207 *
208 * DESCRIPTION:
209 *
210 * Computes y = exp(x*x) while suppressing error amplification
211 * that would ordinarily arise from the inexactness of the
212 * exponential argument x*x.
213 *
214 *
215 *
216 * ACCURACY:
217 *
218 * Relative error:
219 * arithmetic domain # trials peak rms
220 * IEEE -106.566, 106.566 10^5 1.6e-19 4.4e-20
221 *
222 */
223
224#define M 32768.0L
225#define MINV 3.0517578125e-5L
226
227static long double expx2l (long double x)
228{
229 long double u, u1, m, f;
230
231 x = fabsl (x);
232 /* Represent x as an exact multiple of M plus a residual.
233 M is a power of 2 chosen so that exp(m * m) does not overflow
234 or underflow and so that |x - m| is small. */
235 m = MINV * floorl(M * x + 0.5L);
236 f = x - m;
237
238 /* x^2 = m^2 + 2mf + f^2 */
239 u = m * m;
240 u1 = 2 * m * f + f * f;
241
242 if ((u + u1) > MAXLOGL)
243 return (INFINITYL);
244
245 /* u is exact, u1 is small. */
246 u = expl(u) * expl(u1);
247 return (u);
248}
249
250long double erfcl(long double a)
251{
252 long double p, q, x, y, z;
253
254 if (isinf (a))
255 return (signbit(a) ? 2.0 : 0.0);
256
257 if (isnan (a))
258 return (a);
259
260 x = fabsl (a);
261
262 if (x < 1.0L)
263 return (1.0L - erfl(a));
264
265 z = a * a;
266
267 if (z > MAXLOGL)
268 {
269under:
270 mtherr("erfcl", UNDERFLOW);
271 errno = ERANGE;
272 return (signbit(a) ? 2.0 : 0.0);
273 }
274
275 /* Compute z = expl(a * a). */
276 z = expx2l(a);
277 y = 1.0L/x;
278
279 if (x < 8.0L)
280 {
281 p = polevll(y, P, 9);
282 q = p1evll(y, Q, 10);
283 }
284 else
285 {
286 q = y * y;
287 p = y * polevll(q, R, 4);
288 q = p1evll(q, S, 5);
289 }
290 y = p/(q * z);
291
292 if (a < 0.0L)
293 y = 2.0L - y;
294
295 if (y == 0.0L)
296 goto under;
297
298 return (y);
299}
300
301long double erfl(long double x)
302{
303 long double y, z;
304
305 if (x == 0.0L)
306 return (x);
307
308 if (isinf (x))
309 return (signbit(x) ? -1.0L : 1.0L);
310
311 if (fabsl(x) > 1.0L)
312 return (1.0L - erfcl(x));
313
314 z = x * x;
315 y = x * polevll(z, T, 6) / p1evll(z, U, 6);
316 return (y);
317}
318#endif