sysdeps/ieee754/dbl-64/e_exp.c

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