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 <sysdeps/ieee754/ldbl-128/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 TWO1023   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 TWO1023   C[9]
#define TWO15   C[11]
#define TWO8   C[10]

Function Documentation

long double __ieee754_expl ( long double  x)

Definition at line 136 of file e_expl.c.

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

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

      n = __roundl (x*M_1_LN2);
      x = x-n*M_LN2_0;
      xl = n*M_LN2_1;

      tval1 = __roundl (x*TWO8);
      x -= __expl_table[T_EXPL_ARG1+2*tval1];
      xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];

      tval2 = __roundl (x*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 = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
      ex2_u.ieee.exponent += n_i >> unsafe;
      /* Fortunately, there are no subnormal lowpart doubles in
        __expl_table, only normal values and zeros.
        But after scaling it can be subnormal.  */
      exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe);
      if (ex2_u.ieee.exponent2 == 0)
       /* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */;
      else if (exponent2 > 0)
       ex2_u.ieee.exponent2 = exponent2;
      else if (exponent2 <= -54)
       {
         ex2_u.ieee.exponent2 = 0;
         ex2_u.ieee.mantissa2 = 0;
         ex2_u.ieee.mantissa3 = 0;
       }
      else
       {
         static const double
           two54 = 1.80143985094819840000e+16, /* 4350000000000000 */
           twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */
         ex2_u.dd[1] *= two54;
         ex2_u.ieee.exponent2 += n_i >> unsafe;
         ex2_u.dd[1] *= twom54;
       }

      /* 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 TWO1023*x;
}

Here is the call graph for this function:


Variable Documentation

const long double C[] [static]

Definition at line 71 of file e_expl.c.