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