s_expm1f.c revision 50476
175584Sru/* s_expm1f.c -- float version of s_expm1.c. 275584Sru * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 375584Sru */ 475584Sru 575584Sru/* 675584Sru * ==================================================== 775584Sru * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 875584Sru * 975584Sru * Developed at SunPro, a Sun Microsystems, Inc. business. 1075584Sru * Permission to use, copy, modify, and distribute this 1175584Sru * software is freely granted, provided that this notice 1275584Sru * is preserved. 1375584Sru * ==================================================== 1475584Sru */ 1575584Sru 1675584Sru#ifndef lint 1775584Srustatic char rcsid[] = "$FreeBSD: head/lib/msun/src/s_expm1f.c 50476 1999-08-28 00:22:10Z peter $"; 1875584Sru#endif 1975584Sru 2075584Sru#include "math.h" 2175584Sru#include "math_private.h" 2275584Sru 2375584Sru#ifdef __STDC__ 2475584Srustatic const float 2575584Sru#else 2675584Srustatic float 2775584Sru#endif 2875584Sruone = 1.0, 2975584Sruhuge = 1.0e+30, 3075584Srutiny = 1.0e-30, 3175584Sruo_threshold = 8.8721679688e+01,/* 0x42b17180 */ 3275584Sruln2_hi = 6.9313812256e-01,/* 0x3f317180 */ 3375584Sruln2_lo = 9.0580006145e-06,/* 0x3717f7d1 */ 3475584Sruinvln2 = 1.4426950216e+00,/* 0x3fb8aa3b */ 3575584Sru /* scaled coefficients related to expm1 */ 3675584SruQ1 = -3.3333335072e-02, /* 0xbd088889 */ 3775584SruQ2 = 1.5873016091e-03, /* 0x3ad00d01 */ 3875584SruQ3 = -7.9365076090e-05, /* 0xb8a670cd */ 3975584SruQ4 = 4.0082177293e-06, /* 0x36867e54 */ 4075584SruQ5 = -2.0109921195e-07; /* 0xb457edbb */ 4175584Sru 4275584Sru#ifdef __STDC__ 4375584Sru float expm1f(float x) 4475584Sru#else 4575584Sru float expm1f(x) 4675584Sru float x; 4775584Sru#endif 4875584Sru{ 4975584Sru float y,hi,lo,c,t,e,hxs,hfx,r1; 5075584Sru int32_t k,xsb; 5175584Sru u_int32_t hx; 5275584Sru 5375584Sru GET_FLOAT_WORD(hx,x); 5475584Sru xsb = hx&0x80000000; /* sign bit of x */ 5575584Sru if(xsb==0) y=x; else y= -x; /* y = |x| */ 5675584Sru hx &= 0x7fffffff; /* high word of |x| */ 5775584Sru 5875584Sru /* filter out huge and non-finite argument */ 5975584Sru if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ 6075584Sru if(hx >= 0x42b17218) { /* if |x|>=88.721... */ 6175584Sru if(hx>0x7f800000) 6275584Sru return x+x; /* NaN */ 6375584Sru if(hx==0x7f800000) 6475584Sru return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ 6575584Sru if(x > o_threshold) return huge*huge; /* overflow */ 6675584Sru } 6775584Sru if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ 6875584Sru if(x+tiny<(float)0.0) /* raise inexact */ 6975584Sru return tiny-one; /* return -1 */ 7075584Sru } 7175584Sru } 7275584Sru 7375584Sru /* argument reduction */ 7475584Sru if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ 7575584Sru if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ 7675584Sru if(xsb==0) 7775584Sru {hi = x - ln2_hi; lo = ln2_lo; k = 1;} 7875584Sru else 7975584Sru {hi = x + ln2_hi; lo = -ln2_lo; k = -1;} 8075584Sru } else { 8175584Sru k = invln2*x+((xsb==0)?(float)0.5:(float)-0.5); 8275584Sru t = k; 8375584Sru hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ 8475584Sru lo = t*ln2_lo; 8575584Sru } 8675584Sru x = hi - lo; 8775584Sru c = (hi-x)-lo; 8875584Sru } 8975584Sru else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ 9075584Sru t = huge+x; /* return x with inexact flags when x!=0 */ 9175584Sru return x - (t-(huge+x)); 9275584Sru } 9375584Sru else k = 0; 9475584Sru 9575584Sru /* x is now in primary range */ 9675584Sru hfx = (float)0.5*x; 9775584Sru hxs = x*hfx; 9875584Sru r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))); 9975584Sru t = (float)3.0-r1*hfx; 10075584Sru e = hxs*((r1-t)/((float)6.0 - x*t)); 10175584Sru if(k==0) return x - (x*e-hxs); /* c is 0 */ 10275584Sru else { 10375584Sru e = (x*(e-c)-c); 10475584Sru e -= hxs; 10575584Sru if(k== -1) return (float)0.5*(x-e)-(float)0.5; 10675584Sru if(k==1) 10775584Sru if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); 10875584Sru else return one+(float)2.0*(x-e); 10975584Sru if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ 11075584Sru int32_t i; 11175584Sru y = one-(e-x); 11275584Sru GET_FLOAT_WORD(i,y); 11375584Sru SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ 11475584Sru return y-one; 11575584Sru } 11675584Sru t = one; 11775584Sru if(k<23) { 11875584Sru int32_t i; 11975584Sru SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ 12075584Sru y = t-(e-x); 12175584Sru GET_FLOAT_WORD(i,y); 12275584Sru SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ 12375584Sru } else { 12475584Sru int32_t i; 12575584Sru SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ 12675584Sru y = x-(e+t); 12775584Sru y += one; 12875584Sru GET_FLOAT_WORD(i,y); 12975584Sru SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ 13075584Sru } 13175584Sru } 13275584Sru return y; 13375584Sru} 13475584Sru