gonzui

    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:             if (k <= -2 || k>56) {   /* suffice to return exp(x)-1 */
   98:                 int32_t i;
   99:                 y = one-(e-x);
  100:                 GET_FLOAT_WORD(i,y);
  101:                 SET_FLOAT_WORD(y,i+(k<<23));  /* add k to y's exponent */
  102:                 return y-one;
  103:             }
  104:             t = one;
  105:             if(k<23) {
  106:                 int32_t i;
  107:                 SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
  108:                        y = t-(e-x);
  109:                 GET_FLOAT_WORD(i,y);
  110:                 SET_FLOAT_WORD(y,i+(k<<23));  /* add k to y's exponent */
  111:            } else {
  112:                 int32_t i;
  113:                 SET_FLOAT_WORD(t,((0x7f-k)<<23));     /* 2^-k */
  114:                        y = x-(e+t);
  115:                        y += one;
  116:                 GET_FLOAT_WORD(i,y);
  117:                 SET_FLOAT_WORD(y,i+(k<<23));  /* add k to y's exponent */
  118:             }
  119:         }
  120:         return y;
  121: }