Back to index

glibc  2.9
e_log.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 /*                                                                   */
00022 /*      MODULE_NAME:ulog.c                                           */
00023 /*                                                                   */
00024 /*      FUNCTION:ulog                                                */
00025 /*                                                                   */
00026 /*      FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h           */
00027 /*                    mpexp.c mplog.c mpa.c                          */
00028 /*                    ulog.tbl                                       */
00029 /*                                                                   */
00030 /* An ultimate log routine. Given an IEEE double machine number x    */
00031 /* it computes the correctly rounded (to nearest) value of log(x).   */
00032 /* Assumption: Machine arithmetic operations are performed in        */
00033 /* round to nearest mode of IEEE 754 standard.                       */
00034 /*                                                                   */
00035 /*********************************************************************/
00036 
00037 
00038 #include "endian.h"
00039 #include "dla.h"
00040 #include "mpa.h"
00041 #include "MathLib.h"
00042 #include "math_private.h"
00043 
00044 void __mplog(mp_no *, mp_no *, int);
00045 
00046 /*********************************************************************/
00047 /* An ultimate log routine. Given an IEEE double machine number x     */
00048 /* it computes the correctly rounded (to nearest) value of log(x).   */
00049 /*********************************************************************/
00050 double __ieee754_log(double x) {
00051 #define M 4
00052   static const int pr[M]={8,10,18,32};
00053   int i,j,n,ux,dx,p;
00054 #if 0
00055   int k;
00056 #endif
00057   double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
00058          sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
00059          t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww,
00060          a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
00061   number num;
00062   mp_no mpx,mpy,mpy1,mpy2,mperr;
00063 
00064 #include "ulog.tbl"
00065 #include "ulog.h"
00066 
00067   /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */
00068 
00069   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
00070   n=0;
00071   if (ux < 0x00100000) {
00072     if (((ux & 0x7fffffff) | dx) == 0)  return MHALF/ZERO; /* return -INF */
00073     if (ux < 0) return (x-x)/ZERO;                         /* return NaN  */
00074     n -= 54;    x *= two54.d;                              /* scale x     */
00075     num.d = x;
00076   }
00077   if (ux >= 0x7ff00000) return x+x;                        /* INF or NaN  */
00078 
00079   /* Regular values of x */
00080 
00081   w = x-ONE;
00082   if (ABS(w) > U03) { goto case_03; }
00083 
00084 
00085   /*--- Stage I, the case abs(x-1) < 0.03 */
00086 
00087   t8 = MHALF*w;
00088   EMULV(t8,w,a,aa,t1,t2,t3,t4,t5)
00089   EADD(w,a,b,bb)
00090 
00091   /* Evaluate polynomial II */
00092   polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+
00093           w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
00094   c = (aa+bb)+polII;
00095 
00096   /* End stage I, case abs(x-1) < 0.03 */
00097   if ((y=b+(c+b*E2)) == b+(c-b*E2))  return y;
00098 
00099   /*--- Stage II, the case abs(x-1) < 0.03 */
00100 
00101   a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+
00102             w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
00103   EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
00104   ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
00105   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00106   ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
00107   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00108   ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
00109   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00110   ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
00111   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00112   ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
00113   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00114   ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
00115   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00116   ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
00117   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00118   ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
00119   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00120   ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
00121   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00122   MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
00123   ADD2(w,ZERO,    s3,ss3, b, bb,t1,t2)
00124 
00125   /* End stage II, case abs(x-1) < 0.03 */
00126   if ((y=b+(bb+b*E4)) == b+(bb-b*E4))  return y;
00127   goto stage_n;
00128 
00129   /*--- Stage I, the case abs(x-1) > 0.03 */
00130   case_03:
00131 
00132   /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
00133   n += (num.i[HIGH_HALF] >> 20) - 1023;
00134   num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
00135   if (num.d > SQRT_2) { num.d *= HALF;  n++; }
00136   u = num.d;  dbl_n = (double) n;
00137 
00138   /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
00139   num.d += h1.d;
00140   i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;
00141 
00142   /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
00143   num.d = u*Iu[i].d + h2.d;
00144   j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
00145 
00146   /* Compute w=(u-ui*vj)/(ui*vj) */
00147   p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
00148   q=u-p0;   r0=Iu[i].d*Iv[j].d;   w=q*r0;
00149 
00150   /* Evaluate polynomial I */
00151   polI = w+(a2.d+a3.d*w)*w*w;
00152 
00153   /* Add up everything */
00154   nln2a = dbl_n*LN2A;
00155   luai  = Lu[i][0].d;   lubi  = Lu[i][1].d;
00156   lvaj  = Lv[j][0].d;   lvbj  = Lv[j][1].d;
00157   EADD(luai,lvaj,sij,ssij)
00158   EADD(nln2a,sij,A  ,ttij)
00159   B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B;
00160   B  = polI+B0;
00161 
00162   /* End stage I, case abs(x-1) >= 0.03 */
00163   if ((y=A+(B+E1)) == A+(B-E1))  return y;
00164 
00165 
00166   /*--- Stage II, the case abs(x-1) > 0.03 */
00167 
00168   /* Improve the accuracy of r0 */
00169   EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
00170   t=r0*((ONE-sa)-sb);
00171   EADD(r0,t,ra,rb)
00172 
00173   /* Compute w */
00174   MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
00175 
00176   EADD(A,B0,a0,aa0)
00177 
00178   /* Evaluate polynomial III */
00179   s1 = (c3.d+(c4.d+c5.d*w)*w)*w;
00180   EADD(c2.d,s1,s2,ss2)
00181   MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
00182   MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00183   ADD2(s2,ss2,w,ww,s3,ss3,t1,t2)
00184   ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2)
00185 
00186   /* End stage II, case abs(x-1) >= 0.03 */
00187   if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y;
00188 
00189 
00190   /* Final stages. Use multi-precision arithmetic. */
00191   stage_n:
00192 
00193   for (i=0; i<M; i++) {
00194     p = pr[i];
00195     __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
00196     __mplog(&mpx,&mpy,p);
00197     __dbl_mp(e[i].d,&mperr,p);
00198     __add(&mpy,&mperr,&mpy1,p);  __sub(&mpy,&mperr,&mpy2,p);
00199     __mp_dbl(&mpy1,&y1,p);       __mp_dbl(&mpy2,&y2,p);
00200     if (y1==y2)   return y1;
00201   }
00202   return y1;
00203 }