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 #pragma weak __expm1l = expm1l
31
32 #if !defined(__sparc)
33 #error Unsupported architecture
34 #endif
35
36 /*
37 * expm1l(x)
38 *
39 * Table driven method
40 * Written by K.C. Ng, June 1995.
41 * Algorithm :
42 * 1. expm1(x) = x if x<2**-114
43 * 2. if |x| <= 0.0625 = 1/16, use approximation
44 * expm1(x) = x + x*P/(2-P)
83 *
84 * Accuracy:
85 * according to an error analysis, the error is always less than
86 * 2 ulp (unit in the last place).
87 *
88 * Misc. info.
89 * For 113 bit long double
90 * if x > 1.135652340629414394949193107797076342845e+4
91 * then expm1l(x) overflow;
92 *
93 * Constants:
94 * Only decimal values are given. We assume that the compiler will convert
95 * from decimal to binary accurately enough to produce the correct
96 * hexadecimal values.
97 */
98
99 #include "libm.h"
100
101 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
102 extern const long double _TBL_expm1lx[], _TBL_expm1l[];
103
104 static const long double
105 zero = +0.0L,
106 one = +1.0L,
107 two = +2.0L,
108 ln2_64 = +1.083042469624914545964425189778400898568e-2L,
109 ovflthreshold = +1.135652340629414394949193107797076342845e+4L,
110 invln2_32 = +4.616624130844682903551758979206054839765e+1L,
111 ln2_32hi = +2.166084939249829091928849858592451515688e-2L,
112 ln2_32lo = +5.209643502595475652782654157501186731779e-27L,
113 huge = +1.0e4000L,
114 tiny = +1.0e-4000L,
115 P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
116 P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
117 P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
118 P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
119 P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
120 P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
121 P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
122 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
123 T1 = +1.666666666666666666666666666660876387437e-1L,
124 T2 = -2.777777777777777777777707812093173478756e-3L,
125 T3 = +6.613756613756613482074280932874221202424e-5L,
126 T4 = -1.653439153392139954169609822742235851120e-6L,
127 T5 = +4.175314851769539751387852116610973796053e-8L;
128
129 long double
130 expm1l(long double x) {
131 int hx, ix, j, k, m;
132 long double t, r, s, w;
133
134 hx = ((int *) &x)[HIXWORD];
135 ix = hx & ~0x80000000;
136 if (ix >= 0x7fff0000) {
137 if (x != x)
138 return (x + x); /* NaN */
139 if (x < zero)
140 return (-one); /* -inf */
141 return (x); /* +inf */
142 }
143 if (ix < 0x3fff4000) { /* |x| < 1.25 */
144 if (ix < 0x3ffb0000) { /* |x| < 0.0625 */
145 if (ix < 0x3f8d0000) {
146 if ((int) x == 0)
147 return (x); /* |x|<2^-114 */
148 }
149 t = x * x;
150 r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
151 (P5 + t * (P6 + t * P7)))))));
152 return (x + (x * r) / (two - r));
153 }
154 /* compute i = [64*x] */
155 m = 0x4009 - (ix >> 16);
156 j = ((ix & 0x0000ffff) | 0x10000) >> m; /* j=4,...,67 */
157 if (hx < 0)
158 j += 82; /* negative */
159 s = x - _TBL_expm1lx[j];
160 t = s * s;
161 r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
162 r = (s + s) / (two - r);
163 w = _TBL_expm1l[j];
164 return (w + (w + one) * r);
165 }
166 if (hx > 0) {
167 if (x > ovflthreshold)
168 return (huge * huge);
169 k = (int) (invln2_32 * (x + ln2_64));
170 } else {
171 if (x < -80.0)
172 return (tiny - x / x);
173 k = (int) (invln2_32 * (x - ln2_64));
174 }
175 j = k & 0x1f;
176 m = k >> 5;
177 t = (long double) k;
178 x = (x - t * ln2_32hi) - t * ln2_32lo;
179 t = x * x;
180 r = (x - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))))) - two;
181 x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
182 _TBL_expl_lo[j]);
183 return (scalbnl(x, m) - one);
184 }
|
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 #pragma weak __expm1l = expm1l
32
33 #if !defined(__sparc)
34 #error Unsupported architecture
35 #endif
36
37 /*
38 * expm1l(x)
39 *
40 * Table driven method
41 * Written by K.C. Ng, June 1995.
42 * Algorithm :
43 * 1. expm1(x) = x if x<2**-114
44 * 2. if |x| <= 0.0625 = 1/16, use approximation
45 * expm1(x) = x + x*P/(2-P)
84 *
85 * Accuracy:
86 * according to an error analysis, the error is always less than
87 * 2 ulp (unit in the last place).
88 *
89 * Misc. info.
90 * For 113 bit long double
91 * if x > 1.135652340629414394949193107797076342845e+4
92 * then expm1l(x) overflow;
93 *
94 * Constants:
95 * Only decimal values are given. We assume that the compiler will convert
96 * from decimal to binary accurately enough to produce the correct
97 * hexadecimal values.
98 */
99
100 #include "libm.h"
101
102 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
103 extern const long double _TBL_expm1lx[], _TBL_expm1l[];
104 static const long double zero = +0.0L,
105 one = +1.0L,
106 two = +2.0L,
107 ln2_64 = +1.083042469624914545964425189778400898568e-2L,
108 ovflthreshold = +1.135652340629414394949193107797076342845e+4L,
109 invln2_32 = +4.616624130844682903551758979206054839765e+1L,
110 ln2_32hi = +2.166084939249829091928849858592451515688e-2L,
111 ln2_32lo = +5.209643502595475652782654157501186731779e-27L,
112 huge = +1.0e4000L,
113 tiny = +1.0e-4000L,
114 P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
115 P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
116 P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
117 P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
118 P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
119 P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
120 P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
121 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
122 T1 = +1.666666666666666666666666666660876387437e-1L,
123 T2 = -2.777777777777777777777707812093173478756e-3L,
124 T3 = +6.613756613756613482074280932874221202424e-5L,
125 T4 = -1.653439153392139954169609822742235851120e-6L,
126 T5 = +4.175314851769539751387852116610973796053e-8L;
127
128 long double
129 expm1l(long double x)
130 {
131 int hx, ix, j, k, m;
132 long double t, r, s, w;
133
134 hx = ((int *)&x)[HIXWORD];
135 ix = hx & ~0x80000000;
136
137 if (ix >= 0x7fff0000) {
138 if (x != x)
139 return (x + x); /* NaN */
140
141 if (x < zero)
142 return (-one); /* -inf */
143
144 return (x); /* +inf */
145 }
146
147 if (ix < 0x3fff4000) { /* |x| < 1.25 */
148 if (ix < 0x3ffb0000) { /* |x| < 0.0625 */
149 if (ix < 0x3f8d0000) {
150 if ((int)x == 0)
151 return (x); /* |x|<2^-114 */
152 }
153
154 t = x * x;
155 r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
156 (P5 + t * (P6 + t * P7)))))));
157 return (x + (x * r) / (two - r));
158 }
159
160 /* compute i = [64*x] */
161 m = 0x4009 - (ix >> 16);
162 j = ((ix & 0x0000ffff) | 0x10000) >> m; /* j=4,...,67 */
163
164 if (hx < 0)
165 j += 82; /* negative */
166
167 s = x - _TBL_expm1lx[j];
168 t = s * s;
169 r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
170 r = (s + s) / (two - r);
171 w = _TBL_expm1l[j];
172 return (w + (w + one) * r);
173 }
174
175 if (hx > 0) {
176 if (x > ovflthreshold)
177 return (huge * huge);
178
179 k = (int)(invln2_32 * (x + ln2_64));
180 } else {
181 if (x < -80.0)
182 return (tiny - x / x);
183
184 k = (int)(invln2_32 * (x - ln2_64));
185 }
186
187 j = k & 0x1f;
188 m = k >> 5;
189 t = (long double)k;
190 x = (x - t * ln2_32hi) - t * ln2_32lo;
191 t = x * x;
192 r = (x - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))))) - two;
193 x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
194 _TBL_expl_lo[j]);
195 return (scalbnl(x, m) - one);
196 }
|