| 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 |