#include "math_libm.h"
#include "math_private.h"

Include dependency graph for e_exp.c:

Functions
double	__ieee754_exp (double x)

Variables
static const double	one = 1.0
static const double	halF [2] = {0.5,-0.5,}
static const double	huge = 1.0e+300
static const double	twom1000 = 9.33263618503218878990e-302
static const double	o_threshold = 7.09782712893383973096e+02
static const double	u_threshold = -7.45133219101941108420e+02
static const double	ln2HI [2]
static const double	ln2LO [2]
static const double	invln2 = 1.44269504088896338700e+00
static const double	P1 = 1.66666666666666019037e-01
static const double	P2 = -2.77777777770155933842e-03
static const double	P3 = 6.61375632143793436117e-05
static const double	P4 = -1.65339022054652515390e-06
static const double	P5 = 4.13813679705723846039e-08

Function Documentation

double __ieee754_exp ( double x )

Definition at line 96 of file e_exp.c.

References GET_HIGH_WORD, GET_LOW_WORD, halF, huge, invln2, k, ln2HI, ln2LO, o_threshold, one, P1, P2, P3, P4, P5, SET_HIGH_WORD, twom1000, and u_threshold.

{
        double y;
        double hi = 0.0;
        double lo = 0.0;
        double c;
        double t;
        int32_t k=0;
        int32_t xsb;
        u_int32_t hx;
        GET_HIGH_WORD(hx,x);
        xsb = (hx>>31)&1;               /* sign bit of x */
        hx &= 0x7fffffff;               /* high word of |x| */
    /* filter out non-finite argument */
        if(hx >= 0x40862E42) {                  /* if |x|>=709.78... */
            if(hx>=0x7ff00000) {
                u_int32_t lx;
                GET_LOW_WORD(lx,x);
                if(((hx&0xfffff)|lx)!=0)
                     return x+x;                /* NaN */
                else return (xsb==0)? x:0.0;    /* exp(+-inf)={inf,0} */
            }
                #if 1
                if(x > o_threshold) return huge*huge; /* overflow */
                #else  /* !!! FIXME: check this: "huge * huge" is a compiler warning, maybe they wanted +Inf? */
                if(x > o_threshold) return INFINITY; /* overflow */
                #endif
            if(x < u_threshold) return twom1000*twom1000; /* underflow */
        }
    /* argument reduction */
        if(hx > 0x3fd62e42) {           /* if  |x| > 0.5 ln2 */
            if(hx < 0x3FF0A2B2) {       /* and |x| < 1.5 ln2 */
                hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
            } else {
                k  = (int32_t) (invln2*x+halF[xsb]);
                t  = k;
                hi = x - t*ln2HI[0];    /* t*ln2HI is exact here */
                lo = t*ln2LO[0];
            }
            x  = hi - lo;
        }
        else if(hx < 0x3e300000)  {     /* when |x|<2**-28 */
            if(huge+x>one) return one+x;/* trigger inexact */
        }
        else k = 0;
    /* x is now in primary range */
        t  = x*x;
        c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
        if(k==0)        return one-((x*c)/(c-2.0)-x);
        else            y = one-((lo-(x*c)/(2.0-c))-hi);
        if(k >= -1021) {
            u_int32_t hy;
            GET_HIGH_WORD(hy,y);
            SET_HIGH_WORD(y,hy+(k<<20));        /* add k to y's exponent */
            return y;
        } else {
            u_int32_t hy;
            GET_HIGH_WORD(hy,y);
            SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */
            return y*twom1000;
        }
}