libgcc-math/dbl-64/e_log.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 /*                                                                   */
  22 /*      MODULE_NAME:ulog.c                                           */
  23 /*                                                                   */
  24 /*      FUNCTION:ulog                                                */
  25 /*                                                                   */
  26 /*      FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h           */
  27 /*                    mpexp.c mplog.c mpa.c                          */
  28 /*                    ulog.tbl                                       */
  29 /*                                                                   */
  30 /* An ultimate log routine. Given an IEEE double machine number x    */
  31 /* it computes the correctly rounded (to nearest) value of log(x).   */
  32 /* Assumption: Machine arithmetic operations are performed in        */
  33 /* round to nearest mode of IEEE 754 standard.                       */
  34 /*                                                                   */
  35 /*********************************************************************/
  36
  37
  38 #include "endian.h"
  39 #include "dla.h"
  40 #include "mpa.h"
  41 #include "MathLib.h"
  42 #include "math_private.h"
  43
  44 void __mplog(mp_no *, mp_no *, int);
  45
  46 /*********************************************************************/
  47 /* An ultimate log routine. Given an IEEE double machine number x     */
  48 /* it computes the correctly rounded (to nearest) value of log(x).   */
  49 /*********************************************************************/
  50 double __ieee754_log(double x) {
  51 #define M 4
  52   static const int pr[M]={8,10,18,32};
  53   int i,j,n,ux,dx,p;
  54 #if 0
  55   int k;
  56 #endif
  57   double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj,
  58          sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb,
  59          t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww,
  60          a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c;
  61   number num;
  62   mp_no mpx,mpy,mpy1,mpy2,mperr;
  63
  64 #include "ulog.tbl"
  65 #include "ulog.h"
  66
  67   /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */
  68
  69   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
  70   n=0;
  71   if (ux < 0x00100000) {
  72     if (((ux & 0x7fffffff) | dx) == 0)  return MHALF/ZERO; /* return -INF */
  73     if (ux < 0) return (x-x)/ZERO;                         /* return NaN  */
  74     n -= 54;    x *= two54.d;                              /* scale x     */
  75     num.d = x;
  76   }
  77   if (ux >= 0x7ff00000) return x+x;                        /* INF or NaN  */
  78
  79   /* Regular values of x */
  80
  81   w = x-ONE;
  82   if (ABS(w) > U03) { goto case_03; }
  83
  84
  85   /*--- Stage I, the case abs(x-1) < 0.03 */
  86
  87   t8 = MHALF*w;
  88   EMULV(t8,w,a,aa,t1,t2,t3,t4,t5)
  89   EADD(w,a,b,bb)
  90
  91   /* Evaluate polynomial II */
  92   polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+
  93           w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w;
  94   c = (aa+bb)+polII;
  95
  96   /* End stage I, case abs(x-1) < 0.03 */
  97   if ((y=b+(c+b*E2)) == b+(c-b*E2))  return y;
  98
  99   /*--- Stage II, the case abs(x-1) < 0.03 */
 100
 101   a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+
 102             w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d))))))));
 103   EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5)
 104   ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2)
 105   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 106   ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2)
 107   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 108   ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2)
 109   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 110   ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2)
 111   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 112   ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2)
 113   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 114   ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2)
 115   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 116   ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2)
 117   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 118   ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2)
 119   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 120   ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2)
 121   MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 122   MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
 123   ADD2(w,ZERO,    s3,ss3, b, bb,t1,t2)
 124
 125   /* End stage II, case abs(x-1) < 0.03 */
 126   if ((y=b+(bb+b*E4)) == b+(bb-b*E4))  return y;
 127   goto stage_n;
 128
 129   /*--- Stage I, the case abs(x-1) > 0.03 */
 130   case_03:
 131
 132   /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
 133   n += (num.i[HIGH_HALF] >> 20) - 1023;
 134   num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
 135   if (num.d > SQRT_2) { num.d *= HALF;  n++; }
 136   u = num.d;  dbl_n = (double) n;
 137
 138   /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
 139   num.d += h1.d;
 140   i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;
 141
 142   /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
 143   num.d = u*Iu[i].d + h2.d;
 144   j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
 145
 146   /* Compute w=(u-ui*vj)/(ui*vj) */
 147   p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V);
 148   q=u-p0;   r0=Iu[i].d*Iv[j].d;   w=q*r0;
 149
 150   /* Evaluate polynomial I */
 151   polI = w+(a2.d+a3.d*w)*w*w;
 152
 153   /* Add up everything */
 154   nln2a = dbl_n*LN2A;
 155   luai  = Lu[i][0].d;   lubi  = Lu[i][1].d;
 156   lvaj  = Lv[j][0].d;   lvbj  = Lv[j][1].d;
 157   EADD(luai,lvaj,sij,ssij)
 158   EADD(nln2a,sij,A  ,ttij)
 159   B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B;
 160   B  = polI+B0;
 161
 162   /* End stage I, case abs(x-1) >= 0.03 */
 163   if ((y=A+(B+E1)) == A+(B-E1))  return y;
 164
 165
 166   /*--- Stage II, the case abs(x-1) > 0.03 */
 167
 168   /* Improve the accuracy of r0 */
 169   EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5)
 170   t=r0*((ONE-sa)-sb);
 171   EADD(r0,t,ra,rb)
 172
 173   /* Compute w */
 174   MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8)
 175
 176   EADD(A,B0,a0,aa0)
 177
 178   /* Evaluate polynomial III */
 179   s1 = (c3.d+(c4.d+c5.d*w)*w)*w;
 180   EADD(c2.d,s1,s2,ss2)
 181   MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8)
 182   MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 183   ADD2(s2,ss2,w,ww,s3,ss3,t1,t2)
 184   ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2)
 185
 186   /* End stage II, case abs(x-1) >= 0.03 */
 187   if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y;
 188
 189
 190   /* Final stages. Use multi-precision arithmetic. */
 191   stage_n:
 192
 193   for (i=0; i<M; i++) {
 194     p = pr[i];
 195     __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
 196     __mplog(&mpx,&mpy,p);
 197     __dbl_mp(e[i].d,&mperr,p);
 198     __add(&mpy,&mperr,&mpy1,p);  __sub(&mpy,&mperr,&mpy2,p);
 199     __mp_dbl(&mpy1,&y1,p);       __mp_dbl(&mpy2,&y2,p);
 200     if (y1==y2)   return y1;
 201   }
 202   return y1;
 203 }