Back to index

glibc  2.9
Defines | Functions | Variables
e_exp2.c File Reference
#include <stdlib.h>
#include <float.h>
#include <ieee754.h>
#include <math.h>
#include <fenv.h>
#include <inttypes.h>
#include <math_private.h>
#include "t_exp2.h"

Go to the source code of this file.

Defines

#define _GNU_SOURCE

Functions

double __ieee754_exp2 (double x)

Variables

static const volatile double TWO1023 = 8.988465674311579539e+307
static const volatile double TWOM1000 = 9.3326361850321887899e-302

Define Documentation

#define _GNU_SOURCE

Definition at line 28 of file e_exp2.c.


Function Documentation

double __ieee754_exp2 ( double  x)

Definition at line 49 of file e_exp2.c.

{
  static const double himark = (double) DBL_MAX_EXP;
  static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);

  /* Check for usual case.  */
  if (isless (x, himark) && isgreaterequal (x, lomark))
    {
      static const double THREEp42 = 13194139533312.0;
      int tval, unsafe;
      double rx, x22, result;
      union ieee754_double ex2_u, scale_u;
      fenv_t oldenv;

      feholdexcept (&oldenv);
#ifdef FE_TONEAREST
      /* If we don't have this, it's too bad.  */
      fesetround (FE_TONEAREST);
#endif

      /* 1. Argument reduction.
        Choose integers ex, -256 <= t < 256, and some real
        -1/1024 <= x1 <= 1024 so that
        x = ex + t/512 + x1.

        First, calculate rx = ex + t/512.  */
      rx = x + THREEp42;
      rx -= THREEp42;
      x -= rx;  /* Compute x=x1. */
      /* Compute tval = (ex*512 + t)+256.
        Now, t = (tval mod 512)-256 and ex=tval/512  [that's mod, NOT %; and
        /-round-to-nearest not the usual c integer /].  */
      tval = (int) (rx * 512.0 + 256.0);

      /* 2. Adjust for accurate table entry.
        Find e so that
        x = ex + t/512 + e + x2
        where -1e6 < e < 1e6, and
        (double)(2^(t/512+e))
        is accurate to one part in 2^-64.  */

      /* 'tval & 511' is the same as 'tval%512' except that it's always
        positive.
        Compute x = x2.  */
      x -= exp2_deltatable[tval & 511];

      /* 3. Compute ex2 = 2^(t/512+e+ex).  */
      ex2_u.d = exp2_accuratetable[tval & 511];
      tval >>= 9;
      unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
      ex2_u.ieee.exponent += tval >> unsafe;
      scale_u.d = 1.0;
      scale_u.ieee.exponent += tval - (tval >> unsafe);

      /* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
        with maximum error in [-2^-10-2^-30,2^-10+2^-30]
        less than 10^-19.  */

      x22 = (((.0096181293647031180
              * x + .055504110254308625)
             * x + .240226506959100583)
            * x + .69314718055994495) * ex2_u.d;

      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
      fesetenv (&oldenv);

      result = x22 * x + ex2_u.d;

      if (!unsafe)
       return result;
      else
       return result * scale_u.d;
    }
  /* Exceptional cases:  */
  else if (isless (x, himark))
    {
      if (__isinf (x))
       /* e^-inf == 0, with no error.  */
       return 0;
      else
       /* Underflow */
       return TWOM1000 * TWOM1000;
    }
  else
    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
    return TWO1023*x;
}

Here is the call graph for this function:

Here is the caller graph for this function:


Variable Documentation

const volatile double TWO1023 = 8.988465674311579539e+307 [static]

Definition at line 45 of file e_exp2.c.

const volatile double TWOM1000 = 9.3326361850321887899e-302 [static]

Definition at line 46 of file e_exp2.c.