2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
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.
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.
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/>.
19 /***************************************************************************/
20 /* MODULE_NAME:uexp.c */
25 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */
26 /* mpa.c mpexp.x slowexp.c */
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. */
33 /***************************************************************************/
41 #include <math_private.h>
49 double __slowexp (double);
51 /* An ultimate exp routine. Given an IEEE double machine number x it computes
52 the correctly rounded (to nearest) value of e^x. */
55 __ieee754_exp (double x
)
57 double bexp
, t
, eps
, del
, base
, y
, al
, bet
, res
, rem
, cor
;
58 mynumber junk1
, junk2
, binexp
= {{0, 0}};
63 SET_RESTORE_ROUND (FE_TONEAREST
);
66 m
= junk1
.i
[HIGH_HALF
];
69 if (n
> smallint
&& n
< bigint
)
71 y
= x
* log2e
.x
+ three51
.x
;
72 bexp
= y
- three51
.x
; /* multiply the result by 2**bexp */
76 eps
= bexp
* ln_two2
.x
; /* x = bexp*ln(2) + t - eps */
77 t
= x
- bexp
* ln_two1
.x
;
80 base
= y
- three33
.x
; /* t rounded to a multiple of 2**-18 */
82 del
= (t
- base
) - eps
; /* x = bexp*ln(2) + base + del */
83 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
85 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 1023) << 20;
87 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
88 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
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]);
94 rem
= (bet
+ bet
* eps
) + al
* eps
;
96 cor
= (al
- res
) + rem
;
97 if (res
== (res
+ cor
* err_0
))
99 retval
= res
* binexp
.x
;
104 retval
= __slowexp (x
);
106 } /*if error is over bound */
129 /* x is finite, cause either overflow or underflow */
130 if (junk1
.i
[LOW_HALF
] != 0)
135 retval
= (x
> 0) ? inf
.x
: zero
; /* |x| = inf; return either inf or 0 */
139 y
= x
* log2e
.x
+ three51
.x
;
140 bexp
= y
- three51
.x
;
142 eps
= bexp
* ln_two2
.x
;
143 t
= x
- bexp
* ln_two1
.x
;
145 base
= y
- three33
.x
;
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
;
156 cor
= (al
- res
) + rem
;
159 ex
= junk1
.i
[LOW_HALF
];
168 binexp
.i
[HIGH_HALF
] = (1023 + ex
) << 20;
169 if (res
== (res
+ cor
* err_0
))
171 retval
= res
* binexp
.x
;
176 retval
= __slowexp (x
);
177 goto check_uflow_ret
;
178 } /*if error is over bound */
181 binexp
.i
[HIGH_HALF
] = (1023 - ex
) << 20;
184 eps
= 1.0000000001 + err_0
* binexp
.x
;
186 y
= ((1.0 - t
) + res
) + cor
;
189 if (res
== (res
+ eps
* cor
))
191 binexp
.i
[HIGH_HALF
] = 0x00100000;
192 retval
= (res
- 1.0) * binexp
.x
;
193 goto check_uflow_ret
;
197 retval
= __slowexp (x
);
198 goto check_uflow_ret
;
199 } /* if error is over bound */
201 if (retval
< DBL_MIN
)
203 #if FLT_EVAL_METHOD != 0
206 double force_underflow
= tiny
* tiny
;
207 math_force_eval (force_underflow
);
215 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 767) << 20;
216 if (res
== (res
+ cor
* err_0
))
217 retval
= res
* binexp
.x
* t256
.x
;
219 retval
= __slowexp (x
);
230 return hhuge
* hhuge
;
235 #ifndef __ieee754_exp
236 strong_alias (__ieee754_exp
, __exp_finite
)
239 /* Compute e^(x+xx). The routine also receives bound of error of previous
240 calculation. If after computing exp the error exceeds the allowed bounds,
241 the routine returns a non-positive number. Otherwise it returns the
242 computed result, which is always positive. */
245 __exp1 (double x
, double xx
, double error
)
247 double bexp
, t
, eps
, del
, base
, y
, al
, bet
, res
, rem
, cor
;
248 mynumber junk1
, junk2
, binexp
= {{0, 0}};
252 m
= junk1
.i
[HIGH_HALF
];
253 n
= m
& hugeint
; /* no sign */
255 if (n
> smallint
&& n
< bigint
)
257 y
= x
* log2e
.x
+ three51
.x
;
258 bexp
= y
- three51
.x
; /* multiply the result by 2**bexp */
262 eps
= bexp
* ln_two2
.x
; /* x = bexp*ln(2) + t - eps */
263 t
= x
- bexp
* ln_two1
.x
;
266 base
= y
- three33
.x
; /* t rounded to a multiple of 2**-18 */
268 del
= (t
- base
) + (xx
- eps
); /* x = bexp*ln(2) + base + del */
269 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
271 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 1023) << 20;
273 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
274 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
276 al
= coar
.x
[i
] * fine
.x
[j
];
277 bet
= ((coar
.x
[i
] * fine
.x
[j
+ 1] + coar
.x
[i
+ 1] * fine
.x
[j
])
278 + coar
.x
[i
+ 1] * fine
.x
[j
+ 1]);
280 rem
= (bet
+ bet
* eps
) + al
* eps
;
282 cor
= (al
- res
) + rem
;
283 if (res
== (res
+ cor
* (1.0 + error
+ err_1
)))
284 return res
* binexp
.x
;
290 return 1.0; /* if x->0 e^x=1 */
295 return (zero
/ zero
); /* x is NaN, return invalid */
297 return ((x
> 0) ? (hhuge
* hhuge
) : (tiny
* tiny
));
298 /* x is finite, cause either overflow or underflow */
299 if (junk1
.i
[LOW_HALF
] != 0)
300 return (zero
/ zero
); /* x is NaN */
301 return ((x
> 0) ? inf
.x
: zero
); /* |x| = inf; return either inf or 0 */
304 y
= x
* log2e
.x
+ three51
.x
;
305 bexp
= y
- three51
.x
;
307 eps
= bexp
* ln_two2
.x
;
308 t
= x
- bexp
* ln_two1
.x
;
310 base
= y
- three33
.x
;
312 del
= (t
- base
) + (xx
- eps
);
313 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
314 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
315 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
316 al
= coar
.x
[i
] * fine
.x
[j
];
317 bet
= ((coar
.x
[i
] * fine
.x
[j
+ 1] + coar
.x
[i
+ 1] * fine
.x
[j
])
318 + coar
.x
[i
+ 1] * fine
.x
[j
+ 1]);
319 rem
= (bet
+ bet
* eps
) + al
* eps
;
321 cor
= (al
- res
) + rem
;
324 ex
= junk1
.i
[LOW_HALF
];
333 binexp
.i
[HIGH_HALF
] = (1023 + ex
) << 20;
334 if (res
== (res
+ cor
* (1.0 + error
+ err_1
)))
335 return res
* binexp
.x
;
340 binexp
.i
[HIGH_HALF
] = (1023 - ex
) << 20;
343 eps
= 1.00000000001 + (error
+ err_1
) * binexp
.x
;
345 y
= ((1.0 - t
) + res
) + cor
;
348 if (res
== (res
+ eps
* cor
))
350 binexp
.i
[HIGH_HALF
] = 0x00100000;
351 return (res
- 1.0) * binexp
.x
;
358 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 767) << 20;
359 if (res
== (res
+ cor
* (1.0 + error
+ err_1
)))
360 return res
* binexp
.x
* t256
.x
;