2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2014 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 /***************************************************************************/
40 #include <math_private.h>
48 double __slowexp (double);
50 /* An ultimate exp routine. Given an IEEE double machine number x it computes
51 the correctly rounded (to nearest) value of e^x. */
54 __ieee754_exp (double x
)
56 double bexp
, t
, eps
, del
, base
, y
, al
, bet
, res
, rem
, cor
;
57 mynumber junk1
, junk2
, binexp
= {{0, 0}};
61 SET_RESTORE_ROUND (FE_TONEAREST
);
64 m
= junk1
.i
[HIGH_HALF
];
67 if (n
> smallint
&& n
< bigint
)
69 y
= x
* log2e
.x
+ three51
.x
;
70 bexp
= y
- three51
.x
; /* multiply the result by 2**bexp */
74 eps
= bexp
* ln_two2
.x
; /* x = bexp*ln(2) + t - eps */
75 t
= x
- bexp
* ln_two1
.x
;
78 base
= y
- three33
.x
; /* t rounded to a multiple of 2**-18 */
80 del
= (t
- base
) - eps
; /* x = bexp*ln(2) + base + del */
81 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
83 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 1023) << 20;
85 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
86 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
88 al
= coar
.x
[i
] * fine
.x
[j
];
89 bet
= ((coar
.x
[i
] * fine
.x
[j
+ 1] + coar
.x
[i
+ 1] * fine
.x
[j
])
90 + coar
.x
[i
+ 1] * fine
.x
[j
+ 1]);
92 rem
= (bet
+ bet
* eps
) + al
* eps
;
94 cor
= (al
- res
) + rem
;
95 if (res
== (res
+ cor
* err_0
))
97 retval
= res
* binexp
.x
;
102 retval
= __slowexp (x
);
104 } /*if error is over bound */
122 retval
= (x
> 0) ? (hhuge
* hhuge
) : (tiny
* tiny
);
125 /* x is finite, cause either overflow or underflow */
126 if (junk1
.i
[LOW_HALF
] != 0)
131 retval
= (x
> 0) ? inf
.x
: zero
; /* |x| = inf; return either inf or 0 */
135 y
= x
* log2e
.x
+ three51
.x
;
136 bexp
= y
- three51
.x
;
138 eps
= bexp
* ln_two2
.x
;
139 t
= x
- bexp
* ln_two1
.x
;
141 base
= y
- three33
.x
;
143 del
= (t
- base
) - eps
;
144 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
145 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
146 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
147 al
= coar
.x
[i
] * fine
.x
[j
];
148 bet
= ((coar
.x
[i
] * fine
.x
[j
+ 1] + coar
.x
[i
+ 1] * fine
.x
[j
])
149 + coar
.x
[i
+ 1] * fine
.x
[j
+ 1]);
150 rem
= (bet
+ bet
* eps
) + al
* eps
;
152 cor
= (al
- res
) + rem
;
155 ex
= junk1
.i
[LOW_HALF
];
164 binexp
.i
[HIGH_HALF
] = (1023 + ex
) << 20;
165 if (res
== (res
+ cor
* err_0
))
167 retval
= res
* binexp
.x
;
172 retval
= __slowexp (x
);
173 goto check_uflow_ret
;
174 } /*if error is over bound */
177 binexp
.i
[HIGH_HALF
] = (1023 - ex
) << 20;
180 eps
= 1.0000000001 + err_0
* binexp
.x
;
182 y
= ((1.0 - t
) + res
) + cor
;
185 if (res
== (res
+ eps
* cor
))
187 binexp
.i
[HIGH_HALF
] = 0x00100000;
188 retval
= (res
- 1.0) * binexp
.x
;
189 goto check_uflow_ret
;
193 retval
= __slowexp (x
);
194 goto check_uflow_ret
;
195 } /* if error is over bound */
197 if (retval
< DBL_MIN
)
199 #if FLT_EVAL_METHOD != 0
202 double force_underflow
= tiny
* tiny
;
203 math_force_eval (force_underflow
);
209 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 767) << 20;
210 if (res
== (res
+ cor
* err_0
))
212 retval
= res
* binexp
.x
* t256
.x
;
217 retval
= __slowexp (x
);
224 #ifndef __ieee754_exp
225 strong_alias (__ieee754_exp
, __exp_finite
)
228 /* Compute e^(x+xx). The routine also receives bound of error of previous
229 calculation. If after computing exp the error exceeds the allowed bounds,
230 the routine returns a non-positive number. Otherwise it returns the
231 computed result, which is always positive. */
234 __exp1 (double x
, double xx
, double error
)
236 double bexp
, t
, eps
, del
, base
, y
, al
, bet
, res
, rem
, cor
;
237 mynumber junk1
, junk2
, binexp
= {{0, 0}};
241 m
= junk1
.i
[HIGH_HALF
];
242 n
= m
& hugeint
; /* no sign */
244 if (n
> smallint
&& n
< bigint
)
246 y
= x
* log2e
.x
+ three51
.x
;
247 bexp
= y
- three51
.x
; /* multiply the result by 2**bexp */
251 eps
= bexp
* ln_two2
.x
; /* x = bexp*ln(2) + t - eps */
252 t
= x
- bexp
* ln_two1
.x
;
255 base
= y
- three33
.x
; /* t rounded to a multiple of 2**-18 */
257 del
= (t
- base
) + (xx
- eps
); /* x = bexp*ln(2) + base + del */
258 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
260 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 1023) << 20;
262 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
263 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
265 al
= coar
.x
[i
] * fine
.x
[j
];
266 bet
= ((coar
.x
[i
] * fine
.x
[j
+ 1] + coar
.x
[i
+ 1] * fine
.x
[j
])
267 + coar
.x
[i
+ 1] * fine
.x
[j
+ 1]);
269 rem
= (bet
+ bet
* eps
) + al
* eps
;
271 cor
= (al
- res
) + rem
;
272 if (res
== (res
+ cor
* (1.0 + error
+ err_1
)))
273 return res
* binexp
.x
;
279 return 1.0; /* if x->0 e^x=1 */
284 return (zero
/ zero
); /* x is NaN, return invalid */
286 return ((x
> 0) ? (hhuge
* hhuge
) : (tiny
* tiny
));
287 /* x is finite, cause either overflow or underflow */
288 if (junk1
.i
[LOW_HALF
] != 0)
289 return (zero
/ zero
); /* x is NaN */
290 return ((x
> 0) ? inf
.x
: zero
); /* |x| = inf; return either inf or 0 */
293 y
= x
* log2e
.x
+ three51
.x
;
294 bexp
= y
- three51
.x
;
296 eps
= bexp
* ln_two2
.x
;
297 t
= x
- bexp
* ln_two1
.x
;
299 base
= y
- three33
.x
;
301 del
= (t
- base
) + (xx
- eps
);
302 eps
= del
+ del
* del
* (p3
.x
* del
+ p2
.x
);
303 i
= ((junk2
.i
[LOW_HALF
] >> 8) & 0xfffffffe) + 356;
304 j
= (junk2
.i
[LOW_HALF
] & 511) << 1;
305 al
= coar
.x
[i
] * fine
.x
[j
];
306 bet
= ((coar
.x
[i
] * fine
.x
[j
+ 1] + coar
.x
[i
+ 1] * fine
.x
[j
])
307 + coar
.x
[i
+ 1] * fine
.x
[j
+ 1]);
308 rem
= (bet
+ bet
* eps
) + al
* eps
;
310 cor
= (al
- res
) + rem
;
313 ex
= junk1
.i
[LOW_HALF
];
322 binexp
.i
[HIGH_HALF
] = (1023 + ex
) << 20;
323 if (res
== (res
+ cor
* (1.0 + error
+ err_1
)))
324 return res
* binexp
.x
;
329 binexp
.i
[HIGH_HALF
] = (1023 - ex
) << 20;
332 eps
= 1.00000000001 + (error
+ err_1
) * binexp
.x
;
334 y
= ((1.0 - t
) + res
) + cor
;
337 if (res
== (res
+ eps
* cor
))
339 binexp
.i
[HIGH_HALF
] = 0x00100000;
340 return (res
- 1.0) * binexp
.x
;
347 binexp
.i
[HIGH_HALF
] = (junk1
.i
[LOW_HALF
] + 767) << 20;
348 if (res
== (res
+ cor
* (1.0 + error
+ err_1
)))
349 return res
* binexp
.x
* t256
.x
;