1 /* s_expm1f.c -- float version of s_expm1.c. 2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 */ 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 "math.h" 17 #include "math_private.h" 18 19 static const volatile float huge = 1.0e+30, tiny = 1.0e-30; 20 21 static const float 22 one = 1.0, 23 o_threshold = 8.8721679688e+01,/* 0x42b17180 */ 24 ln2_hi = 6.9313812256e-01,/* 0x3f317180 */ 25 ln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */ 26 invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */ 27 /* scaled coefficients related to expm1 */ 28 Q1 = -3.3333335072e-02, /* 0xbd088889 */ 29 Q2 = 1.5873016091e-03, /* 0x3ad00d01 */ 30 Q3 = -7.9365076090e-05, /* 0xb8a670cd */ 31 Q4 = 4.0082177293e-06, /* 0x36867e54 */ 32 Q5 = -2.0109921195e-07; /* 0xb457edbb */ 33 34 float 35 expm1f(float x) 36 { 37 float y,hi,lo,c,t,e,hxs,hfx,r1; 38 int32_t k,xsb; 39 u_int32_t hx; 40 41 GET_FLOAT_WORD(hx,x); 42 xsb = hx&0x80000000; /* sign bit of x */ 43 if(xsb==0) y=x; else y= -x; /* y = |x| */ 44 hx &= 0x7fffffff; /* high word of |x| */ 45 46 /* filter out huge and non-finite argument */ 47 if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ 48 if(hx >= 0x42b17218) { /* if |x|>=88.721... */ 49 if(hx>0x7f800000) 50 return x+x; /* NaN */ 51 if(hx==0x7f800000) 52 return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ 53 if(x > o_threshold) return huge*huge; /* overflow */ 54 } 55 if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ 56 if(x+tiny<(float)0.0) /* raise inexact */ 57 return tiny-one; /* return -1 */ 58 } 59 } 60 61 /* argument reduction */ 62 if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ 63 if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ 64 if(xsb==0) 65 {hi = x - ln2_hi; lo = ln2_lo; k = 1;} 66 else 67 {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} 68 } else { 69 k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5); 70 t = k; 71 hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ 72 lo = t*ln2_lo; 73 } 74 x = hi - lo; 75 c = (hi-x)-lo; 76 } 77 else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ 78 t = huge+x; /* return x with inexact flags when x!=0 */ 79 return x - (t-(huge+x)); 80 } 81 else k = 0; 82 83 /* x is now in primary range */ 84 hfx = (float)0.5*x; 85 hxs = x*hfx; 86 r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); 87 t = (float)3.0-r1*hfx; 88 e = hxs*((r1-t)/((float)6.0 - x*t)); 89 if(k==0) return x - (x*e-hxs); /* c is 0 */ 90 else { 91 e = (x*(e-c)-c); 92 e -= hxs; 93 if(k== -1) return (float)0.5*(x-e)-(float)0.5; 94 if(k==1) { 95 if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); 96 else return one+(float)2.0*(x-e); 97 } 98 if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ 99 int32_t i; 100 y = one-(e-x); 101 GET_FLOAT_WORD(i,y); 102 SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ 103 return y-one; 104 } 105 t = one; 106 if(k<23) { 107 int32_t i; 108 SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ 109 y = t-(e-x); 110 GET_FLOAT_WORD(i,y); 111 SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ 112 } else { 113 int32_t i; 114 SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ 115 y = x-(e+t); 116 y += one; 117 GET_FLOAT_WORD(i,y); 118 SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ 119 } 120 } 121 return y; 122 } 123