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 | |
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 |
263 | |