5 * Common Development and Distribution License (the "License").
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
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
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
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
67 * hexadecimal values.
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 }
|
5 * Common Development and Distribution License (the "License").
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
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 /*
32 * expl(x)
33 * Table driven method
34 * Written by K.C. Ng, November 1988.
35 * Algorithm :
36 * 1. Argument Reduction: given the input x, find r and integer k
37 * and j such that
38 * x = (32k+j)*ln2 + r, |r| <= (1/64)*ln2 .
39 *
40 * 2. expl(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r))
41 * Note:
42 * a. expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
43 * b. 2^(j/32) is represented as
44 * _TBL_expl_hi[j]+_TBL_expl_lo[j]
45 * where
56 * an ulp (unit in the last place).
57 *
58 * Misc. info.
59 * For 113 bit long double
60 * if x > 1.135652340629414394949193107797076342845e+4
61 * then expl(x) overflow;
62 * if x < -1.143346274333629787883724384345262150341e+4
63 * then expl(x) underflow
64 *
65 * Constants:
66 * Only decimal values are given. We assume that the compiler will convert
67 * from decimal to binary accurately enough to produce the correct
68 * hexadecimal values.
69 */
70
71 #pragma weak __expl = expl
72
73 #include "libm.h"
74
75 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
76 static const long double one = 1.0L,
77 two = 2.0L,
78 ln2_64 = 1.083042469624914545964425189778400898568e-2L,
79 ovflthreshold = 1.135652340629414394949193107797076342845e+4L,
80 unflthreshold = -1.143346274333629787883724384345262150341e+4L,
81 invln2_32 = 4.616624130844682903551758979206054839765e+1L,
82 ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
83 ln2_32lo = 5.209643502595475652782654157501186731779e-27L;
84
85 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
86 static const long double t1 = 1.666666666666666666666666666660876387437e-1L,
87 t2 = -2.777777777777777777777707812093173478756e-3L,
88 t3 = 6.613756613756613482074280932874221202424e-5L,
89 t4 = -1.653439153392139954169609822742235851120e-6L,
90 t5 = 4.175314851769539751387852116610973796053e-8L;
91
92 long double
93 expl(long double x)
94 {
95 int *px = (int *)&x, ix, j, k, m;
96 long double t, r;
97
98 ix = px[0]; /* high word of x */
99
100 if (ix >= 0x7fff0000)
101 return (x + x); /* NaN of +inf */
102
103 if (((unsigned)ix) >= 0xffff0000)
104 return (-one / x); /* NaN or -inf */
105
106 if ((ix & 0x7fffffff) < 0x3fc30000) {
107 if ((int)x < 1)
108 return (one + x); /* |x|<2^-60 */
109 }
110
111 if (ix > 0) {
112 if (x > ovflthreshold)
113 return (scalbnl(x, 20000));
114
115 k = (int)(invln2_32 * (x + ln2_64));
116 } else {
117 if (x < unflthreshold)
118 return (scalbnl(-x, -40000));
119
120 k = (int)(invln2_32 * (x - ln2_64));
121 }
122
123 j = k & 0x1f;
124 m = k >> 5;
125 t = (long double)k;
126 x = (x - t * ln2_32hi) - t * ln2_32lo;
127 t = x * x;
128 r = (x - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) - two;
129 x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
130 _TBL_expl_lo[j]);
131 return (scalbnl(x, m));
132 }
|