Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/Q/expm1l.c
+++ new/usr/src/lib/libm/common/Q/expm1l.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
|
↓ open down ↓ |
14 lines elided |
↑ open up ↑ |
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 21
22 22 /*
23 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 +
25 26 /*
26 27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 28 * Use is subject to license terms.
28 29 */
29 30
30 31 #pragma weak __expm1l = expm1l
31 32
32 33 #if !defined(__sparc)
33 34 #error Unsupported architecture
34 35 #endif
35 36
|
↓ open down ↓ |
1 lines elided |
↑ open up ↑ |
36 37 /*
37 38 * expm1l(x)
38 39 *
39 40 * Table driven method
40 41 * Written by K.C. Ng, June 1995.
41 42 * Algorithm :
42 43 * 1. expm1(x) = x if x<2**-114
43 44 * 2. if |x| <= 0.0625 = 1/16, use approximation
44 45 * expm1(x) = x + x*P/(2-P)
45 46 * where
46 - * P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
47 + * P = x - z*(P1+z*(P2+z*(P3+z*(P4+z*(P5+z*P6+z*P7))))), z = x*x;
47 48 * (this formula is derived from
48 49 * 2-P+x = R = x*(exp(x)+1)/(exp(x)-1) ~ 2 + x*x/6 - x^4/360 + ...)
49 50 *
50 51 * P1 = 1.66666666666666666666666666666638500528074603030e-0001
51 52 * P2 = -2.77777777777777777777777759668391122822266551158e-0003
52 53 * P3 = 6.61375661375661375657437408890138814721051293054e-0005
53 54 * P4 = -1.65343915343915303310185228411892601606669528828e-0006
54 55 * P5 = 4.17535139755122945763580609663414647067443411178e-0008
55 56 * P6 = -1.05683795988668526689182102605260986731620026832e-0009
56 57 * P7 = 2.67544168821852702827123344217198187229611470514e-0011
57 58 *
58 59 * Accuracy: |R-x*(exp(x)+1)/(exp(x)-1)|<=2**-119.13
59 60 *
60 61 * 3. For 1/16 < |x| < 1.125, choose x(+-i) ~ +-(i+4.5)/64, i=0,..,67
61 62 * since
62 63 * exp(x) = exp(xi+(x-xi))= exp(xi)*exp((x-xi))
63 64 * we have
64 65 * expm1(x) = expm1(xi)+(exp(xi))*(expm1(x-xi))
65 66 * where
66 67 * |s=x-xi| <= 1/128
67 68 * and
68 69 * expm1(s)=2s/(2-R), R= s-s^2*(T1+s^2*(T2+s^2*(T3+s^2*(T4+s^2*T5))))
69 70 *
70 71 * T1 = 1.666666666666666666666666666660876387437e-1L,
71 72 * T2 = -2.777777777777777777777707812093173478756e-3L,
72 73 * T3 = 6.613756613756613482074280932874221202424e-5L,
73 74 * T4 = -1.653439153392139954169609822742235851120e-6L,
74 75 * T5 = 4.175314851769539751387852116610973796053e-8L;
75 76 *
76 77 * 4. For |x| >= 1.125, return exp(x)-1.
77 78 * (see algorithm for exp)
78 79 *
79 80 * Special cases:
80 81 * expm1l(INF) is INF, expm1l(NaN) is NaN;
81 82 * expm1l(-INF)= -1;
82 83 * for finite argument, only expm1l(0)=0 is exact.
83 84 *
84 85 * Accuracy:
85 86 * according to an error analysis, the error is always less than
86 87 * 2 ulp (unit in the last place).
87 88 *
88 89 * Misc. info.
89 90 * For 113 bit long double
90 91 * if x > 1.135652340629414394949193107797076342845e+4
91 92 * then expm1l(x) overflow;
92 93 *
|
↓ open down ↓ |
36 lines elided |
↑ open up ↑ |
93 94 * Constants:
94 95 * Only decimal values are given. We assume that the compiler will convert
95 96 * from decimal to binary accurately enough to produce the correct
96 97 * hexadecimal values.
97 98 */
98 99
99 100 #include "libm.h"
100 101
101 102 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
102 103 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,
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,
115 114 P1 = +1.66666666666666666666666666666638500528074603030e-0001L,
116 115 P2 = -2.77777777777777777777777759668391122822266551158e-0003L,
117 116 P3 = +6.61375661375661375657437408890138814721051293054e-0005L,
118 117 P4 = -1.65343915343915303310185228411892601606669528828e-0006L,
119 118 P5 = +4.17535139755122945763580609663414647067443411178e-0008L,
120 119 P6 = -1.05683795988668526689182102605260986731620026832e-0009L,
121 120 P7 = +2.67544168821852702827123344217198187229611470514e-0011L,
122 121 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
123 122 T1 = +1.666666666666666666666666666660876387437e-1L,
124 123 T2 = -2.777777777777777777777707812093173478756e-3L,
125 124 T3 = +6.613756613756613482074280932874221202424e-5L,
126 125 T4 = -1.653439153392139954169609822742235851120e-6L,
127 126 T5 = +4.175314851769539751387852116610973796053e-8L;
128 127
129 128 long double
130 -expm1l(long double x) {
129 +expm1l(long double x)
130 +{
131 131 int hx, ix, j, k, m;
132 132 long double t, r, s, w;
133 133
134 - hx = ((int *) &x)[HIXWORD];
134 + hx = ((int *)&x)[HIXWORD];
135 135 ix = hx & ~0x80000000;
136 +
136 137 if (ix >= 0x7fff0000) {
137 138 if (x != x)
138 - return (x + x); /* NaN */
139 + return (x + x); /* NaN */
140 +
139 141 if (x < zero)
140 - return (-one); /* -inf */
141 - return (x); /* +inf */
142 + return (-one); /* -inf */
143 +
144 + return (x); /* +inf */
142 145 }
143 - if (ix < 0x3fff4000) { /* |x| < 1.25 */
144 - if (ix < 0x3ffb0000) { /* |x| < 0.0625 */
146 +
147 + if (ix < 0x3fff4000) { /* |x| < 1.25 */
148 + if (ix < 0x3ffb0000) { /* |x| < 0.0625 */
145 149 if (ix < 0x3f8d0000) {
146 - if ((int) x == 0)
150 + if ((int)x == 0)
147 151 return (x); /* |x|<2^-114 */
148 152 }
153 +
149 154 t = x * x;
150 155 r = (x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t *
151 - (P5 + t * (P6 + t * P7)))))));
156 + (P5 + t * (P6 + t * P7)))))));
152 157 return (x + (x * r) / (two - r));
153 158 }
159 +
154 160 /* compute i = [64*x] */
155 161 m = 0x4009 - (ix >> 16);
156 162 j = ((ix & 0x0000ffff) | 0x10000) >> m; /* j=4,...,67 */
163 +
157 164 if (hx < 0)
158 165 j += 82; /* negative */
166 +
159 167 s = x - _TBL_expm1lx[j];
160 168 t = s * s;
161 169 r = s - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))));
162 170 r = (s + s) / (two - r);
163 171 w = _TBL_expm1l[j];
164 172 return (w + (w + one) * r);
165 173 }
174 +
166 175 if (hx > 0) {
167 176 if (x > ovflthreshold)
168 177 return (huge * huge);
169 - k = (int) (invln2_32 * (x + ln2_64));
178 +
179 + k = (int)(invln2_32 * (x + ln2_64));
170 180 } else {
171 181 if (x < -80.0)
172 182 return (tiny - x / x);
173 - k = (int) (invln2_32 * (x - ln2_64));
183 +
184 + k = (int)(invln2_32 * (x - ln2_64));
174 185 }
186 +
175 187 j = k & 0x1f;
176 188 m = k >> 5;
177 - t = (long double) k;
189 + t = (long double)k;
178 190 x = (x - t * ln2_32hi) - t * ln2_32lo;
179 191 t = x * x;
180 192 r = (x - t * (T1 + t * (T2 + t * (T3 + t * (T4 + t * T5))))) - two;
181 193 x = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (x + x)) / r -
182 - _TBL_expl_lo[j]);
194 + _TBL_expl_lo[j]);
183 195 return (scalbnl(x, m) - one);
184 196 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX