```   1 /*
3  *
4  * The contents of this file are subject to the terms of the
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
18  *
20  */
21
22 /*
24  */
25 /*
27  * Use is subject to license terms.
28  */
29
30 /*
31  * expl(x)
32  * Table driven method
33  * Written by K.C. Ng, November 1988.
34  * Algorithm :
35  *      1. Argument Reduction: given the input x, find r and integer k
36  *         and j such that
37  *                   x = (32k+j)*ln2 + r,  |r| <= (1/64)*ln2 .
38  *
39  *      2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
40  *         Note:
41  *         a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
42  *         b. 2^(j/32) is represented as
43  *                      _TBL_expl_hi[j]+_TBL_expl_lo[j]
44  *         where
45  *              _TBL_expl_hi[j] = 2^(j/32) rounded
46  *              _TBL_expl_lo[j] = 2^(j/32) - _TBL_expl_hi[j].
47  *
48  * Special cases:
49  *      expl(INF) is INF, expl(NaN) is NaN;
50  *      expl(-INF)=  0;
51  *      for finite argument, only expl(0)=1 is exact.
52  *
53  * Accuracy:
54  *      according to an error analysis, the error is always less than
55  *      an ulp (unit in the last place).
56  *
57  * Misc. info.
58  *      For 113 bit long double
59  *              if x >  1.135652340629414394949193107797076342845e+4
60  *      then expl(x) overflow;
61  *              if x < -1.143346274333629787883724384345262150341e+4
62  *      then expl(x) underflow
63  *
64  * Constants:
65  * Only decimal values are given. We assume that the compiler will convert
66  * from decimal to binary accurately enough to produce the correct
68  */
69
70 #pragma weak __expl = expl
71
72 #include "libm.h"
73
74 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
75
76 static const long double
77 one             =  1.0L,
78 two             =  2.0L,
79 ln2_64          =  1.083042469624914545964425189778400898568e-2L,
80 ovflthreshold   =  1.135652340629414394949193107797076342845e+4L,
81 unflthreshold   = -1.143346274333629787883724384345262150341e+4L,
82 invln2_32       =  4.616624130844682903551758979206054839765e+1L,
83 ln2_32hi        =  2.166084939249829091928849858592451515688e-2L,
84 ln2_32lo        =  5.209643502595475652782654157501186731779e-27L;
85
86 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
87 static const long double
88 t1 =   1.666666666666666666666666666660876387437e-1L,
89 t2 =  -2.777777777777777777777707812093173478756e-3L,
90 t3 =   6.613756613756613482074280932874221202424e-5L,
91 t4 =  -1.653439153392139954169609822742235851120e-6L,
92 t5 =   4.175314851769539751387852116610973796053e-8L;
93
94 long double
95 expl(long double x) {
96         int *px = (int *) &x, ix, j, k, m;
97         long double t, r;
98
99         ix = px[0];                             /* high word of x */
100         if (ix >= 0x7fff0000)
101                 return (x + x);                 /* NaN of +inf */
102         if (((unsigned) ix) >= 0xffff0000)
103                 return (-one / x);              /* NaN or -inf */
104         if ((ix & 0x7fffffff) < 0x3fc30000) {
105                 if ((int) x < 1)
106                         return (one + x);       /* |x|<2^-60 */
107         }
108         if (ix > 0) {
109                 if (x > ovflthreshold)
110                         return (scalbnl(x, 20000));
111                 k = (int) (invln2_32 * (x + ln2_64));
112         } else {
113                 if (x < unflthreshold)
114                         return (scalbnl(-x, -40000));
115                 k = (int) (invln2_32 * (x - ln2_64));
116         }
117         j  = k&0x1f;
118         m  = k>>5;
119         t  = (long double) k;
120         x  = (x - t * ln2_32hi) - t * ln2_32lo;
121         t  = x * x;
122         r  = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
123         x  = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
124                 _TBL_expl_lo[j]);
125         return (scalbnl(x, m));
126 }
```