Back to index

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

Go to the source code of this file.

Defines

#define _GNU_SOURCE
#define himark   C[0]
#define lomark   C[1]
#define THREEp96   C[2]
#define THREEp103   C[3]
#define THREEp111   C[4]
#define M_1_LN2   C[5]
#define M_LN2_0   C[6]
#define M_LN2_1   C[7]
#define TINY   C[8]
#define TWO16383   C[9]
#define TWO8   C[10]
#define TWO15   C[11]
#define P1   C[12]
#define P2   C[13]
#define P3   C[14]
#define P4   C[15]
#define P5   C[16]
#define P6   C[17]

Functions

long double __ieee754_expl (long double x)

Variables

static const long double C []

Define Documentation

#define _GNU_SOURCE

Definition at line 61 of file e_expl.c.

#define himark   C[0]
#define lomark   C[1]
#define M_1_LN2   C[5]
#define M_LN2_0   C[6]
#define M_LN2_1   C[7]
#define P1   C[12]
#define P2   C[13]
#define P3   C[14]
#define P4   C[15]
#define P5   C[16]
#define P6   C[17]
#define THREEp103   C[3]
#define THREEp111   C[4]
#define THREEp96   C[2]
#define TINY   C[8]
#define TWO15   C[11]
#define TWO16383   C[9]
#define TWO8   C[10]

Function Documentation

long double __ieee754_expl ( long double  x)

Definition at line 137 of file e_expl.c.

{
  /* Check for usual case.  */
  if (isless (x, himark) && isgreater (x, lomark))
    {
      int tval1, tval2, unsafe, n_i;
      long double x22, n, t, result, xl;
      union ieee854_long_double ex2_u, scale_u;
      fenv_t oldenv;

      feholdexcept (&oldenv);
#ifdef FE_TONEAREST
      fesetround (FE_TONEAREST);
#endif

      /* Calculate n.  */
      n = x * M_1_LN2 + THREEp111;
      n -= THREEp111;
      x = x - n * M_LN2_0;
      xl = n * M_LN2_1;

      /* Calculate t/256.  */
      t = x + THREEp103;
      t -= THREEp103;

      /* Compute tval1 = t.  */
      tval1 = (int) (t * TWO8);

      x -= __expl_table[T_EXPL_ARG1+2*tval1];
      xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];

      /* Calculate t/32768.  */
      t = x + THREEp96;
      t -= THREEp96;

      /* Compute tval2 = t.  */
      tval2 = (int) (t * TWO15);

      x -= __expl_table[T_EXPL_ARG2+2*tval2];
      xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];

      x = x + xl;

      /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]).  */
      ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
              * __expl_table[T_EXPL_RES2 + tval2];
      n_i = (int)n;
      /* 'unsafe' is 1 iff n_1 != 0.  */
      unsafe = abs(n_i) >= -LDBL_MIN_EXP - 1;
      ex2_u.ieee.exponent += n_i >> unsafe;

      /* Compute scale = 2^n_1.  */
      scale_u.d = 1.0L;
      scale_u.ieee.exponent += n_i - (n_i >> unsafe);

      /* Approximate e^x2 - 1, using a seventh-degree polynomial,
        with maximum error in [-2^-16-2^-53,2^-16+2^-53]
        less than 4.8e-39.  */
      x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));

      /* Return result.  */
      fesetenv (&oldenv);

      result = x22 * ex2_u.d + ex2_u.d;

      /* Now we can test whether the result is ultimate or if we are unsure.
        In the later case we should probably call a mpn based routine to give
        the ultimate result.
        Empirically, this routine is already ultimate in about 99.9986% of
        cases, the test below for the round to nearest case will be false
        in ~ 99.9963% of cases.
        Without proc2 routine maximum error which has been seen is
        0.5000262 ulp.

         union ieee854_long_double ex3_u;

         #ifdef FE_TONEAREST
           fesetround (FE_TONEAREST);
         #endif
         ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
         ex2_u.d = result;
         ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
                             - ex2_u.ieee.exponent;
         n_i = abs (ex3_u.d);
         n_i = (n_i + 1) / 2;
         fesetenv (&oldenv);
         #ifdef FE_TONEAREST
         if (fegetround () == FE_TONEAREST)
           n_i -= 0x4000;
         #endif
         if (!n_i) {
           return __ieee754_expl_proc2 (origx);
         }
       */
      if (!unsafe)
       return result;
      else
       return result * scale_u.d;
    }
  /* Exceptional cases:  */
  else if (isless (x, himark))
    {
      if (__isinfl (x))
       /* e^-inf == 0, with no error.  */
       return 0;
      else
       /* Underflow */
       return TINY * TINY;
    }
  else
    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
    return TWO16383*x;
}

Here is the call graph for this function:


Variable Documentation

const long double C[] [static]

Definition at line 72 of file e_expl.c.