Back to index

glibc  2.9
Functions | Variables
s_expm1f.c File Reference
#include "math.h"
#include "math_private.h"

Go to the source code of this file.

Functions

float __expm1f (float x)

Variables

static const volatile float huge = 1.0e+30
static const volatile float tiny = 1.0e-30
static float one = 1.0
static float o_threshold = 8.8721679688e+01
static float ln2_hi = 6.9313812256e-01
static float ln2_lo = 9.0580006145e-06
static float invln2 = 1.4426950216e+00
static float Q1 = -3.3333335072e-02
static float Q2 = 1.5873016091e-03
static float Q3 = -7.9365076090e-05
static float Q4 = 4.0082177293e-06
static float Q5 = -2.0109921195e-07

Function Documentation

float __expm1f ( float  x)

Definition at line 46 of file s_expm1f.c.

{
       float y,hi,lo,c,t,e,hxs,hfx,r1;
       int32_t k,xsb;
       u_int32_t hx;

       GET_FLOAT_WORD(hx,x);
       xsb = hx&0x80000000;        /* sign bit of x */
       if(xsb==0) y=x; else y= -x; /* y = |x| */
       hx &= 0x7fffffff;           /* high word of |x| */

    /* filter out huge and non-finite argument */
       if(hx >= 0x4195b844) {                    /* if |x|>=27*ln2 */
           if(hx >= 0x42b17218) {         /* if |x|>=88.721... */
                if(hx>0x7f800000)
                  return x+x;       /* NaN */
              if(hx==0x7f800000)
                  return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
               if(x > o_threshold) return huge*huge; /* overflow */
           }
           if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */
              if(x+tiny<(float)0.0)       /* raise inexact */
              return tiny-one;     /* return -1 */
           }
       }

    /* argument reduction */
       if(hx > 0x3eb17218) {              /* if  |x| > 0.5 ln2 */
           if(hx < 0x3F851592) {   /* and |x| < 1.5 ln2 */
              if(xsb==0)
                  {hi = x - ln2_hi; lo =  ln2_lo;  k =  1;}
              else
                  {hi = x + ln2_hi; lo = -ln2_lo;  k = -1;}
           } else {
              k  = invln2*x+((xsb==0)?(float)0.5:(float)-0.5);
              t  = k;
              hi = x - t*ln2_hi;   /* t*ln2_hi is exact here */
              lo = t*ln2_lo;
           }
           x  = hi - lo;
           c  = (hi-x)-lo;
       }
       else if(hx < 0x33000000) {         /* when |x|<2**-25, return x */
           t = huge+x;      /* return x with inexact flags when x!=0 */
           return x - (t-(huge+x));
       }
       else k = 0;

    /* x is now in primary range */
       hfx = (float)0.5*x;
       hxs = x*hfx;
       r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
       t  = (float)3.0-r1*hfx;
       e  = hxs*((r1-t)/((float)6.0 - x*t));
       if(k==0) return x - (x*e-hxs);            /* c is 0 */
       else {
           e  = (x*(e-c)-c);
           e -= hxs;
           if(k== -1) return (float)0.5*(x-e)-(float)0.5;
           if(k==1) {
                     if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5));
                     else         return  one+(float)2.0*(x-e);
           }
           if (k <= -2 || k>56) {   /* suffice to return exp(x)-1 */
               int32_t i;
               y = one-(e-x);
              GET_FLOAT_WORD(i,y);
              SET_FLOAT_WORD(y,i+(k<<23));       /* add k to y's exponent */
               return y-one;
           }
           t = one;
           if(k<23) {
               int32_t i;
               SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */
                     y = t-(e-x);
              GET_FLOAT_WORD(i,y);
              SET_FLOAT_WORD(y,i+(k<<23));       /* add k to y's exponent */
          } else {
               int32_t i;
              SET_FLOAT_WORD(t,((0x7f-k)<<23));  /* 2^-k */
                     y = x-(e+t);
                     y += one;
              GET_FLOAT_WORD(i,y);
              SET_FLOAT_WORD(y,i+(k<<23));       /* add k to y's exponent */
           }
       }
       return y;
}

Here is the caller graph for this function:


Variable Documentation

const volatile float huge = 1.0e+30 [static]

Definition at line 23 of file s_expm1f.c.

float invln2 = 1.4426950216e+00 [static]

Definition at line 35 of file s_expm1f.c.

float ln2_hi = 6.9313812256e-01 [static]

Definition at line 33 of file s_expm1f.c.

float ln2_lo = 9.0580006145e-06 [static]

Definition at line 34 of file s_expm1f.c.

float o_threshold = 8.8721679688e+01 [static]

Definition at line 32 of file s_expm1f.c.

float one = 1.0 [static]

Definition at line 31 of file s_expm1f.c.

float Q1 = -3.3333335072e-02 [static]

Definition at line 37 of file s_expm1f.c.

float Q2 = 1.5873016091e-03 [static]

Definition at line 38 of file s_expm1f.c.

float Q3 = -7.9365076090e-05 [static]

Definition at line 39 of file s_expm1f.c.

float Q4 = 4.0082177293e-06 [static]

Definition at line 40 of file s_expm1f.c.

float Q5 = -2.0109921195e-07 [static]

Definition at line 41 of file s_expm1f.c.

const volatile float tiny = 1.0e-30 [static]

Definition at line 24 of file s_expm1f.c.