sysdeps/ieee754/dbl-64/e_exp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2013 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /***************************************************************************/
  20 /*  MODULE_NAME:uexp.c                                                     */
  21 /*                                                                         */
  22 /*  FUNCTION:uexp                                                          */
  23 /*           exp1                                                          */
  24 /*                                                                         */
  25 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
  26 /*              mpa.c mpexp.x slowexp.c                                    */
  27 /*                                                                         */
  28 /* An ultimate exp routine. Given an IEEE double machine number x          */
  29 /* it computes the correctly rounded (to nearest) value of e^x             */
  30 /* Assumption: Machine arithmetic operations are performed in              */
  31 /* round to nearest mode of IEEE 754 standard.                             */
  32 /*                                                                         */
  33 /***************************************************************************/
  34
  35 #include "endian.h"
  36 #include "uexp.h"
  37 #include "mydefs.h"
  38 #include "MathLib.h"
  39 #include "uexp.tbl"
  40 #include <math_private.h>
  41 #include <fenv.h>
  42
  43 #ifndef SECTION
  44 # define SECTION
  45 #endif
  46
  47 double __slowexp(double);
  48
  49 /***************************************************************************/
  50 /* An ultimate exp routine. Given an IEEE double machine number x          */
  51 /* it computes the correctly rounded (to nearest) value of e^x             */
  52 /***************************************************************************/
  53 double
  54 SECTION
  55 __ieee754_exp(double x) {
  56   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  57   mynumber junk1, junk2, binexp  = {{0,0}};
  58 #if 0
  59   int4 k;
  60 #endif
  61   int4 i,j,m,n,ex;
  62   double retval;
  63
  64   SET_RESTORE_ROUND (FE_TONEAREST);
  65
  66   junk1.x = x;
  67   m = junk1.i[HIGH_HALF];
  68   n = m&hugeint;
  69
  70   if (n > smallint && n < bigint) {
  71
  72     y = x*log2e.x + three51.x;
  73     bexp = y - three51.x;      /*  multiply the result by 2**bexp        */
  74
  75     junk1.x = y;
  76
  77     eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
  78     t = x - bexp*ln_two1.x;
  79
  80     y = t + three33.x;
  81     base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
  82     junk2.x = y;
  83     del = (t - base) - eps;    /*  x = bexp*ln(2) + base + del           */
  84     eps = del + del*del*(p3.x*del + p2.x);
  85
  86     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
  87
  88     i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
  89     j = (junk2.i[LOW_HALF]&511)<<1;
  90
  91     al = coar.x[i]*fine.x[j];
  92     bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
  93
  94     rem=(bet + bet*eps)+al*eps;
  95     res = al + rem;
  96     cor = (al - res) + rem;
  97     if  (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; }
  98     else { retval = __slowexp(x); goto ret; } /*if error is over bound */
  99   }
 100
 101   if (n <= smallint) { retval = 1.0; goto ret; }
 102
 103   if (n >= badint) {
 104     if (n > infint) { retval = x+x; goto ret; }               /* x is NaN */
 105     if (n < infint) { retval = (x>0) ? (hhuge*hhuge) : (tiny*tiny); goto ret; }
 106     /* x is finite,  cause either overflow or underflow  */
 107     if (junk1.i[LOW_HALF] != 0) { retval = x+x; goto ret; } /*  x is NaN  */
 108     retval = (x>0)?inf.x:zero;             /* |x| = inf;  return either inf or 0 */
 109     goto ret;
 110   }
 111
 112   y = x*log2e.x + three51.x;
 113   bexp = y - three51.x;
 114   junk1.x = y;
 115   eps = bexp*ln_two2.x;
 116   t = x - bexp*ln_two1.x;
 117   y = t + three33.x;
 118   base = y - three33.x;
 119   junk2.x = y;
 120   del = (t - base) - eps;
 121   eps = del + del*del*(p3.x*del + p2.x);
 122   i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
 123   j = (junk2.i[LOW_HALF]&511)<<1;
 124   al = coar.x[i]*fine.x[j];
 125   bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
 126   rem=(bet + bet*eps)+al*eps;
 127   res = al + rem;
 128   cor = (al - res) + rem;
 129   if (m>>31) {
 130     ex=junk1.i[LOW_HALF];
 131     if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
 132     if (ex >=-1022) {
 133       binexp.i[HIGH_HALF] = (1023+ex)<<20;
 134       if  (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; }
 135       else { retval = __slowexp(x); goto ret; } /*if error is over bound */
 136     }
 137     ex = -(1022+ex);
 138     binexp.i[HIGH_HALF] = (1023-ex)<<20;
 139     res*=binexp.x;
 140     cor*=binexp.x;
 141     eps=1.0000000001+err_0*binexp.x;
 142     t=1.0+res;
 143     y = ((1.0-t)+res)+cor;
 144     res=t+y;
 145     cor = (t-res)+y;
 146     if (res == (res + eps*cor))
 147     { binexp.i[HIGH_HALF] = 0x00100000;
 148       retval = (res-1.0)*binexp.x;
 149       goto ret;
 150     }
 151     else { retval = __slowexp(x); goto ret; } /*   if error is over bound    */
 152   }
 153   else {
 154     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
 155     if (res == (res+cor*err_0)) { retval = res*binexp.x*t256.x; goto ret; }
 156     else { retval = __slowexp(x); goto ret; }
 157   }
 158  ret:
 159   return retval;
 160 }
 161 #ifndef __ieee754_exp
 162 strong_alias (__ieee754_exp, __exp_finite)
 163 #endif
 164
 165 /************************************************************************/
 166 /* Compute e^(x+xx)(Double-Length number) .The routine also receive     */
 167 /* bound of error of previous calculation .If after computing exp       */
 168 /* error bigger than allows routine return non positive number          */
 169 /*else return   e^(x + xx)   (always positive )                         */
 170 /************************************************************************/
 171
 172 double
 173 SECTION
 174 __exp1(double x, double xx, double error) {
 175   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
 176   mynumber junk1, junk2, binexp  = {{0,0}};
 177 #if 0
 178   int4 k;
 179 #endif
 180   int4 i,j,m,n,ex;
 181
 182   junk1.x = x;
 183   m = junk1.i[HIGH_HALF];
 184   n = m&hugeint;                 /* no sign */
 185
 186   if (n > smallint && n < bigint) {
 187     y = x*log2e.x + three51.x;
 188     bexp = y - three51.x;      /*  multiply the result by 2**bexp        */
 189
 190     junk1.x = y;
 191
 192     eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
 193     t = x - bexp*ln_two1.x;
 194
 195     y = t + three33.x;
 196     base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
 197     junk2.x = y;
 198     del = (t - base) + (xx-eps);    /*  x = bexp*ln(2) + base + del      */
 199     eps = del + del*del*(p3.x*del + p2.x);
 200
 201     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
 202
 203     i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
 204     j = (junk2.i[LOW_HALF]&511)<<1;
 205
 206     al = coar.x[i]*fine.x[j];
 207     bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
 208
 209     rem=(bet + bet*eps)+al*eps;
 210     res = al + rem;
 211     cor = (al - res) + rem;
 212     if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
 213     else return -10.0;
 214   }
 215
 216   if (n <= smallint) return 1.0; /*  if x->0 e^x=1 */
 217
 218   if (n >= badint) {
 219     if (n > infint) return(zero/zero);    /* x is NaN,  return invalid */
 220     if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
 221     /* x is finite,  cause either overflow or underflow  */
 222     if (junk1.i[LOW_HALF] != 0)  return (zero/zero);        /*  x is NaN  */
 223     return ((x>0)?inf.x:zero );   /* |x| = inf;  return either inf or 0 */
 224   }
 225
 226   y = x*log2e.x + three51.x;
 227   bexp = y - three51.x;
 228   junk1.x = y;
 229   eps = bexp*ln_two2.x;
 230   t = x - bexp*ln_two1.x;
 231   y = t + three33.x;
 232   base = y - three33.x;
 233   junk2.x = y;
 234   del = (t - base) + (xx-eps);
 235   eps = del + del*del*(p3.x*del + p2.x);
 236   i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
 237   j = (junk2.i[LOW_HALF]&511)<<1;
 238   al = coar.x[i]*fine.x[j];
 239   bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
 240   rem=(bet + bet*eps)+al*eps;
 241   res = al + rem;
 242   cor = (al - res) + rem;
 243   if (m>>31) {
 244     ex=junk1.i[LOW_HALF];
 245     if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
 246     if (ex >=-1022) {
 247       binexp.i[HIGH_HALF] = (1023+ex)<<20;
 248       if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
 249       else return -10.0;
 250     }
 251     ex = -(1022+ex);
 252     binexp.i[HIGH_HALF] = (1023-ex)<<20;
 253     res*=binexp.x;
 254     cor*=binexp.x;
 255     eps=1.00000000001+(error+err_1)*binexp.x;
 256     t=1.0+res;
 257     y = ((1.0-t)+res)+cor;
 258     res=t+y;
 259     cor = (t-res)+y;
 260     if (res == (res + eps*cor))
 261       {binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
 262     else return -10.0;
 263   }
 264   else {
 265     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
 266     if  (res == (res+cor*(1.0+error+err_1)))
 267       return res*binexp.x*t256.x;
 268     else return -10.0;
 269   }
 270 }