glibc/sysdeps/ieee754/dbl-64/e_exp.c

/*
 * IBM Accurate Mathematical Library
 * written by International Business Machines Corp.
 * Copyright (C) 2001 Free Software Foundation
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
/***************************************************************************/
/*  MODULE_NAME:uexp.c                                                     */
/*                                                                         */
/*  FUNCTION:uexp                                                          */
/*           exp1                                                          */
/*                                                                         */
/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
/*              mpa.c mpexp.x slowexp.c                                    */
/*                                                                         */
/* An ultimate exp routine. Given an IEEE double machine number x          */
/* it computes the correctly rounded (to nearest) value of e^x             */
/* Assumption: Machine arithmetic operations are performed in              */
/* round to nearest mode of IEEE 754 standard.                             */
/*                                                                         */
/***************************************************************************/

#include "endian.h"
#include "uexp.h"
#include "mydefs.h"
#include "MathLib.h"
#include "uexp.tbl"
#include "math_private.h"

double __slowexp(double);

/***************************************************************************/
/* An ultimate exp routine. Given an IEEE double machine number x          */
/* it computes the correctly rounded (to nearest) value of e^x             */
/***************************************************************************/
double __ieee754_exp(double x) {
  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(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) - 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);
  }
}

/************************************************************************/
/* Compute e^(x+xx)(Double-Length number) .The routine also receive     */
/* bound of error of previous calculation .If after computing exp       */
/* error bigger than allows routine return non positive number          */
/*else return   e^(x + xx)   (always positive )                         */
/************************************************************************/

double __exp1(double x, double xx, double error) {
  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;
  }
}