Back to index

glibc  2.9
Functions
e_exp.c File Reference
#include "endian.h"
#include "uexp.h"
#include "mydefs.h"
#include "MathLib.h"
#include "uexp.tbl"
#include "math_private.h"

Go to the source code of this file.

Functions

double __slowexp (double)
double __ieee754_exp (double x)
double __exp1 (double x, double xx, double error)

Function Documentation

double __exp1 ( double  x,
double  xx,
double  error 
)

Definition at line 156 of file e_exp.c.

                                                 {
  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  mynumber junk1, junk2, binexp  = {{0,0}};
#if 0
  int4 k;
#endif
  int4 i,j,m,n,ex;

  junk1.x = x;
  m = junk1.i[HIGH_HALF];
  n = m&hugeint;                 /* no sign */

  if (n > smallint && n < bigint) {
    y = x*log2e.x + three51.x;
    bexp = y - three51.x;      /*  multiply the result by 2**bexp        */

    junk1.x = y;

    eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
    t = x - bexp*ln_two1.x;

    y = t + three33.x;
    base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
    junk2.x = y;
    del = (t - base) + (xx-eps);    /*  x = bexp*ln(2) + base + del      */
    eps = del + del*del*(p3.x*del + p2.x);

    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;

    i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
    j = (junk2.i[LOW_HALF]&511)<<1;

    al = coar.x[i]*fine.x[j];
    bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];

    rem=(bet + bet*eps)+al*eps;
    res = al + rem;
    cor = (al - res) + rem;
    if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
    else return -10.0;
  }

  if (n <= smallint) return 1.0; /*  if x->0 e^x=1 */

  if (n >= badint) {
    if (n > infint) return(zero/zero);    /* x is NaN,  return invalid */
    if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
    /* x is finite,  cause either overflow or underflow  */
    if (junk1.i[LOW_HALF] != 0)  return (zero/zero);        /*  x is NaN  */
    return ((x>0)?inf.x:zero );   /* |x| = inf;  return either inf or 0 */
  }

  y = x*log2e.x + three51.x;
  bexp = y - three51.x;
  junk1.x = y;
  eps = bexp*ln_two2.x;
  t = x - bexp*ln_two1.x;
  y = t + three33.x;
  base = y - three33.x;
  junk2.x = y;
  del = (t - base) + (xx-eps);
  eps = del + del*del*(p3.x*del + p2.x);
  i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
  j = (junk2.i[LOW_HALF]&511)<<1;
  al = coar.x[i]*fine.x[j];
  bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
  rem=(bet + bet*eps)+al*eps;
  res = al + rem;
  cor = (al - res) + rem;
  if (m>>31) {
    ex=junk1.i[LOW_HALF];
    if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
    if (ex >=-1022) {
      binexp.i[HIGH_HALF] = (1023+ex)<<20;
      if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
      else return -10.0;
    }
    ex = -(1022+ex);
    binexp.i[HIGH_HALF] = (1023-ex)<<20;
    res*=binexp.x;
    cor*=binexp.x;
    eps=1.00000000001+(error+err_1)*binexp.x;
    t=1.0+res;
    y = ((1.0-t)+res)+cor;
    res=t+y;
    cor = (t-res)+y;
    if (res == (res + eps*cor))
      {binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
    else return -10.0;
  }
  else {
    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
    if  (res == (res+cor*(1.0+error+err_1)))
      return res*binexp.x*t256.x;
    else return -10.0;
  }
}

Here is the caller graph for this function:

double __ieee754_exp ( double  x)

