1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2018 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
49void __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/*********************************************************************/
55double
56SECTION
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 */
165case_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. */
237stage_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
261strong_alias (__ieee754_log, __log_finite)
262#endif
263