gonzui

    1: /* e_expf.c -- float version of e_exp.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;
   20: 
   21: static const float
   22: one     = 1.0,
   23: halF[2] = {0.5,-0.5,},
   24: twom100 = 7.8886090522e-31,      /* 2**-100=0x0d800000 */
   25: o_threshold=  8.8721679688e+01,  /* 0x42b17180 */
   26: u_threshold= -1.0397208405e+02,  /* 0xc2cff1b5 */
   27: ln2HI[2]   ={ 6.9313812256e-01,         /* 0x3f317180 */
   28:              -6.9313812256e-01,},      /* 0xbf317180 */
   29: ln2LO[2]   ={ 9.0580006145e-06,         /* 0x3717f7d1 */
   30:              -9.0580006145e-06,},      /* 0xb717f7d1 */
   31: invln2 =  1.4426950216e+00,             /* 0x3fb8aa3b */
   32: P1   =  1.6666667163e-01, /* 0x3e2aaaab */
   33: P2   = -2.7777778450e-03, /* 0xbb360b61 */
   34: P3   =  6.6137559770e-05, /* 0x388ab355 */
   35: P4   = -1.6533901999e-06, /* 0xb5ddea0e */
   36: P5   =  4.1381369442e-08; /* 0x3331bb4c */
   37: 
   38: float
   39: expf(float x)   /* default IEEE double exp */
   40: {
   41:         float y,hi,lo,c,t;
   42:         int32_t k,xsb;
   43:         u_int32_t hx;
   44: 
   45:         GET_FLOAT_WORD(hx,x);
   46:         xsb = (hx>>31)&1;              /* sign bit of x */
   47:         hx &= 0x7fffffff;              /* high word of |x| */
   48: 
   49:     /* filter out non-finite argument */
   50:         if(hx >= 0x42b17218) {                 /* if |x|>=88.721... */
   51:             if(hx>0x7f800000)
   52:                  return x+x;                   /* NaN */
   53:             if(hx==0x7f800000)
   54:                 return (xsb==0)? x:0.0;               /* exp(+-inf)={inf,0} */
   55:             if(x > o_threshold) return huge*huge; /* overflow */
   56:             if(x < u_threshold) return twom100*twom100; /* underflow */
   57:         }
   58: 
   59:     /* argument reduction */
   60:         if(hx > 0x3eb17218) {          /* if  |x| > 0.5 ln2 */ 
   61:             if(hx < 0x3F851592) {      /* and |x| < 1.5 ln2 */
   62:                 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
   63:             } else {
   64:                 k  = invln2*x+halF[xsb];
   65:                 t  = k;
   66:                 hi = x - t*ln2HI[0];  /* t*ln2HI is exact here */
   67:                 lo = t*ln2LO[0];
   68:             }
   69:             x  = hi - lo;
   70:         } 
   71:         else if(hx < 0x31800000)  {    /* when |x|<2**-28 */
   72:             if(huge+x>one) return one+x;/* trigger inexact */
   73:         }
   74:         else k = 0;
   75: 
   76:     /* x is now in primary range */
   77:         t  = x*x;
   78:         c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
   79:         if(k==0)       return one-((x*c)/(c-(float)2.0)-x); 
   80:         else           y = one-((lo-(x*c)/((float)2.0-c))-hi);
   81:         if(k >= -125) {
   82:             u_int32_t hy;
   83:             GET_FLOAT_WORD(hy,y);
   84:             SET_FLOAT_WORD(y,hy+(k<<23));      /* add k to y's exponent */
   85:             return y;
   86:         } else {
   87:             u_int32_t hy;
   88:             GET_FLOAT_WORD(hy,y);
   89:             SET_FLOAT_WORD(y,hy+((k+100)<<23));        /* add k to y's exponent */
   90:             return y*twom100;
   91:         }
   92: }