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 | /* MODULE_NAME:uexp.c */ |
21 | /* */ |
22 | /* FUNCTION:uexp */ |
23 | /* exp1 */ |
24 | /* */ |
25 | /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */ |
26 | /* mpa.c mpexp.x slowexp.c */ |
27 | /* */ |
28 | /* An ultimate exp routine. Given an IEEE double machine number x */ |
29 | /* it computes the correctly rounded (to nearest) value of e^x */ |
30 | /* Assumption: Machine arithmetic operations are performed in */ |
31 | /* round to nearest mode of IEEE 754 standard. */ |
32 | /* */ |
33 | /***************************************************************************/ |
34 | |
35 | #include <math.h> |
36 | #include "endian.h" |
37 | #include "uexp.h" |
38 | #include "mydefs.h" |
39 | #include "MathLib.h" |
40 | #include "uexp.tbl" |
41 | #include <math_private.h> |
42 | #include <fenv.h> |
43 | #include <float.h> |
44 | |
45 | #ifndef SECTION |
46 | # define SECTION |
47 | #endif |
48 | |
49 | double __slowexp (double); |
50 | |
51 | /* An ultimate exp routine. Given an IEEE double machine number x it computes |
52 | the correctly rounded (to nearest) value of e^x. */ |
53 | double |
54 | SECTION |
55 | __ieee754_exp (double x) |
56 | { |
57 | double bexp, t, eps, del, base, y, al, bet, res, rem, cor; |
58 | mynumber junk1, junk2, binexp = {{0, 0}}; |
59 | int4 i, j, m, n, ex; |
60 | double retval; |
61 | |
62 | { |
63 | SET_RESTORE_ROUND (FE_TONEAREST); |
64 | |
65 | junk1.x = x; |
66 | m = junk1.i[HIGH_HALF]; |
67 | n = m & hugeint; |
68 | |
69 | if (n > smallint && n < bigint) |
70 | { |
71 | y = x * log2e.x + three51.x; |
72 | bexp = y - three51.x; /* multiply the result by 2**bexp */ |
73 | |
74 | junk1.x = y; |
75 | |
76 | eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */ |
77 | t = x - bexp * ln_two1.x; |
78 | |
79 | y = t + three33.x; |
80 | base = y - three33.x; /* t rounded to a multiple of 2**-18 */ |
81 | junk2.x = y; |
82 | del = (t - base) - eps; /* x = bexp*ln(2) + base + del */ |
83 | eps = del + del * del * (p3.x * del + p2.x); |
84 | |
85 | binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20; |
86 | |
87 | i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; |
88 | j = (junk2.i[LOW_HALF] & 511) << 1; |
89 | |
90 | al = coar.x[i] * fine.x[j]; |
91 | bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) |
92 | + coar.x[i + 1] * fine.x[j + 1]); |
93 | |
94 | rem = (bet + bet * eps) + al * eps; |
95 | res = al + rem; |
96 | cor = (al - res) + rem; |
97 | if (res == (res + cor * err_0)) |
98 | { |
99 | retval = res * binexp.x; |
100 | goto ret; |
101 | } |
102 | else |
103 | { |
104 | retval = __slowexp (x); |
105 | goto ret; |
106 | } /*if error is over bound */ |
107 | } |
108 | |
109 | if (n <= smallint) |
110 | { |
111 | retval = 1.0; |
112 | goto ret; |
113 | } |
114 | |
115 | if (n >= badint) |
116 | { |
117 | if (n > infint) |
118 | { |
119 | retval = x + x; |
120 | goto ret; |
121 | } /* x is NaN */ |
122 | if (n < infint) |
123 | { |
124 | if (x > 0) |
125 | goto ret_huge; |
126 | else |
127 | goto ret_tiny; |
128 | } |
129 | /* x is finite, cause either overflow or underflow */ |
130 | if (junk1.i[LOW_HALF] != 0) |
131 | { |
132 | retval = x + x; |
133 | goto ret; |
134 | } /* x is NaN */ |
135 | retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */ |
136 | goto ret; |
137 | } |
138 | |
139 | y = x * log2e.x + three51.x; |
140 | bexp = y - three51.x; |
141 | junk1.x = y; |
142 | eps = bexp * ln_two2.x; |
143 | t = x - bexp * ln_two1.x; |
144 | y = t + three33.x; |
145 | base = y - three33.x; |
146 | junk2.x = y; |
147 | del = (t - base) - eps; |
148 | eps = del + del * del * (p3.x * del + p2.x); |
149 | i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; |
150 | j = (junk2.i[LOW_HALF] & 511) << 1; |
151 | al = coar.x[i] * fine.x[j]; |
152 | bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) |
153 | + coar.x[i + 1] * fine.x[j + 1]); |
154 | rem = (bet + bet * eps) + al * eps; |
155 | res = al + rem; |
156 | cor = (al - res) + rem; |
157 | if (m >> 31) |
158 | { |
159 | ex = junk1.i[LOW_HALF]; |
160 | if (res < 1.0) |
161 | { |
162 | res += res; |
163 | cor += cor; |
164 | ex -= 1; |
165 | } |
166 | if (ex >= -1022) |
167 | { |
168 | binexp.i[HIGH_HALF] = (1023 + ex) << 20; |
169 | if (res == (res + cor * err_0)) |
170 | { |
171 | retval = res * binexp.x; |
172 | goto ret; |
173 | } |
174 | else |
175 | { |
176 | retval = __slowexp (x); |
177 | goto check_uflow_ret; |
178 | } /*if error is over bound */ |
179 | } |
180 | ex = -(1022 + ex); |
181 | binexp.i[HIGH_HALF] = (1023 - ex) << 20; |
182 | res *= binexp.x; |
183 | cor *= binexp.x; |
184 | eps = 1.0000000001 + err_0 * binexp.x; |
185 | t = 1.0 + res; |
186 | y = ((1.0 - t) + res) + cor; |
187 | res = t + y; |
188 | cor = (t - res) + y; |
189 | if (res == (res + eps * cor)) |
190 | { |
191 | binexp.i[HIGH_HALF] = 0x00100000; |
192 | retval = (res - 1.0) * binexp.x; |
193 | goto check_uflow_ret; |
194 | } |
195 | else |
196 | { |
197 | retval = __slowexp (x); |
198 | goto check_uflow_ret; |
199 | } /* if error is over bound */ |
200 | check_uflow_ret: |
201 | if (retval < DBL_MIN) |
202 | { |
203 | double force_underflow = tiny * tiny; |
204 | math_force_eval (force_underflow); |
205 | } |
206 | if (retval == 0) |
207 | goto ret_tiny; |
208 | goto ret; |
209 | } |
210 | else |
211 | { |
212 | binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20; |
213 | if (res == (res + cor * err_0)) |
214 | retval = res * binexp.x * t256.x; |
215 | else |
216 | retval = __slowexp (x); |
217 | if (isinf (retval)) |
218 | goto ret_huge; |
219 | else |
220 | goto ret; |
221 | } |
222 | } |
223 | ret: |
224 | return retval; |
225 | |
226 | ret_huge: |
227 | return hhuge * hhuge; |
228 | |
229 | ret_tiny: |
230 | return tiny * tiny; |
231 | } |
232 | #ifndef __ieee754_exp |
233 | strong_alias (__ieee754_exp, __exp_finite) |
234 | #endif |
235 | |
236 | /* Compute e^(x+xx). The routine also receives bound of error of previous |
237 | calculation. If after computing exp the error exceeds the allowed bounds, |
238 | the routine returns a non-positive number. Otherwise it returns the |
239 | computed result, which is always positive. */ |
240 | double |
241 | SECTION |
242 | __exp1 (double x, double xx, double error) |
243 | { |
244 | double bexp, t, eps, del, base, y, al, bet, res, rem, cor; |
245 | mynumber junk1, junk2, binexp = {{0, 0}}; |
246 | int4 i, j, m, n, ex; |
247 | |
248 | junk1.x = x; |
249 | m = junk1.i[HIGH_HALF]; |
250 | n = m & hugeint; /* no sign */ |
251 | |
252 | if (n > smallint && n < bigint) |
253 | { |
254 | y = x * log2e.x + three51.x; |
255 | bexp = y - three51.x; /* multiply the result by 2**bexp */ |
256 | |
257 | junk1.x = y; |
258 | |
259 | eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */ |
260 | t = x - bexp * ln_two1.x; |
261 | |
262 | y = t + three33.x; |
263 | base = y - three33.x; /* t rounded to a multiple of 2**-18 */ |
264 | junk2.x = y; |
265 | del = (t - base) + (xx - eps); /* x = bexp*ln(2) + base + del */ |
266 | eps = del + del * del * (p3.x * del + p2.x); |
267 | |
268 | binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20; |
269 | |
270 | i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; |
271 | j = (junk2.i[LOW_HALF] & 511) << 1; |
272 | |
273 | al = coar.x[i] * fine.x[j]; |
274 | bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) |
275 | + coar.x[i + 1] * fine.x[j + 1]); |
276 | |
277 | rem = (bet + bet * eps) + al * eps; |
278 | res = al + rem; |
279 | cor = (al - res) + rem; |
280 | if (res == (res + cor * (1.0 + error + err_1))) |
281 | return res * binexp.x; |
282 | else |
283 | return -10.0; |
284 | } |
285 | |
286 | if (n <= smallint) |
287 | return 1.0; /* if x->0 e^x=1 */ |
288 | |
289 | if (n >= badint) |
290 | { |
291 | if (n > infint) |
292 | return (zero / zero); /* x is NaN, return invalid */ |
293 | if (n < infint) |
294 | return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny)); |
295 | /* x is finite, cause either overflow or underflow */ |
296 | if (junk1.i[LOW_HALF] != 0) |
297 | return (zero / zero); /* x is NaN */ |
298 | return ((x > 0) ? inf.x : zero); /* |x| = inf; return either inf or 0 */ |
299 | } |
300 | |
301 | y = x * log2e.x + three51.x; |
302 | bexp = y - three51.x; |
303 | junk1.x = y; |
304 | eps = bexp * ln_two2.x; |
305 | t = x - bexp * ln_two1.x; |
306 | y = t + three33.x; |
307 | base = y - three33.x; |
308 | junk2.x = y; |
309 | del = (t - base) + (xx - eps); |
310 | eps = del + del * del * (p3.x * del + p2.x); |
311 | i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356; |
312 | j = (junk2.i[LOW_HALF] & 511) << 1; |
313 | al = coar.x[i] * fine.x[j]; |
314 | bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j]) |
315 | + coar.x[i + 1] * fine.x[j + 1]); |
316 | rem = (bet + bet * eps) + al * eps; |
317 | res = al + rem; |
318 | cor = (al - res) + rem; |
319 | if (m >> 31) |
320 | { |
321 | ex = junk1.i[LOW_HALF]; |
322 | if (res < 1.0) |
323 | { |
324 | res += res; |
325 | cor += cor; |
326 | ex -= 1; |
327 | } |
328 | if (ex >= -1022) |
329 | { |
330 | binexp.i[HIGH_HALF] = (1023 + ex) << 20; |
331 | if (res == (res + cor * (1.0 + error + err_1))) |
332 | return res * binexp.x; |
333 | else |
334 | return -10.0; |
335 | } |
336 | ex = -(1022 + ex); |
337 | binexp.i[HIGH_HALF] = (1023 - ex) << 20; |
338 | res *= binexp.x; |
339 | cor *= binexp.x; |
340 | eps = 1.00000000001 + (error + err_1) * binexp.x; |
341 | t = 1.0 + res; |
342 | y = ((1.0 - t) + res) + cor; |
343 | res = t + y; |
344 | cor = (t - res) + y; |
345 | if (res == (res + eps * cor)) |
346 | { |
347 | binexp.i[HIGH_HALF] = 0x00100000; |
348 | return (res - 1.0) * binexp.x; |
349 | } |
350 | else |
351 | return -10.0; |
352 | } |
353 | else |
354 | { |
355 | binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20; |
356 | if (res == (res + cor * (1.0 + error + err_1))) |
357 | return res * binexp.x * t256.x; |
358 | else |
359 | return -10.0; |
360 | } |
361 | } |
362 | |