OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.c */ |
2 /* | 2 /* |
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 */ | 4 */ |
5 /* | 5 /* |
6 * ==================================================== | 6 * ==================================================== |
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 * | 8 * |
9 * Developed at SunPro, a Sun Microsystems, Inc. business. | 9 * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 * Permission to use, copy, modify, and distribute this | 10 * Permission to use, copy, modify, and distribute this |
11 * software is freely granted, provided that this notice | 11 * software is freely granted, provided that this notice |
12 * is preserved. | 12 * is preserved. |
13 * ==================================================== | 13 * ==================================================== |
14 */ | 14 */ |
15 | 15 |
16 #include "libm.h" | 16 #include "libm.h" |
17 | 17 |
18 static const float | 18 static const float o_threshold = 8.8721679688e+01, /* 0x42b17180 */ |
19 o_threshold = 8.8721679688e+01, /* 0x42b17180 */ | 19 ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ |
20 ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ | 20 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ |
21 ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ | 21 invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */ |
22 invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */ | 22 /* |
23 /* | 23 * Domain [-0.34568, 0.34568]
, range ~[-6.694e-10, 6.696e-10]: |
24 * Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]: | 24 * |6 / x * (1 + 2 * (1 / (ex
p(x) - 1) - 1 / x)) - q(x)| < 2**-30.04 |
25 * |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04 | 25 * Scaled coefficients: Qn_he
re = 2**n * Qn_for_q (see s_expm1.c): |
26 * Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c): | 26 */ |
27 */ | 27 Q1 = -3.3333212137e-2, /* -0x888868.0p-28 */ |
28 Q1 = -3.3333212137e-2, /* -0x888868.0p-28 */ | 28 Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */ |
29 Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */ | |
30 | 29 |
31 float expm1f(float x) | 30 float expm1f(float x) { |
32 { | 31 float_t y, hi, lo, c, t, e, hxs, hfx, r1, twopk; |
33 » float_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk; | 32 union { |
34 » union {float f; uint32_t i;} u = {x}; | 33 float f; |
35 » uint32_t hx = u.i & 0x7fffffff; | 34 uint32_t i; |
36 » int k, sign = u.i >> 31; | 35 } u = {x}; |
| 36 uint32_t hx = u.i & 0x7fffffff; |
| 37 int k, sign = u.i >> 31; |
37 | 38 |
38 » /* filter out huge and non-finite argument */ | 39 /* filter out huge and non-finite argument */ |
39 » if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ | 40 if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ |
40 » » if (hx > 0x7f800000) /* NaN */ | 41 if (hx > 0x7f800000) /* NaN */ |
41 » » » return x; | 42 return x; |
42 » » if (sign) | 43 if (sign) |
43 » » » return -1; | 44 return -1; |
44 » » if (x > o_threshold) { | 45 if (x > o_threshold) { |
45 » » » x *= 0x1p127f; | 46 x *= 0x1p127f; |
46 » » » return x; | 47 return x; |
47 » » } | 48 } |
48 » } | 49 } |
49 | 50 |
50 » /* argument reduction */ | 51 /* argument reduction */ |
51 » if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ | 52 if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ |
52 » » if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ | 53 if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ |
53 » » » if (!sign) { | 54 if (!sign) { |
54 » » » » hi = x - ln2_hi; | 55 hi = x - ln2_hi; |
55 » » » » lo = ln2_lo; | 56 lo = ln2_lo; |
56 » » » » k = 1; | 57 k = 1; |
57 » » » } else { | 58 } else { |
58 » » » » hi = x + ln2_hi; | 59 hi = x + ln2_hi; |
59 » » » » lo = -ln2_lo; | 60 lo = -ln2_lo; |
60 » » » » k = -1; | 61 k = -1; |
61 » » » } | 62 } |
62 » » } else { | 63 } else { |
63 » » » k = invln2*x + (sign ? -0.5f : 0.5f); | 64 k = invln2 * x + (sign ? -0.5f : 0.5f); |
64 » » » t = k; | 65 t = k; |
65 » » » hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ | 66 hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ |
66 » » » lo = t*ln2_lo; | 67 lo = t * ln2_lo; |
67 » » } | 68 } |
68 » » x = hi-lo; | 69 x = hi - lo; |
69 » » c = (hi-x)-lo; | 70 c = (hi - x) - lo; |
70 » } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */ | 71 } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */ |
71 » » if (hx < 0x00800000) | 72 if (hx < 0x00800000) |
72 » » » FORCE_EVAL(x*x); | 73 FORCE_EVAL(x * x); |
73 » » return x; | 74 return x; |
74 » } else | 75 } else |
75 » » k = 0; | 76 k = 0; |
76 | 77 |
77 » /* x is now in primary range */ | 78 /* x is now in primary range */ |
78 » hfx = 0.5f*x; | 79 hfx = 0.5f * x; |
79 » hxs = x*hfx; | 80 hxs = x * hfx; |
80 » r1 = 1.0f+hxs*(Q1+hxs*Q2); | 81 r1 = 1.0f + hxs * (Q1 + hxs * Q2); |
81 » t = 3.0f - r1*hfx; | 82 t = 3.0f - r1 * hfx; |
82 » e = hxs*((r1-t)/(6.0f - x*t)); | 83 e = hxs * ((r1 - t) / (6.0f - x * t)); |
83 » if (k == 0) /* c is 0 */ | 84 if (k == 0) /* c is 0 */ |
84 » » return x - (x*e-hxs); | 85 return x - (x * e - hxs); |
85 » e = x*(e-c) - c; | 86 e = x * (e - c) - c; |
86 » e -= hxs; | 87 e -= hxs; |
87 » /* exp(x) ~ 2^k (x_reduced - e + 1) */ | 88 /* exp(x) ~ 2^k (x_reduced - e + 1) */ |
88 » if (k == -1) | 89 if (k == -1) |
89 » » return 0.5f*(x-e) - 0.5f; | 90 return 0.5f * (x - e) - 0.5f; |
90 » if (k == 1) { | 91 if (k == 1) { |
91 » » if (x < -0.25f) | 92 if (x < -0.25f) |
92 » » » return -2.0f*(e-(x+0.5f)); | 93 return -2.0f * (e - (x + 0.5f)); |
93 » » return 1.0f + 2.0f*(x-e); | 94 return 1.0f + 2.0f * (x - e); |
94 » } | 95 } |
95 » u.i = (0x7f+k)<<23; /* 2^k */ | 96 u.i = (0x7f + k) << 23; /* 2^k */ |
96 » twopk = u.f; | 97 twopk = u.f; |
97 » if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ | 98 if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ |
98 » » y = x - e + 1.0f; | 99 y = x - e + 1.0f; |
99 » » if (k == 128) | 100 if (k == 128) |
100 » » » y = y*2.0f*0x1p127f; | 101 y = y * 2.0f * 0x1p127f; |
101 » » else | 102 else |
102 » » » y = y*twopk; | 103 y = y * twopk; |
103 » » return y - 1.0f; | 104 return y - 1.0f; |
104 » } | 105 } |
105 » u.i = (0x7f-k)<<23; /* 2^-k */ | 106 u.i = (0x7f - k) << 23; /* 2^-k */ |
106 » if (k < 23) | 107 if (k < 23) |
107 » » y = (x-e+(1-u.f))*twopk; | 108 y = (x - e + (1 - u.f)) * twopk; |
108 » else | 109 else |
109 » » y = (x-(e+u.f)+1)*twopk; | 110 y = (x - (e + u.f) + 1) * twopk; |
110 » return y; | 111 return y; |
111 } | 112 } |
OLD | NEW |