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/* */
27/* An ultimate exp routine. Given an IEEE double machine number x */
28/* it computes an almost correctly rounded (to nearest) value of e^x */
29/* Assumption: Machine arithmetic operations are performed in */
30/* round to nearest mode of IEEE 754 standard. */
31/* */
32/***************************************************************************/
33
34#include <math.h>
35#include "endian.h"
36#include "uexp.h"
37#include "mydefs.h"
38#include "MathLib.h"
39#include "uexp.tbl"
40#include <math-barriers.h>
41#include <math_private.h>
42#include <fenv.h>
43#include <float.h>
44#include "eexp.tbl"
45
46#ifndef SECTION
47# define SECTION
48#endif
49
50double
51SECTION
52__ieee754_exp (double x)
53{
54 double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
55 double z;
56 mynumber junk1, junk2, binexp = {{0, 0}};
57 int4 i, j, m, n, ex;
58 int4 k;
59 double retval;
60
61 {
62 SET_RESTORE_ROUND (FE_TONEAREST);
63
64 junk1.x = x;
65 m = junk1.i[HIGH_HALF];
66 n = m & hugeint;
67
68 if (n < 0x3ff0a2b2) /* |x| < 1.03972053527832 */
69 {
70 if (n < 0x3f862e42) /* |x| < 3/2 ln 2 */
71 {
72 if (n < 0x3ed00000) /* |x| < 1/64 ln 2 */
73 {
74 if (n < 0x3e300000) /* |x| < 2^18 */
75 {
76 retval = one + junk1.x;
77 goto ret;
78 }
79 retval = one + junk1.x * (one + half * junk1.x);
80 goto ret;
81 }
82 t = junk1.x * junk1.x;
83 retval = junk1.x + (t * (half + junk1.x * t2) +
84 (t * t) * (t3 + junk1.x * t4 + t * t5));
85 retval = one + retval;
86 goto ret;
87 }
88
89 /* Find the multiple of 2^-6 nearest x. */
90 k = n >> 20;
91 j = (0x00100000 | (n & 0x000fffff)) >> (0x40c - k);
92 j = (j - 1) & ~1;
93 if (m < 0)
94 j += 134;
95 z = junk1.x - TBL2[j];
96 t = z * z;
97 retval = z + (t * (half + (z * t2))
98 + (t * t) * (t3 + z * t4 + t * t5));
99 retval = TBL2[j + 1] + TBL2[j + 1] * retval;
100 goto ret;
101 }
102
103 if (n < bigint) /* && |x| >= 1.03972053527832 */
104 {
105 y = x * log2e.x + three51.x;
106 bexp = y - three51.x; /* multiply the result by 2**bexp */
107
108 junk1.x = y;
109
110 eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
111 t = x - bexp * ln_two1.x;
112
113 y = t + three33.x;
114 base = y - three33.x; /* t rounded to a multiple of 2**-18 */
115 junk2.x = y;
116 del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
117 eps = del + del * del * (p3.x * del + p2.x);
118
119 binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
120
121 i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
122 j = (junk2.i[LOW_HALF] & 511) << 1;
123
124 al = coar.x[i] * fine.x[j];
125 bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
126 + coar.x[i + 1] * fine.x[j + 1]);
127
128 rem = (bet + bet * eps) + al * eps;
129 res = al + rem;
130 /* Maximum relative error is 7.8e-22 (70.1 bits).
131 Maximum ULP error is 0.500007. */
132 retval = res * binexp.x;
133 goto ret;
134 }
135
136 if (n >= badint)
137 {
138 if (n > infint)
139 {
140 retval = x + x;
141 goto ret;
142 } /* x is NaN */
143 if (n < infint)
144 {
145 if (x > 0)
146 goto ret_huge;
147 else
148 goto ret_tiny;
149 }
150 /* x is finite, cause either overflow or underflow */
151 if (junk1.i[LOW_HALF] != 0)
152 {
153 retval = x + x;
154 goto ret;
155 } /* x is NaN */
156 retval = (x > 0) ? inf.x : zero; /* |x| = inf; return either inf or 0 */
157 goto ret;
158 }
159
160 y = x * log2e.x + three51.x;
161 bexp = y - three51.x;
162 junk1.x = y;
163 eps = bexp * ln_two2.x;
164 t = x - bexp * ln_two1.x;
165 y = t + three33.x;
166 base = y - three33.x;
167 junk2.x = y;
168 del = (t - base) - eps;
169 eps = del + del * del * (p3.x * del + p2.x);
170 i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
171 j = (junk2.i[LOW_HALF] & 511) << 1;
172 al = coar.x[i] * fine.x[j];
173 bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
174 + coar.x[i + 1] * fine.x[j + 1]);
175 rem = (bet + bet * eps) + al * eps;
176 res = al + rem;
177 cor = (al - res) + rem;
178 if (m >> 31)
179 {
180 ex = junk1.i[LOW_HALF];
181 if (res < 1.0)
182 {
183 res += res;
184 cor += cor;
185 ex -= 1;
186 }
187 if (ex >= -1022)
188 {
189 binexp.i[HIGH_HALF] = (1023 + ex) << 20;
190 /* Does not underflow: res >= 1.0, binexp >= 0x1p-1022
191 Maximum relative error is 7.8e-22 (70.1 bits).
192 Maximum ULP error is 0.500007. */
193 retval = res * binexp.x;
194 goto ret;
195 }
196 ex = -(1022 + ex);
197 binexp.i[HIGH_HALF] = (1023 - ex) << 20;
198 res *= binexp.x;
199 cor *= binexp.x;
200 t = 1.0 + res;
201 y = ((1.0 - t) + res) + cor;
202 res = t + y;
203 /* Maximum ULP error is 0.5000035. */
204 binexp.i[HIGH_HALF] = 0x00100000;
205 retval = (res - 1.0) * binexp.x;
206 if (retval < DBL_MIN)
207 {
208 double force_underflow = tiny * tiny;
209 math_force_eval (force_underflow);
210 }
211 if (retval == 0)
212 goto ret_tiny;
213 goto ret;
214 }
215 else
216 {
217 binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
218 /* Maximum relative error is 7.8e-22 (70.1 bits).
219 Maximum ULP error is 0.500007. */
220 retval = res * binexp.x * t256.x;
221 if (isinf (retval))
222 goto ret_huge;
223 else
224 goto ret;
225 }
226 }
227ret:
228 return retval;
229
230 ret_huge:
231 return hhuge * hhuge;
232
233 ret_tiny:
234 return tiny * tiny;
235}
236#ifndef __ieee754_exp
237strong_alias (__ieee754_exp, __exp_finite)
238#endif
239
240/* Compute e^(x+xx). */
241double
242SECTION
243__exp1 (double x, double xx)
244{
245 double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
246 mynumber junk1, junk2, binexp = {{0, 0}};
247 int4 i, j, m, n, ex;
248
249 junk1.x = x;
250 m = junk1.i[HIGH_HALF];
251 n = m & hugeint; /* no sign */
252
253 /* fabs (x) > 5.551112e-17 and fabs (x) < 7.080010e+02. */
254 if (n > smallint && n < bigint)
255 {
256 y = x * log2e.x + three51.x;
257 bexp = y - three51.x; /* multiply the result by 2**bexp */
258
259 junk1.x = y;
260
261 eps = bexp * ln_two2.x; /* x = bexp*ln(2) + t - eps */
262 t = x - bexp * ln_two1.x;
263
264 y = t + three33.x;
265 base = y - three33.x; /* t rounded to a multiple of 2**-18 */
266 junk2.x = y;
267 del = (t - base) + (xx - eps); /* x = bexp*ln(2) + base + del */
268 eps = del + del * del * (p3.x * del + p2.x);
269
270 binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
271
272 i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
273 j = (junk2.i[LOW_HALF] & 511) << 1;
274
275 al = coar.x[i] * fine.x[j];
276 bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
277 + coar.x[i + 1] * fine.x[j + 1]);
278
279 rem = (bet + bet * eps) + al * eps;
280 res = al + rem;
281 /* Maximum relative error before rounding is 8.8e-22 (69.9 bits).
282 Maximum ULP error is 0.500008. */
283 return res * binexp.x;
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 /* x < 0. */
322 ex = junk1.i[LOW_HALF];
323 if (res < 1.0)
324 {
325 res += res;
326 cor += cor;
327 ex -= 1;
328 }
329 if (ex >= -1022)
330 {
331 binexp.i[HIGH_HALF] = (1023 + ex) << 20;
332 /* Maximum ULP error is 0.500008. */
333 return res * binexp.x;
334 }
335 /* Denormal case - ex < -1022. */
336 ex = -(1022 + ex);
337 binexp.i[HIGH_HALF] = (1023 - ex) << 20;
338 res *= binexp.x;
339 cor *= binexp.x;
340 t = 1.0 + res;
341 y = ((1.0 - t) + res) + cor;
342 res = t + y;
343 binexp.i[HIGH_HALF] = 0x00100000;
344 /* Maximum ULP error is 0.500004. */
345 return (res - 1.0) * binexp.x;
346 }
347 else
348 {
349 binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
350 /* Maximum ULP error is 0.500008. */
351 return res * binexp.x * t256.x;
352 }
353}
354