OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.c */ |
2 /* | 2 /* |
3 * ==================================================== | 3 * ==================================================== |
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 * | 5 * |
6 * Developed at SunPro, a Sun Microsystems, Inc. business. | 6 * Developed at SunPro, a Sun Microsystems, Inc. business. |
7 * Permission to use, copy, modify, and distribute this | 7 * Permission to use, copy, modify, and distribute this |
8 * software is freely granted, provided that this notice | 8 * software is freely granted, provided that this notice |
9 * is preserved. | 9 * is preserved. |
10 * ==================================================== | 10 * ==================================================== |
(...skipping 88 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
99 * | 99 * |
100 * Constants: | 100 * Constants: |
101 * The hexadecimal values are the intended ones for the following | 101 * The hexadecimal values are the intended ones for the following |
102 * constants. The decimal values may be used, provided that the | 102 * constants. The decimal values may be used, provided that the |
103 * compiler will convert from decimal to binary accurately enough | 103 * compiler will convert from decimal to binary accurately enough |
104 * to produce the hexadecimal values shown. | 104 * to produce the hexadecimal values shown. |
105 */ | 105 */ |
106 | 106 |
107 #include "libm.h" | 107 #include "libm.h" |
108 | 108 |
109 static const double | 109 static const double o_threshold = |
110 o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ | 110 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ |
111 ln2_hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ | 111 ln2_hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ |
112 ln2_lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ | 112 ln2_lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ |
113 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ | 113 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ |
114 /* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */ | 114 /* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: |
115 Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */ | 115 */ |
116 Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ | 116 Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */ |
117 Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ | 117 Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ |
118 Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ | 118 Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ |
119 Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ | 119 Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ |
| 120 Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */ |
120 | 121 |
121 double expm1(double x) | 122 double expm1(double x) { |
122 { | 123 double_t y, hi, lo, c, t, e, hxs, hfx, r1, twopk; |
123 » double_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk; | 124 union { |
124 » union {double f; uint64_t i;} u = {x}; | 125 double f; |
125 » uint32_t hx = u.i>>32 & 0x7fffffff; | 126 uint64_t i; |
126 » int k, sign = u.i>>63; | 127 } u = {x}; |
| 128 uint32_t hx = u.i >> 32 & 0x7fffffff; |
| 129 int k, sign = u.i >> 63; |
127 | 130 |
128 » /* filter out huge and non-finite argument */ | 131 /* filter out huge and non-finite argument */ |
129 » if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ | 132 if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ |
130 » » if (isnan(x)) | 133 if (isnan(x)) |
131 » » » return x; | 134 return x; |
132 » » if (sign) | 135 if (sign) |
133 » » » return -1; | 136 return -1; |
134 » » if (x > o_threshold) { | 137 if (x > o_threshold) { |
135 » » » x *= 0x1p1023; | 138 x *= 0x1p1023; |
136 » » » return x; | 139 return x; |
137 » » } | 140 } |
138 » } | 141 } |
139 | 142 |
140 » /* argument reduction */ | 143 /* argument reduction */ |
141 » if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ | 144 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ |
142 » » if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ | 145 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ |
143 » » » if (!sign) { | 146 if (!sign) { |
144 » » » » hi = x - ln2_hi; | 147 hi = x - ln2_hi; |
145 » » » » lo = ln2_lo; | 148 lo = ln2_lo; |
146 » » » » k = 1; | 149 k = 1; |
147 » » » } else { | 150 } else { |
148 » » » » hi = x + ln2_hi; | 151 hi = x + ln2_hi; |
149 » » » » lo = -ln2_lo; | 152 lo = -ln2_lo; |
150 » » » » k = -1; | 153 k = -1; |
151 » » » } | 154 } |
152 » » } else { | 155 } else { |
153 » » » k = invln2*x + (sign ? -0.5 : 0.5); | 156 k = invln2 * x + (sign ? -0.5 : 0.5); |
154 » » » t = k; | 157 t = k; |
155 » » » hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ | 158 hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ |
156 » » » lo = t*ln2_lo; | 159 lo = t * ln2_lo; |
157 » » } | 160 } |
158 » » x = hi-lo; | 161 x = hi - lo; |
159 » » c = (hi-x)-lo; | 162 c = (hi - x) - lo; |
160 » } else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */ | 163 } else if (hx < 0x3c900000) { /* |x| < 2**-54, return x */ |
161 » » if (hx < 0x00100000) | 164 if (hx < 0x00100000) |
162 » » » FORCE_EVAL((float)x); | 165 FORCE_EVAL((float)x); |
163 » » return x; | 166 return x; |
164 » } else | 167 } else |
165 » » k = 0; | 168 k = 0; |
166 | 169 |
167 » /* x is now in primary range */ | 170 /* x is now in primary range */ |
168 » hfx = 0.5*x; | 171 hfx = 0.5 * x; |
169 » hxs = x*hfx; | 172 hxs = x * hfx; |
170 » r1 = 1.0+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); | 173 r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); |
171 » t = 3.0-r1*hfx; | 174 t = 3.0 - r1 * hfx; |
172 » e = hxs*((r1-t)/(6.0 - x*t)); | 175 e = hxs * ((r1 - t) / (6.0 - x * t)); |
173 » if (k == 0) /* c is 0 */ | 176 if (k == 0) /* c is 0 */ |
174 » » return x - (x*e-hxs); | 177 return x - (x * e - hxs); |
175 » e = x*(e-c) - c; | 178 e = x * (e - c) - c; |
176 » e -= hxs; | 179 e -= hxs; |
177 » /* exp(x) ~ 2^k (x_reduced - e + 1) */ | 180 /* exp(x) ~ 2^k (x_reduced - e + 1) */ |
178 » if (k == -1) | 181 if (k == -1) |
179 » » return 0.5*(x-e) - 0.5; | 182 return 0.5 * (x - e) - 0.5; |
180 » if (k == 1) { | 183 if (k == 1) { |
181 » » if (x < -0.25) | 184 if (x < -0.25) |
182 » » » return -2.0*(e-(x+0.5)); | 185 return -2.0 * (e - (x + 0.5)); |
183 » » return 1.0+2.0*(x-e); | 186 return 1.0 + 2.0 * (x - e); |
184 » } | 187 } |
185 » u.i = (uint64_t)(0x3ff + k)<<52; /* 2^k */ | 188 u.i = (uint64_t)(0x3ff + k) << 52; /* 2^k */ |
186 » twopk = u.f; | 189 twopk = u.f; |
187 » if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ | 190 if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ |
188 » » y = x - e + 1.0; | 191 y = x - e + 1.0; |
189 » » if (k == 1024) | 192 if (k == 1024) |
190 » » » y = y*2.0*0x1p1023; | 193 y = y * 2.0 * 0x1p1023; |
191 » » else | 194 else |
192 » » » y = y*twopk; | 195 y = y * twopk; |
193 » » return y - 1.0; | 196 return y - 1.0; |
194 » } | 197 } |
195 » u.i = (uint64_t)(0x3ff - k)<<52; /* 2^-k */ | 198 u.i = (uint64_t)(0x3ff - k) << 52; /* 2^-k */ |
196 » if (k < 20) | 199 if (k < 20) |
197 » » y = (x-e+(1-u.f))*twopk; | 200 y = (x - e + (1 - u.f)) * twopk; |
198 » else | 201 else |
199 » » y = (x-(e+u.f)+1)*twopk; | 202 y = (x - (e + u.f) + 1) * twopk; |
200 » return y; | 203 return y; |
201 } | 204 } |
OLD | NEW |