Back to index

glibc  2.9
e_exp.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /***************************************************************************/
00021 /*  MODULE_NAME:uexp.c                                                     */
00022 /*                                                                         */
00023 /*  FUNCTION:uexp                                                          */
00024 /*           exp1                                                          */
00025 /*                                                                         */
00026 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
00027 /*              mpa.c mpexp.x slowexp.c                                    */
00028 /*                                                                         */
00029 /* An ultimate exp routine. Given an IEEE double machine number x          */
00030 /* it computes the correctly rounded (to nearest) value of e^x             */
00031 /* Assumption: Machine arithmetic operations are performed in              */
00032 /* round to nearest mode of IEEE 754 standard.                             */
00033 /*                                                                         */
00034 /***************************************************************************/
00035 
00036 #include "endian.h"
00037 #include "uexp.h"
00038 #include "mydefs.h"
00039 #include "MathLib.h"
00040 #include "uexp.tbl"
00041 #include "math_private.h"
00042 
00043 double __slowexp(double);
00044 
00045 /***************************************************************************/
00046 /* An ultimate exp routine. Given an IEEE double machine number x          */
00047 /* it computes the correctly rounded (to nearest) value of e^x             */
00048 /***************************************************************************/
00049 double __ieee754_exp(double x) {
00050   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
00051   mynumber junk1, junk2, binexp  = {{0,0}};
00052 #if 0
00053   int4 k;
00054 #endif
00055   int4 i,j,m,n,ex;
00056 
00057   junk1.x = x;
00058   m = junk1.i[HIGH_HALF];
00059   n = m&hugeint;
00060 
00061   if (n > smallint && n < bigint) {
00062 
00063     y = x*log2e.x + three51.x;
00064     bexp = y - three51.x;      /*  multiply the result by 2**bexp        */
00065 
00066     junk1.x = y;
00067 
00068     eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
00069     t = x - bexp*ln_two1.x;
00070 
00071     y = t + three33.x;
00072     base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
00073     junk2.x = y;
00074     del = (t - base) - eps;    /*  x = bexp*ln(2) + base + del           */
00075     eps = del + del*del*(p3.x*del + p2.x);
00076 
00077     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
00078 
00079     i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
00080     j = (junk2.i[LOW_HALF]&511)<<1;
00081 
00082     al = coar.x[i]*fine.x[j];
00083     bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
00084 
00085     rem=(bet + bet*eps)+al*eps;
00086     res = al + rem;
00087     cor = (al - res) + rem;
00088     if  (res == (res+cor*err_0)) return res*binexp.x;
00089     else return __slowexp(x); /*if error is over bound */
00090   }
00091 
00092   if (n <= smallint) return 1.0;
00093 
00094   if (n >= badint) {
00095     if (n > infint) return(x+x);               /* x is NaN */
00096     if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
00097     /* x is finite,  cause either overflow or underflow  */
00098     if (junk1.i[LOW_HALF] != 0)  return (x+x);                /*  x is NaN  */
00099     return ((x>0)?inf.x:zero );             /* |x| = inf;  return either inf or 0 */
00100   }
00101 
00102   y = x*log2e.x + three51.x;
00103   bexp = y - three51.x;
00104   junk1.x = y;
00105   eps = bexp*ln_two2.x;
00106   t = x - bexp*ln_two1.x;
00107   y = t + three33.x;
00108   base = y - three33.x;
00109   junk2.x = y;
00110   del = (t - base) - eps;
00111   eps = del + del*del*(p3.x*del + p2.x);
00112   i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
00113   j = (junk2.i[LOW_HALF]&511)<<1;
00114   al = coar.x[i]*fine.x[j];
00115   bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
00116   rem=(bet + bet*eps)+al*eps;
00117   res = al + rem;
00118   cor = (al - res) + rem;
00119   if (m>>31) {
00120     ex=junk1.i[LOW_HALF];
00121     if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
00122     if (ex >=-1022) {
00123       binexp.i[HIGH_HALF] = (1023+ex)<<20;
00124       if  (res == (res+cor*err_0)) return res*binexp.x;
00125       else return __slowexp(x); /*if error is over bound */
00126     }
00127     ex = -(1022+ex);
00128     binexp.i[HIGH_HALF] = (1023-ex)<<20;
00129     res*=binexp.x;
00130     cor*=binexp.x;
00131     eps=1.0000000001+err_0*binexp.x;
00132     t=1.0+res;
00133     y = ((1.0-t)+res)+cor;
00134     res=t+y;
00135     cor = (t-res)+y;
00136     if (res == (res + eps*cor))
00137     { binexp.i[HIGH_HALF] = 0x00100000;
00138       return (res-1.0)*binexp.x;
00139     }
00140     else return __slowexp(x); /*   if error is over bound    */
00141   }
00142   else {
00143     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
00144     if  (res == (res+cor*err_0)) return res*binexp.x*t256.x;
00145     else return __slowexp(x);
00146   }
00147 }
00148 
00149 /************************************************************************/
00150 /* Compute e^(x+xx)(Double-Length number) .The routine also receive     */
00151 /* bound of error of previous calculation .If after computing exp       */
00152 /* error bigger than allows routine return non positive number          */
00153 /*else return   e^(x + xx)   (always positive )                         */
00154 /************************************************************************/
00155 
00156 double __exp1(double x, double xx, double error) {
00157   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
00158   mynumber junk1, junk2, binexp  = {{0,0}};
00159 #if 0
00160   int4 k;
00161 #endif
00162   int4 i,j,m,n,ex;
00163 
00164   junk1.x = x;
00165   m = junk1.i[HIGH_HALF];
00166   n = m&hugeint;                 /* no sign */
00167 
00168   if (n > smallint && n < bigint) {
00169     y = x*log2e.x + three51.x;
00170     bexp = y - three51.x;      /*  multiply the result by 2**bexp        */
00171 
00172     junk1.x = y;
00173 
00174     eps = bexp*ln_two2.x;      /* x = bexp*ln(2) + t - eps               */
00175     t = x - bexp*ln_two1.x;
00176 
00177     y = t + three33.x;
00178     base = y - three33.x;      /* t rounded to a multiple of 2**-18      */
00179     junk2.x = y;
00180     del = (t - base) + (xx-eps);    /*  x = bexp*ln(2) + base + del      */
00181     eps = del + del*del*(p3.x*del + p2.x);
00182 
00183     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
00184 
00185     i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
00186     j = (junk2.i[LOW_HALF]&511)<<1;
00187 
00188     al = coar.x[i]*fine.x[j];
00189     bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
00190 
00191     rem=(bet + bet*eps)+al*eps;
00192     res = al + rem;
00193     cor = (al - res) + rem;
00194     if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
00195     else return -10.0;
00196   }
00197 
00198   if (n <= smallint) return 1.0; /*  if x->0 e^x=1 */
00199 
00200   if (n >= badint) {
00201     if (n > infint) return(zero/zero);    /* x is NaN,  return invalid */
00202     if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
00203     /* x is finite,  cause either overflow or underflow  */
00204     if (junk1.i[LOW_HALF] != 0)  return (zero/zero);        /*  x is NaN  */
00205     return ((x>0)?inf.x:zero );   /* |x| = inf;  return either inf or 0 */
00206   }
00207 
00208   y = x*log2e.x + three51.x;
00209   bexp = y - three51.x;
00210   junk1.x = y;
00211   eps = bexp*ln_two2.x;
00212   t = x - bexp*ln_two1.x;
00213   y = t + three33.x;
00214   base = y - three33.x;
00215   junk2.x = y;
00216   del = (t - base) + (xx-eps);
00217   eps = del + del*del*(p3.x*del + p2.x);
00218   i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
00219   j = (junk2.i[LOW_HALF]&511)<<1;
00220   al = coar.x[i]*fine.x[j];
00221   bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
00222   rem=(bet + bet*eps)+al*eps;
00223   res = al + rem;
00224   cor = (al - res) + rem;
00225   if (m>>31) {
00226     ex=junk1.i[LOW_HALF];
00227     if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
00228     if (ex >=-1022) {
00229       binexp.i[HIGH_HALF] = (1023+ex)<<20;
00230       if  (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
00231       else return -10.0;
00232     }
00233     ex = -(1022+ex);
00234     binexp.i[HIGH_HALF] = (1023-ex)<<20;
00235     res*=binexp.x;
00236     cor*=binexp.x;
00237     eps=1.00000000001+(error+err_1)*binexp.x;
00238     t=1.0+res;
00239     y = ((1.0-t)+res)+cor;
00240     res=t+y;
00241     cor = (t-res)+y;
00242     if (res == (res + eps*cor))
00243       {binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
00244     else return -10.0;
00245   }
00246   else {
00247     binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
00248     if  (res == (res+cor*(1.0+error+err_1)))
00249       return res*binexp.x*t256.x;
00250     else return -10.0;
00251   }
00252 }