sysdeps/ieee754/dbl-64/e_log.c

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