| s_expl.c (251325) | s_expl.c (251327) |
|---|---|
| 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 251325 2013-06-03 17:40:52Z kargl $"); | 30__FBSDID("$FreeBSD: head/lib/msun/ld80/s_expl.c 251327 2013-06-03 17:51:08Z 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 * --- 6 unchanged lines hidden (view full) --- 45#include <ieeefp.h> 46#endif 47 48#include "fpmath.h" 49#include "math.h" 50#include "math_private.h" 51 52#define INTERVALS 128 | 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 * --- 6 unchanged lines hidden (view full) --- 45#include <ieeefp.h> 46#endif 47 48#include "fpmath.h" 49#include "math.h" 50#include "math_private.h" 51 52#define INTERVALS 128 |
| 53#define LOG2_INTERVALS 7 |
|
| 53#define BIAS (LDBL_MAX_EXP - 1) 54 55static const long double 56huge = 0x1p10000L, 57twom10000 = 0x1p-10000L; 58/* XXX Prevent gcc from erroneously constant folding this: */ 59static volatile const long double tiny = 0x1p-10000L; 60 --- 203 unchanged lines hidden (view full) --- 264#if defined(HAVE_EFFICIENT_IRINTL) 265 n = irintl(fn); 266#elif defined(HAVE_EFFICIENT_IRINT) 267 n = irint(fn); 268#else 269 n = (int)fn; 270#endif 271 n2 = (unsigned)n % INTERVALS; | 54#define BIAS (LDBL_MAX_EXP - 1) 55 56static const long double 57huge = 0x1p10000L, 58twom10000 = 0x1p-10000L; 59/* XXX Prevent gcc from erroneously constant folding this: */ 60static volatile const long double tiny = 0x1p-10000L; 61 --- 203 unchanged lines hidden (view full) --- 265#if defined(HAVE_EFFICIENT_IRINTL) 266 n = irintl(fn); 267#elif defined(HAVE_EFFICIENT_IRINT) 268 n = irint(fn); 269#else 270 n = (int)fn; 271#endif 272 n2 = (unsigned)n % INTERVALS; |
| 272 k = (n - n2) / INTERVALS; | 273 /* Depend on the sign bit being propagated: */ 274 k = n >> LOG2_INTERVALS; |
| 273 r1 = x - fn * L1; 274 r2 = -fn * L2; 275 276 /* Prepare scale factors. */ 277 v.xbits.man = 1ULL << 63; 278 if (k >= LDBL_MIN_EXP) { 279 v.xbits.expsign = BIAS + k; 280 twopk = v.e; --- 23 unchanged lines hidden --- | 275 r1 = x - fn * L1; 276 r2 = -fn * L2; 277 278 /* Prepare scale factors. */ 279 v.xbits.man = 1ULL << 63; 280 if (k >= LDBL_MIN_EXP) { 281 v.xbits.expsign = BIAS + k; 282 twopk = v.e; --- 23 unchanged lines hidden --- |