Definition at line 49 of file e_exp.c.

                               {
  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  mynumber junk1, junk2, binexp  = {{0,0}};
#if 0
  int4 k;
#endif
  int4 i,j,m,n,ex;

  junk1.x = x;
  m = junk1.i[HIGH_HALF];
  n = m&hugeint;

  if (n > smallint && n < bigint) {

    y = x*log2e.x + three51.x;
    bexp = y - three51.x;      /*  multiply the result by 2**bexp        */

    junk1.x = y;

    eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
    t = x - bexp*ln_two1.x;

    y = t + three33.x;
    base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
    junk2.x = y;
    del = (t - base) - eps;    /*  x = bexp*ln(2) + base + del           */
    eps = del + del*del*(p3.x*del + p2.x);

    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;

    i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
    j = (junk2.i[LOW_HALF]&511)<<1;

    al = coar.x[i]*fine.x[j];
    bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];

    rem=(bet + bet*eps)+al*eps;
    res = al + rem;
    cor = (al - res) + rem;
    if  (res == (res+cor*err_0)) return res*binexp.x;
    else return __slowexp(x); /*if error is over bound */
  }

  if (n <= smallint) return 1.0;

  if (n >= badint) {
    if (n > infint) return(x+x);               /* x is NaN */
    if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
    /* x is finite,  cause either overflow or underflow  */
    if (junk1.i[LOW_HALF] != 0)  return (x+x);                /*  x is NaN  */
    return ((x>0)?inf.x:zero );             /* |x| = inf;  return either inf or 0 */
  }

  y = x*log2e.x + three51.x;
  bexp = y - three51.x;
  junk1.x = y;
  eps = bexp*ln_two2.x;
  t = x - bexp*ln_two1.x;
  y = t + three33.x;
  base = y - three33.x;
  junk2.x = y;
  del = (t - base) - eps;
  eps = del + del*del*(p3.x*del + p2.x);
  i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
  j = (junk2.i[LOW_HALF]&511)<<1;
  al = coar.x[i]*fine.x[j];
  bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
  rem=(bet + bet*eps)+al*eps;
  res = al + rem;
  cor = (al - res) + rem;
  if (m>>31) {
    ex=junk1.i[LOW_HALF];
    if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
    if (ex >=-1022) {
      binexp.i[HIGH_HALF] = (1023+ex)<<20;
      if  (res == (res+cor*err_0)) return res*binexp.x;
      else return __slowexp(x); /*if error is over bound */
    }
    ex = -(1022+ex);
    binexp.i[HIGH_HALF] = (1023-ex)<<20;
    res*=binexp.x;
    cor*=binexp.x;
    eps=1.0000000001+err_0*binexp.x;
    t=1.0+res;
    y = ((1.0-t)+res)+cor;
    res=t+y;
    cor = (t-res)+y;
    if (res == (res + eps*cor))
    { binexp.i[HIGH_HALF] = 0x00100000;
      return (res-1.0)*binexp.x;
    }
    else return __slowexp(x); /*   if error is over bound    */
  }
  else {
    binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
    if  (res == (res+cor*err_0)) return res*binexp.x*t256.x;
    else return __slowexp(x);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

double __slowexp ( double  )

Definition at line 37 of file slowexp.c.

                           {
  double w,z,res,eps=3.0e-26;
#if 0
  double y;
#endif
  int p;
#if 0
  int orig,i;
#endif
  mp_no mpx, mpy, mpz,mpw,mpeps,mpcor;

  p=6;
  __dbl_mp(x,&mpx,p); /* Convert a double precision number  x               */
                    /* into a multiple precision number mpx with prec. p. */
  __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */
  __dbl_mp(eps,&mpeps,p);
  __mul(&mpeps,&mpy,&mpcor,p);
  __add(&mpy,&mpcor,&mpw,p);
  __sub(&mpy,&mpcor,&mpz,p);
  __mp_dbl(&mpw, &w, p);
  __mp_dbl(&mpz, &z, p);
  if (w == z) return w;
  else  {                   /* if calculating is not exactly   */
    p = 32;
    __dbl_mp(x,&mpx,p);
    __mpexp(&mpx, &mpy, p);
    __mp_dbl(&mpy, &res, p);
    return res;
  }
}

Here is the call graph for this function: