(__expm1l, expm1l): Remove NO_LONG_DOUBLE aliases.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-96 / s_erfl.c
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 /* Long double expansions contributed by
13    Stephen L. Moshier <moshier@na-net.ornl.gov>  */
14
15 /* double erf(double x)
16  * double erfc(double x)
17  *                           x
18  *                    2      |\
19  *     erf(x)  =  ---------  | exp(-t*t)dt
20  *                 sqrt(pi) \|
21  *                           0
22  *
23  *     erfc(x) =  1-erf(x)
24  *  Note that
25  *              erf(-x) = -erf(x)
26  *              erfc(-x) = 2 - erfc(x)
27  *
28  * Method:
29  *      1. For |x| in [0, 0.84375]
30  *          erf(x)  = x + x*R(x^2)
31  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
32  *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
33  *         Remark. The formula is derived by noting
34  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
35  *         and that
36  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
37  *         is close to one. The interval is chosen because the fix
38  *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
39  *         near 0.6174), and by some experiment, 0.84375 is chosen to
40  *         guarantee the error is less than one ulp for erf.
41  *
42  *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
43  *         c = 0.84506291151 rounded to single (24 bits)
44  *      erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
45  *      erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
46  *                        1+(c+P1(s)/Q1(s))    if x < 0
47  *         Remark: here we use the taylor series expansion at x=1.
48  *              erf(1+s) = erf(1) + s*Poly(s)
49  *                       = 0.845.. + P1(s)/Q1(s)
50  *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
51  *
52  *      3. For x in [1.25,1/0.35(~2.857143)],
53  *      erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
54  *              z=1/x^2
55  *      erf(x)  = 1 - erfc(x)
56  *
57  *      4. For x in [1/0.35,107]
58  *      erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
59  *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
60  *                             if -6.666<x<0
61  *                      = 2.0 - tiny            (if x <= -6.666)
62  *              z=1/x^2
63  *      erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
64  *      erf(x)  = sign(x)*(1.0 - tiny)
65  *      Note1:
66  *         To compute exp(-x*x-0.5625+R/S), let s be a single
67  *         precision number and s := x; then
68  *              -x*x = -s*s + (s-x)*(s+x)
69  *              exp(-x*x-0.5626+R/S) =
70  *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
71  *      Note2:
72  *         Here 4 and 5 make use of the asymptotic series
73  *                        exp(-x*x)
74  *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
75  *                        x*sqrt(pi)
76  *
77  *      5. For inf > x >= 107
78  *      erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
79  *      erfc(x) = tiny*tiny (raise underflow) if x > 0
80  *                      = 2 - tiny if x<0
81  *
82  *      7. Special case:
83  *      erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
84  *      erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
85  *              erfc/erf(NaN) is NaN
86  */
87
88
89 #include "math.h"
90 #include "math_private.h"
91
92 #ifdef __STDC__
93 static const long double
94 #else
95 static long double
96 #endif
97 tiny = 1e-4931L,
98   half = 0.5L,
99   one = 1.0L,
100   two = 2.0L,
101         /* c = (float)0.84506291151 */
102   erx = 0.845062911510467529296875L,
103 /*
104  * Coefficients for approximation to  erf on [0,0.84375]
105  */
106   /* 2/sqrt(pi) - 1 */
107   efx = 1.2837916709551257389615890312154517168810E-1L,
108   /* 8 * (2/sqrt(pi) - 1) */
109   efx8 = 1.0270333367641005911692712249723613735048E0L,
110
111   pp[6] = {
112     1.122751350964552113068262337278335028553E6L,
113     -2.808533301997696164408397079650699163276E6L,
114     -3.314325479115357458197119660818768924100E5L,
115     -6.848684465326256109712135497895525446398E4L,
116     -2.657817695110739185591505062971929859314E3L,
117     -1.655310302737837556654146291646499062882E2L,
118   },
119
120   qq[6] = {
121     8.745588372054466262548908189000448124232E6L,
122     3.746038264792471129367533128637019611485E6L,
123     7.066358783162407559861156173539693900031E5L,
124     7.448928604824620999413120955705448117056E4L,
125     4.511583986730994111992253980546131408924E3L,
126     1.368902937933296323345610240009071254014E2L,
127     /* 1.000000000000000000000000000000000000000E0 */
128   },
129
130 /*
131  * Coefficients for approximation to  erf  in [0.84375,1.25]
132  */
133 /* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
134    -0.15625 <= x <= +.25
135    Peak relative error 8.5e-22  */
136
137   pa[8] = {
138     -1.076952146179812072156734957705102256059E0L,
139      1.884814957770385593365179835059971587220E2L,
140     -5.339153975012804282890066622962070115606E1L,
141      4.435910679869176625928504532109635632618E1L,
142      1.683219516032328828278557309642929135179E1L,
143     -2.360236618396952560064259585299045804293E0L,
144      1.852230047861891953244413872297940938041E0L,
145      9.394994446747752308256773044667843200719E-2L,
146   },
147
148   qa[7] =  {
149     4.559263722294508998149925774781887811255E2L,
150     3.289248982200800575749795055149780689738E2L,
151     2.846070965875643009598627918383314457912E2L,
152     1.398715859064535039433275722017479994465E2L,
153     6.060190733759793706299079050985358190726E1L,
154     2.078695677795422351040502569964299664233E1L,
155     4.641271134150895940966798357442234498546E0L,
156     /* 1.000000000000000000000000000000000000000E0 */
157   },
158
159 /*
160  * Coefficients for approximation to  erfc in [1.25,1/0.35]
161  */
162 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
163    1/2.85711669921875 < 1/x < 1/1.25
164    Peak relative error 3.1e-21  */
165
166     ra[] = {
167       1.363566591833846324191000679620738857234E-1L,
168       1.018203167219873573808450274314658434507E1L,
169       1.862359362334248675526472871224778045594E2L,
170       1.411622588180721285284945138667933330348E3L,
171       5.088538459741511988784440103218342840478E3L,
172       8.928251553922176506858267311750789273656E3L,
173       7.264436000148052545243018622742770549982E3L,
174       2.387492459664548651671894725748959751119E3L,
175       2.220916652813908085449221282808458466556E2L,
176     },
177
178     sa[] = {
179       -1.382234625202480685182526402169222331847E1L,
180       -3.315638835627950255832519203687435946482E2L,
181       -2.949124863912936259747237164260785326692E3L,
182       -1.246622099070875940506391433635999693661E4L,
183       -2.673079795851665428695842853070996219632E4L,
184       -2.880269786660559337358397106518918220991E4L,
185       -1.450600228493968044773354186390390823713E4L,
186       -2.874539731125893533960680525192064277816E3L,
187       -1.402241261419067750237395034116942296027E2L,
188       /* 1.000000000000000000000000000000000000000E0 */
189     },
190 /*
191  * Coefficients for approximation to  erfc in [1/.35,107]
192  */
193 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
194    1/6.6666259765625 < 1/x < 1/2.85711669921875
195    Peak relative error 4.2e-22  */
196     rb[] = {
197       -4.869587348270494309550558460786501252369E-5L,
198       -4.030199390527997378549161722412466959403E-3L,
199       -9.434425866377037610206443566288917589122E-2L,
200       -9.319032754357658601200655161585539404155E-1L,
201       -4.273788174307459947350256581445442062291E0L,
202       -8.842289940696150508373541814064198259278E0L,
203       -7.069215249419887403187988144752613025255E0L,
204       -1.401228723639514787920274427443330704764E0L,
205     },
206
207     sb[] = {
208       4.936254964107175160157544545879293019085E-3L,
209       1.583457624037795744377163924895349412015E-1L,
210       1.850647991850328356622940552450636420484E0L,
211       9.927611557279019463768050710008450625415E0L,
212       2.531667257649436709617165336779212114570E1L,
213       2.869752886406743386458304052862814690045E1L,
214       1.182059497870819562441683560749192539345E1L,
215       /* 1.000000000000000000000000000000000000000E0 */
216     },
217 /* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
218    1/107 <= 1/x <= 1/6.6666259765625
219    Peak relative error 1.1e-21  */
220     rc[] = {
221       -8.299617545269701963973537248996670806850E-5L,
222       -6.243845685115818513578933902532056244108E-3L,
223       -1.141667210620380223113693474478394397230E-1L,
224       -7.521343797212024245375240432734425789409E-1L,
225       -1.765321928311155824664963633786967602934E0L,
226       -1.029403473103215800456761180695263439188E0L,
227     },
228
229     sc[] = {
230       8.413244363014929493035952542677768808601E-3L,
231       2.065114333816877479753334599639158060979E-1L,
232       1.639064941530797583766364412782135680148E0L,
233       4.936788463787115555582319302981666347450E0L,
234       5.005177727208955487404729933261347679090E0L,
235       /* 1.000000000000000000000000000000000000000E0 */
236     };
237
238 #ifdef __STDC__
239 long double
240 __erfl (long double x)
241 #else
242 long double
243 __erfl (x)
244      long double x;
245 #endif
246 {
247   long double R, S, P, Q, s, y, z, r;
248   int32_t ix, i;
249   u_int32_t se, i0, i1;
250
251   GET_LDOUBLE_WORDS (se, i0, i1, x);
252   ix = se & 0x7fff;
253
254   if (ix >= 0x7fff)
255     {                           /* erf(nan)=nan */
256       i = ((se & 0xffff) >> 15) << 1;
257       return (long double) (1 - i) + one / x;   /* erf(+-inf)=+-1 */
258     }
259
260   ix = (ix << 16) | (i0 >> 16);
261   if (ix < 0x3ffed800) /* |x|<0.84375 */
262     {
263       if (ix < 0x3fde8000) /* |x|<2**-33 */
264         {
265           if (ix < 0x00080000)
266             return 0.125 * (8.0 * x + efx8 * x);        /*avoid underflow */
267           return x + efx * x;
268         }
269       z = x * x;
270       r = pp[0] + z * (pp[1]
271           + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
272       s = qq[0] + z * (qq[1]
273           + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
274       y = r / s;
275       return x + x * y;
276     }
277   if (ix < 0x3fffa000) /* 1.25 */
278     {                           /* 0.84375 <= |x| < 1.25 */
279       s = fabsl (x) - one;
280       P = pa[0] + s * (pa[1] + s * (pa[2]
281         + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
282       Q = qa[0] + s * (qa[1] + s * (qa[2]
283         + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
284       if ((se & 0x8000) == 0)
285         return erx + P / Q;
286       else
287         return -erx - P / Q;
288     }
289   if (ix >= 0x4001d555) /* 6.6666259765625 */
290     {                           /* inf>|x|>=6.666 */
291       if ((se & 0x8000) == 0)
292         return one - tiny;
293       else
294         return tiny - one;
295     }
296   x = fabsl (x);
297   s = one / (x * x);
298   if (ix < 0x4000b6db) /* 2.85711669921875 */
299     {
300       R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
301           s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
302       S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
303           s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
304     }
305   else
306     {                           /* |x| >= 1/0.35 */
307       R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
308          s * (rb[5] + s * (rb[6] + s * rb[7]))))));
309       S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
310          s * (sb[5] + s * (sb[6] + s))))));
311     }
312   z = x;
313   GET_LDOUBLE_WORDS (i, i0, i1, z);
314   i1 = 0;
315   SET_LDOUBLE_WORDS (z, i, i0, i1);
316   r =
317     __ieee754_expl (-z * z - 0.5625) * __ieee754_expl ((z - x) * (z + x) +
318                                                      R / S);
319   if ((se & 0x8000) == 0)
320     return one - r / x;
321   else
322     return r / x - one;
323 }
324
325 weak_alias (__erfl, erfl)
326 #ifdef NO_LONG_DOUBLE
327 strong_alias (__erf, __erfl)
328 weak_alias (__erf, erfl)
329 #endif
330 #ifdef __STDC__
331      long double
332      __erfcl (long double x)
333 #else
334      long double
335      __erfcl (x)
336      long double x;
337 #endif
338 {
339   int32_t hx, ix;
340   long double R, S, P, Q, s, y, z, r;
341   u_int32_t se, i0, i1;
342
343   GET_LDOUBLE_WORDS (se, i0, i1, x);
344   ix = se & 0x7fff;
345   if (ix >= 0x7fff)
346     {                           /* erfc(nan)=nan */
347       /* erfc(+-inf)=0,2 */
348       return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
349     }
350
351   ix = (ix << 16) | (i0 >> 16);
352   if (ix < 0x3ffed800) /* |x|<0.84375 */
353     {
354       if (ix < 0x3fbe0000) /* |x|<2**-65 */
355         return one - x;
356       z = x * x;
357       r = pp[0] + z * (pp[1]
358           + z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
359       s = qq[0] + z * (qq[1]
360           + z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
361       y = r / s;
362       if (ix < 0x3ffd8000) /* x<1/4 */
363         {
364           return one - (x + x * y);
365         }
366       else
367         {
368           r = x * y;
369           r += (x - half);
370           return half - r;
371         }
372     }
373   if (ix < 0x3fffa000) /* 1.25 */
374     {                           /* 0.84375 <= |x| < 1.25 */
375       s = fabsl (x) - one;
376       P = pa[0] + s * (pa[1] + s * (pa[2]
377         + s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
378       Q = qa[0] + s * (qa[1] + s * (qa[2]
379         + s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
380       if ((se & 0x8000) == 0)
381         {
382           z = one - erx;
383           return z - P / Q;
384         }
385       else
386         {
387           z = erx + P / Q;
388           return one + z;
389         }
390     }
391   if (ix < 0x4005d600) /* 107 */
392     {                           /* |x|<107 */
393       x = fabsl (x);
394       s = one / (x * x);
395       if (ix < 0x4000b6db) /* 2.85711669921875 */
396         {                       /* |x| < 1/.35 ~ 2.857143 */
397           R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
398               s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
399           S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
400               s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
401         }
402       else if (ix < 0x4001d555) /* 6.6666259765625 */
403         {                       /* 6.666 > |x| >= 1/.35 ~ 2.857143 */
404           R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
405               s * (rb[5] + s * (rb[6] + s * rb[7]))))));
406           S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
407               s * (sb[5] + s * (sb[6] + s))))));
408         }
409       else
410         {                       /* |x| >= 6.666 */
411           if (se & 0x8000)
412             return two - tiny;  /* x < -6.666 */
413
414           R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
415                                                     s * (rc[4] + s * rc[5]))));
416           S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
417                                                     s * (sc[4] + s))));
418         }
419       z = x;
420       GET_LDOUBLE_WORDS (hx, i0, i1, z);
421       i1 = 0;
422       i0 &= 0xffffff00;
423       SET_LDOUBLE_WORDS (z, hx, i0, i1);
424       r = __ieee754_expl (-z * z - 0.5625) *
425         __ieee754_expl ((z - x) * (z + x) + R / S);
426       if ((se & 0x8000) == 0)
427         return r / x;
428       else
429         return two - r / x;
430     }
431   else
432     {
433       if ((se & 0x8000) == 0)
434         return tiny * tiny;
435       else
436         return two - tiny;
437     }
438 }
439
440 weak_alias (__erfcl, erfcl)
441 #ifdef NO_LONG_DOUBLE
442 strong_alias (__erfc, __erfcl)
443 weak_alias (__erfc, erfcl)
444 #endif