lib/msun/src/e_expf.c

   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 #ifndef lint
  17 static char rcsid[] = "$FreeBSD$";
  18 #endif
  19
  20 #include "math.h"
  21 #include "math_private.h"
  22
  23 static const float
  24 one     = 1.0,
  25 halF[2] = {0.5,-0.5,},
  26 huge    = 1.0e+30,
  27 twom100 = 7.8886090522e-31,      /* 2**-100=0x0d800000 */
  28 o_threshold=  8.8721679688e+01,  /* 0x42b17180 */
  29 u_threshold= -1.0397208405e+02,  /* 0xc2cff1b5 */
  30 invln2 =  1.4426950216e+00,             /* 0x3fb8aa3b */
  31 P1   =  1.6666667163e-01, /* 0x3e2aaaab */
  32 P2   = -2.7777778450e-03, /* 0xbb360b61 */
  33 P3   =  6.6137559770e-05, /* 0x388ab355 */
  34 P4   = -1.6533901999e-06, /* 0xb5ddea0e */
  35 P5   =  4.1381369442e-08; /* 0x3331bb4c */
  36 double ln2[2] = { 6.93147180369123816490e-01, -6.93147180369123816490e-01 };
  37
  38 float
  39 __ieee754_expf(float x) /* IEEE float exp */
  40 {
  41         float y,c,t;
  42         int32_t k=0,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;               /* |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                 x = x-ln2[xsb]; k = 1-xsb-xsb;
  63             } else {
  64                 k  = invln2*x+halF[xsb];
  65                 t  = k;
  66                 x = x - t*ln2[0];
  67             }
  68         }
  69         else if(hx < 0x31800000)  {     /* when |x|<2**-28 */
  70             if(huge+x>one) return one+x;/* trigger inexact */
  71         }
  72         else k = 0;
  73
  74     /* x is now in primary range */
  75         t  = x*x;
  76         c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
  77         y = one-(((double)x*c)/(c-2.0)-x);
  78         if(k==0)        return y;
  79         if(k >= -125) {
  80             u_int32_t hy;
  81             GET_FLOAT_WORD(hy,y);
  82             SET_FLOAT_WORD(y,hy+(k<<23));       /* add k to y's exponent */
  83             return y;
  84         } else {
  85             u_int32_t hy;
  86             GET_FLOAT_WORD(hy,y);
  87             SET_FLOAT_WORD(y,hy+((k+100)<<23)); /* add k to y's exponent */
  88             return y*twom100;
  89         }
  90 }