| s_expl.c (251334) | s_expl.c (251335) |
|---|---|
| 1/*- 2 * Copyright (c) 2009-2013 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-2013 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 251334 2013-06-03 18:40:00Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 251335 2013-06-03 18:51:34Z 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 * --- 202 unchanged lines hidden (view full) --- 241 uint16_t hx, ix; 242 243 /* Filter out exceptional cases. */ 244 u.e = x; 245 hx = u.xbits.expsign; 246 ix = hx & 0x7fff; 247 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 248 if (ix == BIAS + LDBL_MAX_EXP) { | 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 * --- 202 unchanged lines hidden (view full) --- 241 uint16_t hx, ix; 242 243 /* Filter out exceptional cases. */ 244 u.e = x; 245 hx = u.xbits.expsign; 246 ix = hx & 0x7fff; 247 if (ix >= BIAS + 13) { /* |x| >= 8192 or x is NaN */ 248 if (ix == BIAS + LDBL_MAX_EXP) { |
| 249 if (hx & 0x8000 && u.xbits.man == 1ULL << 63) 250 return (0.0L); /* x is -Inf */ 251 return (x + x); /* x is +Inf, NaN or unsupported */ | 249 if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */ 250 return (-1 / x); 251 return (x + x); /* x is +Inf, +NaN or unsupported */ |
| 252 } 253 if (x > o_threshold) 254 return (huge * huge); 255 if (x < u_threshold) 256 return (tiny * tiny); | 252 } 253 if (x > o_threshold) 254 return (huge * huge); 255 if (x < u_threshold) 256 return (tiny * tiny); |
| 257 } else if (ix < BIAS - 66) { /* |x| < 0x1p-66 */ 258 /* includes pseudo-denormals */ 259 if (huge + x > 1.0L) /* trigger inexact iff x != 0 */ 260 return (1.0L + x); | 257 } else if (ix < BIAS - 65) { /* |x| < 0x1p-65 (includes pseudos) */ 258 return (1 + x); /* 1 with inexact iff x != 0 */ |
| 261 } 262 263 ENTERI(); 264 265 /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */ 266 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 267 fn = x * INV_L + 0x1.8p63 - 0x1.8p63; 268 r = x - fn * L1 - fn * L2; /* r = r1 + r2 done independently. */ --- 38 unchanged lines hidden --- | 259 } 260 261 ENTERI(); 262 263 /* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */ 264 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ 265 fn = x * INV_L + 0x1.8p63 - 0x1.8p63; 266 r = x - fn * L1 - fn * L2; /* r = r1 + r2 done independently. */ --- 38 unchanged lines hidden --- |