sysdeps/ieee754/dbl-64/e_log.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2016 Free Software Foundation, Inc.
   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.h>
  42 #include <math_private.h>
  43 #include <stap-probe.h>
  44
  45 #ifndef SECTION
  46 # define SECTION
  47 #endif
  48
  49 void __mplog (mp_no *, mp_no *, int);
  50
  51 /*********************************************************************/
  52 /* An ultimate log routine. Given an IEEE double machine number x     */
  53 /* it computes the correctly rounded (to nearest) value of log(x).   */
  54 /*********************************************************************/
  55 double
  56 SECTION
  57 __ieee754_log (double x)
  58 {
  59 #define M 4
  60   static const int pr[M] = { 8, 10, 18, 32 };
  61   int i, j, n, ux, dx, p;
  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;
  78   ux = num.i[HIGH_HALF];
  79   dx = num.i[LOW_HALF];
  80   n = 0;
  81   if (__glibc_unlikely (ux < 0x00100000))
  82     {
  83       if (__glibc_unlikely (((ux & 0x7fffffff) | dx) == 0))
  84         return MHALF / 0.0;     /* return -INF */
  85       if (__glibc_unlikely (ux < 0))
  86         return (x - x) / 0.0;   /* return NaN  */
  87       n -= 54;
  88       x *= two54.d;             /* scale x     */
  89       num.d = x;
  90     }
  91   if (__glibc_unlikely (ux >= 0x7ff00000))
  92     return x + x;               /* INF or NaN  */
  93
  94   /* Regular values of x */
  95
  96   w = x - 1;
  97   if (__glibc_likely (fabs (w) > U03))
  98     goto case_03;
  99
 100   /* log (1) is +0 in all rounding modes.  */
 101   if (w == 0.0)
 102     return 0.0;
 103
 104   /*--- Stage I, the case abs(x-1) < 0.03 */
 105
 106   t8 = MHALF * w;
 107   EMULV (t8, w, a, aa, t1, t2, t3, t4, t5);
 108   EADD (w, a, b, bb);
 109   /* Evaluate polynomial II */
 110   polII = b7.d + w * b8.d;
 111   polII = b6.d + w * polII;
 112   polII = b5.d + w * polII;
 113   polII = b4.d + w * polII;
 114   polII = b3.d + w * polII;
 115   polII = b2.d + w * polII;
 116   polII = b1.d + w * polII;
 117   polII = b0.d + w * polII;
 118   polII *= w * w * w;
 119   c = (aa + bb) + polII;
 120
 121   /* End stage I, case abs(x-1) < 0.03 */
 122   if ((y = b + (c + b * E2)) == b + (c - b * E2))
 123     return y;
 124
 125   /*--- Stage II, the case abs(x-1) < 0.03 */
 126
 127   a = d19.d + w * d20.d;
 128   a = d18.d + w * a;
 129   a = d17.d + w * a;
 130   a = d16.d + w * a;
 131   a = d15.d + w * a;
 132   a = d14.d + w * a;
 133   a = d13.d + w * a;
 134   a = d12.d + w * a;
 135   a = d11.d + w * a;
 136
 137   EMULV (w, a, s2, ss2, t1, t2, t3, t4, t5);
 138   ADD2 (d10.d, dd10.d, s2, ss2, s3, ss3, t1, t2);
 139   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 140   ADD2 (d9.d, dd9.d, s2, ss2, s3, ss3, t1, t2);
 141   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 142   ADD2 (d8.d, dd8.d, s2, ss2, s3, ss3, t1, t2);
 143   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 144   ADD2 (d7.d, dd7.d, s2, ss2, s3, ss3, t1, t2);
 145   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 146   ADD2 (d6.d, dd6.d, s2, ss2, s3, ss3, t1, t2);
 147   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 148   ADD2 (d5.d, dd5.d, s2, ss2, s3, ss3, t1, t2);
 149   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 150   ADD2 (d4.d, dd4.d, s2, ss2, s3, ss3, t1, t2);
 151   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 152   ADD2 (d3.d, dd3.d, s2, ss2, s3, ss3, t1, t2);
 153   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 154   ADD2 (d2.d, dd2.d, s2, ss2, s3, ss3, t1, t2);
 155   MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 156   MUL2 (w, 0, s2, ss2, s3, ss3, t1, t2, t3, t4, t5, t6, t7, t8);
 157   ADD2 (w, 0, s3, ss3, b, bb, t1, t2);
 158
 159   /* End stage II, case abs(x-1) < 0.03 */
 160   if ((y = b + (bb + b * E4)) == b + (bb - b * E4))
 161     return y;
 162   goto stage_n;
 163
 164   /*--- Stage I, the case abs(x-1) > 0.03 */
 165 case_03:
 166
 167   /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
 168   n += (num.i[HIGH_HALF] >> 20) - 1023;
 169   num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000;
 170   if (num.d > SQRT_2)
 171     {
 172       num.d *= HALF;
 173       n++;
 174     }
 175   u = num.d;
 176   dbl_n = (double) n;
 177
 178   /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */
 179   num.d += h1.d;
 180   i = (num.i[HIGH_HALF] & 0x000fffff) >> 12;
 181
 182   /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */
 183   num.d = u * Iu[i].d + h2.d;
 184   j = (num.i[HIGH_HALF] & 0x000fffff) >> 4;
 185
 186   /* Compute w=(u-ui*vj)/(ui*vj) */
 187   p0 = (1 + (i - 75) * DEL_U) * (1 + (j - 180) * DEL_V);
 188   q = u - p0;
 189   r0 = Iu[i].d * Iv[j].d;
 190   w = q * r0;
 191
 192   /* Evaluate polynomial I */
 193   polI = w + (a2.d + a3.d * w) * w * w;
 194
 195   /* Add up everything */
 196   nln2a = dbl_n * LN2A;
 197   luai = Lu[i][0].d;
 198   lubi = Lu[i][1].d;
 199   lvaj = Lv[j][0].d;
 200   lvbj = Lv[j][1].d;
 201   EADD (luai, lvaj, sij, ssij);
 202   EADD (nln2a, sij, A, ttij);
 203   B0 = (((lubi + lvbj) + ssij) + ttij) + dbl_n * LN2B;
 204   B = polI + B0;
 205
 206   /* End stage I, case abs(x-1) >= 0.03 */
 207   if ((y = A + (B + E1)) == A + (B - E1))
 208     return y;
 209
 210
 211   /*--- Stage II, the case abs(x-1) > 0.03 */
 212
 213   /* Improve the accuracy of r0 */
 214   EMULV (p0, r0, sa, sb, t1, t2, t3, t4, t5);
 215   t = r0 * ((1 - sa) - sb);
 216   EADD (r0, t, ra, rb);
 217
 218   /* Compute w */
 219   MUL2 (q, 0, ra, rb, w, ww, t1, t2, t3, t4, t5, t6, t7, t8);
 220
 221   EADD (A, B0, a0, aa0);
 222
 223   /* Evaluate polynomial III */
 224   s1 = (c3.d + (c4.d + c5.d * w) * w) * w;
 225   EADD (c2.d, s1, s2, ss2);
 226   MUL2 (s2, ss2, w, ww, s3, ss3, t1, t2, t3, t4, t5, t6, t7, t8);
 227   MUL2 (s3, ss3, w, ww, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 228   ADD2 (s2, ss2, w, ww, s3, ss3, t1, t2);
 229   ADD2 (s3, ss3, a0, aa0, a1, aa1, t1, t2);
 230
 231   /* End stage II, case abs(x-1) >= 0.03 */
 232   if ((y = a1 + (aa1 + E3)) == a1 + (aa1 - E3))
 233     return y;
 234
 235
 236   /* Final stages. Use multi-precision arithmetic. */
 237 stage_n:
 238
 239   for (i = 0; i < M; i++)
 240     {
 241       p = pr[i];
 242       __dbl_mp (x, &mpx, p);
 243       __dbl_mp (y, &mpy, p);
 244       __mplog (&mpx, &mpy, p);
 245       __dbl_mp (e[i].d, &mperr, p);
 246       __add (&mpy, &mperr, &mpy1, p);
 247       __sub (&mpy, &mperr, &mpy2, p);
 248       __mp_dbl (&mpy1, &y1, p);
 249       __mp_dbl (&mpy2, &y2, p);
 250       if (y1 == y2)
 251         {
 252           LIBC_PROBE (slowlog, 3, &p, &x, &y1);
 253           return y1;
 254         }
 255     }
 256   LIBC_PROBE (slowlog_inexact, 3, &p, &x, &y1);
 257   return y1;
 258 }
 259
 260 #ifndef __ieee754_log
 261 strong_alias (__ieee754_log, __log_finite)
 262 #endif