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/* 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
49double __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. */
53double
54SECTION
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 }
223ret:
224 return retval;
225
226 ret_huge:
227 return hhuge * hhuge;
228
229 ret_tiny:
230 return tiny * tiny;
231}
232#ifndef __ieee754_exp
233strong_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. */
240double
241SECTION
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