e_expf.c revision 175468
12116Sjkh/* e_expf.c -- float version of e_exp.c. 22116Sjkh * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 32116Sjkh */ 42116Sjkh 52116Sjkh/* 62116Sjkh * ==================================================== 72116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 82116Sjkh * 92116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 102116Sjkh * Permission to use, copy, modify, and distribute this 118870Srgrimes * software is freely granted, provided that this notice 122116Sjkh * is preserved. 132116Sjkh * ==================================================== 142116Sjkh */ 152116Sjkh 162116Sjkh#ifndef lint 1750476Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/e_expf.c 175468 2008-01-18 22:19:04Z das $"; 182116Sjkh#endif 192116Sjkh 202116Sjkh#include "math.h" 212116Sjkh#include "math_private.h" 222116Sjkh 232116Sjkhstatic const float 242116Sjkhone = 1.0, 252116SjkhhalF[2] = {0.5,-0.5,}, 262116Sjkhhuge = 1.0e+30, 272116Sjkho_threshold= 8.8721679688e+01, /* 0x42b17180 */ 282116Sjkhu_threshold= -1.0397208405e+02, /* 0xc2cff1b5 */ 29152947Sbdeln2HI[2] ={ 6.9314575195e-01, /* 0x3f317200 */ 30152947Sbde -6.9314575195e-01,}, /* 0xbf317200 */ 31152947Sbdeln2LO[2] ={ 1.4286067653e-06, /* 0x35bfbe8e */ 32152947Sbde -1.4286067653e-06,}, /* 0xb5bfbe8e */ 332116Sjkhinvln2 = 1.4426950216e+00, /* 0x3fb8aa3b */ 342116SjkhP1 = 1.6666667163e-01, /* 0x3e2aaaab */ 352116SjkhP2 = -2.7777778450e-03, /* 0xbb360b61 */ 362116SjkhP3 = 6.6137559770e-05, /* 0x388ab355 */ 372116SjkhP4 = -1.6533901999e-06, /* 0xb5ddea0e */ 382116SjkhP5 = 4.1381369442e-08; /* 0x3331bb4c */ 392116Sjkh 40175468Sdasstatic volatile float twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */ 41175468Sdas 4297407Salfredfloat 43142369Sdas__ieee754_expf(float x) /* default IEEE double exp */ 442116Sjkh{ 45142369Sdas float y,hi=0.0,lo=0.0,c,t; 4617141Sjkh int32_t k=0,xsb; 472116Sjkh u_int32_t hx; 482116Sjkh 492116Sjkh GET_FLOAT_WORD(hx,x); 502116Sjkh xsb = (hx>>31)&1; /* sign bit of x */ 51142369Sdas hx &= 0x7fffffff; /* high word of |x| */ 522116Sjkh 532116Sjkh /* filter out non-finite argument */ 542116Sjkh if(hx >= 0x42b17218) { /* if |x|>=88.721... */ 552116Sjkh if(hx>0x7f800000) 562116Sjkh return x+x; /* NaN */ 572116Sjkh if(hx==0x7f800000) 582116Sjkh return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ 592116Sjkh if(x > o_threshold) return huge*huge; /* overflow */ 602116Sjkh if(x < u_threshold) return twom100*twom100; /* underflow */ 612116Sjkh } 622116Sjkh 632116Sjkh /* argument reduction */ 648870Srgrimes if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ 652116Sjkh if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ 66142369Sdas hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; 672116Sjkh } else { 682116Sjkh k = invln2*x+halF[xsb]; 692116Sjkh t = k; 70142369Sdas hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */ 71142369Sdas lo = t*ln2LO[0]; 722116Sjkh } 73142369Sdas x = hi - lo; 748870Srgrimes } 752116Sjkh else if(hx < 0x31800000) { /* when |x|<2**-28 */ 762116Sjkh if(huge+x>one) return one+x;/* trigger inexact */ 772116Sjkh } 782116Sjkh else k = 0; 792116Sjkh 802116Sjkh /* x is now in primary range */ 812116Sjkh t = x*x; 822116Sjkh c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); 83142369Sdas if(k==0) return one-((x*c)/(c-(float)2.0)-x); 84142369Sdas else y = one-((lo-(x*c)/((float)2.0-c))-hi); 852116Sjkh if(k >= -125) { 862116Sjkh u_int32_t hy; 872116Sjkh GET_FLOAT_WORD(hy,y); 882116Sjkh SET_FLOAT_WORD(y,hy+(k<<23)); /* add k to y's exponent */ 892116Sjkh return y; 902116Sjkh } else { 912116Sjkh u_int32_t hy; 922116Sjkh GET_FLOAT_WORD(hy,y); 932116Sjkh SET_FLOAT_WORD(y,hy+((k+100)<<23)); /* add k to y's exponent */ 942116Sjkh return y*twom100; 952116Sjkh } 962116Sjkh} 97