| s_expl.c (240864) | s_expl.c (240865) |
|---|---|
| 1/*- 2 * Copyright (c) 2009-2012 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 13 unchanged lines hidden (view full) --- 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 * 26 * Optimized by Bruce D. Evans. 27 */ 28 29#include <sys/cdefs.h> | 1/*- 2 * Copyright (c) 2009-2012 Steven G. Kargl 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 13 unchanged lines hidden (view full) --- 22 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 23 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 24 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 25 * 26 * Optimized by Bruce D. Evans. 27 */ 28 29#include <sys/cdefs.h> |
| 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 240864 2012-09-23 18:06:27Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 240865 2012-09-23 18:13:46Z kargl $"); |
| 31 32/*- 33 * Compute the exponential of x for Intel 80-bit format. This is based on: 34 * 35 * PTP Tang, "Table-driven implementation of the exponential function 36 * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 37 * 144-157 (1989). 38 * --- 207 unchanged lines hidden (view full) --- 246 return (x + x); /* x is +Inf, NaN or unsupported */ 247 } 248 if (x > o_threshold.e) 249 return (huge * huge); 250 if (x < u_threshold.e) 251 return (tiny * tiny); 252 } else if (ix <= BIAS - 34) { /* |x| < 0x1p-33 */ 253 /* includes pseudo-denormals */ | 31 32/*- 33 * Compute the exponential of x for Intel 80-bit format. This is based on: 34 * 35 * PTP Tang, "Table-driven implementation of the exponential function 36 * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 37 * 144-157 (1989). 38 * --- 207 unchanged lines hidden (view full) --- 246 return (x + x); /* x is +Inf, NaN or unsupported */ 247 } 248 if (x > o_threshold.e) 249 return (huge * huge); 250 if (x < u_threshold.e) 251 return (tiny * tiny); 252 } else if (ix <= BIAS - 34) { /* |x| < 0x1p-33 */ 253 /* includes pseudo-denormals */ |
| 254 if (huge + x > 1.0L) /* trigger inexact iff x != 0 */ | 254 if (huge + x > 1.0L) /* trigger inexact iff x != 0 */ |
| 255 return (1.0L + x); 256 } 257 258 ENTERI(); 259 260 /* Reduce x to (k*ln2 + midpoint[n2] + r1 + r2). */ 261 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 262 fn = x * INV_L + 0x1.8p63 - 0x1.8p63; --- 41 unchanged lines hidden --- | 255 return (1.0L + x); 256 } 257 258 ENTERI(); 259 260 /* Reduce x to (k*ln2 + midpoint[n2] + r1 + r2). */ 261 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 262 fn = x * INV_L + 0x1.8p63 - 0x1.8p63; --- 41 unchanged lines hidden --- |