libgcc-math/dbl-64/e_exp.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001 Free Software Foundation
   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, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  19  */
  20 /***************************************************************************/
  21 /*  MODULE_NAME:uexp.c                                                     */
  22 /*                                                                         */
  23 /*  FUNCTION:uexp                                                          */
  24 /*           exp1                                                          */
  25 /*                                                                         */
  26 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
  27 /*              mpa.c mpexp.x slowexp.c                                    */
  28 /*                                                                         */
  29 /* An ultimate exp routine. Given an IEEE double machine number x          */
  30 /* it computes the correctly rounded (to nearest) value of e^x             */
  31 /* Assumption: Machine arithmetic operations are performed in              */
  32 /* round to nearest mode of IEEE 754 standard.                             */
  33 /*                                                                         */
  34 /***************************************************************************/
  35
  36 #include "endian.h"
  37 #include "uexp.h"
  38 #include "mydefs.h"
  39 #include "MathLib.h"
  40 #include "uexp.tbl"
  41 #include "math_private.h"
  42
  43 double __slowexp(double);
  44
  45 /***************************************************************************/
  46 /* An ultimate exp routine. Given an IEEE double machine number x          */
  47 /* it computes the correctly rounded (to nearest) value of e^x             */
  48 /***************************************************************************/
  49 double __ieee754_exp(double x) {
  50   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
  51   mynumber junk1, junk2, binexp  = {{0,0}};
  52 #if 0
  53   int4 k;
  54 #endif
  55   int4 i,j,m,n,ex;
  56
  57   junk1.x = x;
  58   m = junk1.i[HIGH_HALF];
  59   n = m&hugeint;
  60
  61   if (n > smallint && n < bigint) {
  62
  63     y = x*log2e.x + three51.x;
  64     bexp = y - three51.x;      /*  multiply the result by 2**bexp        */
  65
  66     junk1.x = y;
  67
  68     eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
  69     t = x - bexp*ln_two1.x;
  70
  71     y = t + three33.x;
  72     base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
  73     junk2.x = y;
  74     del = (t - base) - eps;    /*  x = bexp*ln(2) + base + del           */
  75     eps = del + del*del*(p3.x*del + p2.x);
  76
  77     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
  78
  79     i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
  80     j = (junk2.i[LOW_HALF]&511)<<1;
  81
  82     al = coar.x[i]*fine.x[j];
  83     bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
  84
  85     rem=(bet + bet*eps)+al*eps;
  86     res = al + rem;
  87     cor = (al - res) + rem;
  88     if  (res == (res+cor*err_0)) return res*binexp.x;
  89     else return __slowexp(x); /*if error is over bound */
  90   }
  91
  92   if (n <= smallint) return 1.0;
  93
  94   if (n >= badint) {
  95     if (n > infint) return(x+x);               /* x is NaN */
  96     if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
  97     /* x is finite,  cause either overflow or underflow  */
  98     if (junk1.i[LOW_HALF] != 0)  return (x+x);                /*  x is NaN  */
  99     return ((x>0)?inf.x:zero );             /* |x| = inf;  return either inf or 0 */
 100   }
 101
 102   y = x*log2e.x + three51.x;
 103   bexp = y - three51.x;
 104   junk1.x = y;
 105   eps = bexp*ln_two2.x;
 106   t = x - bexp*ln_two1.x;
 107   y = t + three33.x;
 108   base = y - three33.x;
 109   junk2.x = y;
 110   del = (t - base) - eps;
 111   eps = del + del*del*(p3.x*del + p2.x);
 112   i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
 113   j = (junk2.i[LOW_HALF]&511)<<1;
 114   al = coar.x[i]*fine.x[j];
 115   bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
 116   rem=(bet + bet*eps)+al*eps;
 117   res = al + rem;
 118   cor = (al - res) + rem;
 119   if (m>>31) {
 120     ex=junk1.i[LOW_HALF];
 121     if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
 122     if (ex >=-1022) {
 123       binexp.i[HIGH_HALF] = (1023+ex)<<20;
 124       if  (res == (res+cor*err_0)) return res*binexp.x;
 125       else return __slowexp(x); /*if error is over bound */
 126     }
 127     ex = -(1022+ex);
 128     binexp.i[HIGH_HALF] = (1023-ex)<<20;
 129     res*=binexp.x;
 130     cor*=binexp.x;
 131     eps=1.0000000001+err_0*binexp.x;
 132     t=1.0+res;
 133     y = ((1.0-t)+res)+cor;
 134     res=t+y;
 135     cor = (t-res)+y;
 136     if (res == (res + eps*cor))
 137     { binexp.i[HIGH_HALF] = 0x00100000;
 138       return (res-1.0)*binexp.x;
 139     }
 140     else return __slowexp(x); /*   if error is over bound    */
 141   }
 142   else {
 143     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
 144     if  (res == (res+cor*err_0)) return res*binexp.x*t256.x;
 145     else return __slowexp(x);
 146   }
 147 }
 148
 149 /************************************************************************/
 150 /* Compute e^(x+xx)(Double-Length number) .The routine also receive     */
 151 /* bound of error of previous calculation .If after computing exp       */
 152 /* error bigger than allows routine return non positive number          */
 153 /*else return   e^(x + xx)   (always positive )                         */
 154 /************************************************************************/
 155
 156 double __exp1(double x, double xx, double error) {
 157   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
 158   mynumber junk1, junk2, binexp  = {{0,0}};
 159 #if 0
 160   int4 k;
 161 #endif
 162   int4 i,j,m,n,ex;
 163
 164   junk1.x = x;
 165   m = junk1.i[HIGH_HALF];
 166   n = m&hugeint;                 /* no sign */
 167
 168   if (n > smallint && n < bigint) {
 169     y = x*log2e.x + three51.x;
 170     bexp = y - three51.x;      /*  multiply the result by 2**bexp        */
 171
 172     junk1.x = y;
 173
 174     eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
 175     t = x - bexp*ln_two1.x;
 176
 177     y = t + three33.x;
 178     base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
 179     junk2.x = y;
 180     del = (t - base) + (xx-eps);    /*  x = bexp*ln(2) + base + del      */
 181     eps = del + del*del*(p3.x*del + p2.x);
 182
 183     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
 184
 185     i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
 186     j = (junk2.i[LOW_HALF]&511)<<1;
 187
 188     al = coar.x[i]*fine.x[j];
 189     bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
 190
 191     rem=(bet + bet*eps)+al*eps;
 192     res = al + rem;
 193     cor = (al - res) + rem;
 194     if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
 195     else return -10.0;
 196   }
 197
 198   if (n <= smallint) return 1.0; /*  if x->0 e^x=1 */
 199
 200   if (n >= badint) {
 201     if (n > infint) return(zero/zero);    /* x is NaN,  return invalid */
 202     if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
 203     /* x is finite,  cause either overflow or underflow  */
 204     if (junk1.i[LOW_HALF] != 0)  return (zero/zero);        /*  x is NaN  */
 205     return ((x>0)?inf.x:zero );   /* |x| = inf;  return either inf or 0 */
 206   }
 207
 208   y = x*log2e.x + three51.x;
 209   bexp = y - three51.x;
 210   junk1.x = y;
 211   eps = bexp*ln_two2.x;
 212   t = x - bexp*ln_two1.x;
 213   y = t + three33.x;
 214   base = y - three33.x;
 215   junk2.x = y;
 216   del = (t - base) + (xx-eps);
 217   eps = del + del*del*(p3.x*del + p2.x);
 218   i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
 219   j = (junk2.i[LOW_HALF]&511)<<1;
 220   al = coar.x[i]*fine.x[j];
 221   bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
 222   rem=(bet + bet*eps)+al*eps;
 223   res = al + rem;
 224   cor = (al - res) + rem;
 225   if (m>>31) {
 226     ex=junk1.i[LOW_HALF];
 227     if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
 228     if (ex >=-1022) {
 229       binexp.i[HIGH_HALF] = (1023+ex)<<20;
 230       if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
 231       else return -10.0;
 232     }
 233     ex = -(1022+ex);
 234     binexp.i[HIGH_HALF] = (1023-ex)<<20;
 235     res*=binexp.x;
 236     cor*=binexp.x;
 237     eps=1.00000000001+(error+err_1)*binexp.x;
 238     t=1.0+res;
 239     y = ((1.0-t)+res)+cor;
 240     res=t+y;
 241     cor = (t-res)+y;
 242     if (res == (res + eps*cor))
 243       {binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
 244     else return -10.0;
 245   }
 246   else {
 247     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
 248     if  (res == (res+cor*(1.0+error+err_1)))
 249       return res*binexp.x*t256.x;
 250     else return -10.0;
 251   }
 252 }