From 78419467a2f88744ae2445fca5eb442877ebb1b0 Mon Sep 17 00:00:00 2001
From: Elliott Hughes <enh@google.com>
Date: Wed, 12 Jun 2013 16:37:58 -0700
Subject: [PATCH] Take upstream libm changes.
Mostly workarounds for GCC and Clang bugs.
Change-Id: I4ef428a42d4ac6d622659053711a8cc416925727
---
libm/upstream-freebsd/lib/msun/src/e_acosh.c | 6 +
libm/upstream-freebsd/lib/msun/src/e_atanh.c | 6 +
libm/upstream-freebsd/lib/msun/src/e_exp.c | 2 +-
libm/upstream-freebsd/lib/msun/src/e_expf.c | 5 +-
libm/upstream-freebsd/lib/msun/src/e_log.c | 9 +-
libm/upstream-freebsd/lib/msun/src/e_log10.c | 9 +-
libm/upstream-freebsd/lib/msun/src/e_log10f.c | 3 +-
libm/upstream-freebsd/lib/msun/src/e_log2.c | 9 +-
libm/upstream-freebsd/lib/msun/src/e_log2f.c | 3 +-
libm/upstream-freebsd/lib/msun/src/e_logf.c | 3 +-
.../lib/msun/src/math_private.h | 292 ++++++++++++++++++
libm/upstream-freebsd/lib/msun/src/s_asinh.c | 6 +
libm/upstream-freebsd/lib/msun/src/s_exp2.c | 5 +-
libm/upstream-freebsd/lib/msun/src/s_exp2f.c | 5 +-
libm/upstream-freebsd/lib/msun/src/s_expm1.c | 7 +-
libm/upstream-freebsd/lib/msun/src/s_expm1f.c | 3 +-
libm/upstream-freebsd/lib/msun/src/s_fma.c | 8 +-
libm/upstream-freebsd/lib/msun/src/s_fmal.c | 8 +-
libm/upstream-freebsd/lib/msun/src/s_log1p.c | 7 +-
libm/upstream-freebsd/lib/msun/src/s_log1pf.c | 3 +-
.../lib/msun/src/s_nearbyint.c | 6 +-
21 files changed, 383 insertions(+), 22 deletions(-)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_acosh.c b/libm/upstream-freebsd/lib/msun/src/e_acosh.c
index a0cc6cb55..358c8bd6a 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_acosh.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_acosh.c
@@ -29,6 +29,8 @@ __FBSDID("$FreeBSD$");
* acosh(NaN) is NaN without signal.
*/
+#include <float.h>
+
#include "math.h"
#include "math_private.h"
@@ -60,3 +62,7 @@ __ieee754_acosh(double x)
return log1p(t+sqrt(2.0*t+t*t));
}
}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(acosh, acoshl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/e_atanh.c b/libm/upstream-freebsd/lib/msun/src/e_atanh.c
index ab8a2e181..422ff2698 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_atanh.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_atanh.c
@@ -33,6 +33,8 @@ __FBSDID("$FreeBSD$");
*
*/
+#include <float.h>
+
#include "math.h"
#include "math_private.h"
@@ -60,3 +62,7 @@ __ieee754_atanh(double x)
t = 0.5*log1p((x+x)/(one-x));
if(hx>=0) return t; else return -t;
}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(atanh, atanhl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/e_exp.c b/libm/upstream-freebsd/lib/msun/src/e_exp.c
index e432bc83c..94c9769da 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_exp.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_exp.c
@@ -84,7 +84,6 @@ __FBSDID("$FreeBSD$");
static const double
one = 1.0,
halF[2] = {0.5,-0.5,},
-huge = 1.0e+300,
o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
@@ -99,6 +98,7 @@ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
static volatile double
+huge = 1.0e+300,
twom1000= 9.33263618503218878990e-302; /* 2**-1000=0x01700000,0*/
double
diff --git a/libm/upstream-freebsd/lib/msun/src/e_expf.c b/libm/upstream-freebsd/lib/msun/src/e_expf.c
index a4790765b..b1fe2c537 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_expf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_expf.c
@@ -24,7 +24,6 @@ __FBSDID("$FreeBSD$");
static const float
one = 1.0,
halF[2] = {0.5,-0.5,},
-huge = 1.0e+30,
o_threshold= 8.8721679688e+01, /* 0x42b17180 */
u_threshold= -1.0397208405e+02, /* 0xc2cff1b5 */
ln2HI[2] ={ 6.9314575195e-01, /* 0x3f317200 */
@@ -39,7 +38,9 @@ invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
P1 = 1.6666625440e-1, /* 0xaaaa8f.0p-26 */
P2 = -2.7667332906e-3; /* -0xb55215.0p-32 */
-static volatile float twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */
+static volatile float
+huge = 1.0e+30,
+twom100 = 7.8886090522e-31; /* 2**-100=0x0d800000 */
float
__ieee754_expf(float x)
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log.c b/libm/upstream-freebsd/lib/msun/src/e_log.c
index 204fb48e4..68bc1070b 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log.c
@@ -65,6 +65,8 @@ __FBSDID("$FreeBSD$");
* to produce the hexadecimal values shown.
*/
+#include <float.h>
+
#include "math.h"
#include "math_private.h"
@@ -81,6 +83,7 @@ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
+static volatile double vzero = 0.0;
double
__ieee754_log(double x)
@@ -94,7 +97,7 @@ __ieee754_log(double x)
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
- return -two54/zero; /* log(+-0)=-inf */
+ return -two54/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx,x);
@@ -138,3 +141,7 @@ __ieee754_log(double x)
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log, logl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log10.c b/libm/upstream-freebsd/lib/msun/src/e_log10.c
index 104d2574c..3c89ed2d8 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log10.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log10.c
@@ -22,6 +22,8 @@ __FBSDID("$FreeBSD$");
* in not-quite-routine extra precision.
*/
+#include <float.h>
+
#include "math.h"
#include "math_private.h"
#include "k_log.h"
@@ -34,6 +36,7 @@ log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
static const double zero = 0.0;
+static volatile double vzero = 0.0;
double
__ieee754_log10(double x)
@@ -47,7 +50,7 @@ __ieee754_log10(double x)
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
- return -two54/zero; /* log(+-0)=-inf */
+ return -two54/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx,x);
@@ -85,3 +88,7 @@ __ieee754_log10(double x)
return val_lo + val_hi;
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log10, log10l);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log10f.c b/libm/upstream-freebsd/lib/msun/src/e_log10f.c
index c8765941c..9856df2e7 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log10f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log10f.c
@@ -28,6 +28,7 @@ log10_2hi = 3.0102920532e-01, /* 0x3e9a2080 */
log10_2lo = 7.9034151668e-07; /* 0x355427db */
static const float zero = 0.0;
+static volatile float vzero = 0.0;
float
__ieee754_log10f(float x)
@@ -40,7 +41,7 @@ __ieee754_log10f(float x)
k=0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff)==0)
- return -two25/zero; /* log(+-0)=-inf */
+ return -two25/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(hx,x);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log2.c b/libm/upstream-freebsd/lib/msun/src/e_log2.c
index 1fc44a5fe..4766cdb49 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log2.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log2.c
@@ -24,6 +24,8 @@ __FBSDID("$FreeBSD$");
* in not-quite-routine extra precision.
*/
+#include <float.h>
+
#include "math.h"
#include "math_private.h"
#include "k_log.h"
@@ -34,6 +36,7 @@ ivln2hi = 1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
ivln2lo = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
static const double zero = 0.0;
+static volatile double vzero = 0.0;
double
__ieee754_log2(double x)
@@ -47,7 +50,7 @@ __ieee754_log2(double x)
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
- return -two54/zero; /* log(+-0)=-inf */
+ return -two54/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54; x *= two54; /* subnormal number, scale up x */
GET_HIGH_WORD(hx,x);
@@ -108,3 +111,7 @@ __ieee754_log2(double x)
return val_lo + val_hi;
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log2, log2l);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/e_log2f.c b/libm/upstream-freebsd/lib/msun/src/e_log2f.c
index 716634664..1794484e6 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_log2f.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_log2f.c
@@ -26,6 +26,7 @@ ivln2hi = 1.4428710938e+00, /* 0x3fb8b000 */
ivln2lo = -1.7605285393e-04; /* 0xb9389ad4 */
static const float zero = 0.0;
+static volatile float vzero = 0.0;
float
__ieee754_log2f(float x)
@@ -38,7 +39,7 @@ __ieee754_log2f(float x)
k=0;
if (hx < 0x00800000) { /* x < 2**-126 */
if ((hx&0x7fffffff)==0)
- return -two25/zero; /* log(+-0)=-inf */
+ return -two25/vzero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(hx,x);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_logf.c b/libm/upstream-freebsd/lib/msun/src/e_logf.c
index c3be6eda5..ec3985fcb 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_logf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_logf.c
@@ -30,6 +30,7 @@ Lg3 = 0x91e9ee.0p-25, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
static const float zero = 0.0;
+static volatile float vzero = 0.0;
float
__ieee754_logf(float x)
@@ -42,7 +43,7 @@ __ieee754_logf(float x)
k=0;
if (ix < 0x00800000) { /* x < 2**-126 */
if ((ix&0x7fffffff)==0)
- return -two25/zero; /* log(+-0)=-inf */
+ return -two25/vzero; /* log(+-0)=-inf */
if (ix<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 25; x *= two25; /* subnormal number, scale up x */
GET_FLOAT_WORD(ix,x);
diff --git a/libm/upstream-freebsd/lib/msun/src/math_private.h b/libm/upstream-freebsd/lib/msun/src/math_private.h
index 5662df01a..8ebc7fbf3 100644
--- a/libm/upstream-freebsd/lib/msun/src/math_private.h
+++ b/libm/upstream-freebsd/lib/msun/src/math_private.h
@@ -188,6 +188,33 @@ do { \
(d) = sf_u.value; \
} while (0)
+/*
+ * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
+ * double.
+ */
+
+#define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \
+do { \
+ union IEEEl2bits ew_u; \
+ ew_u.e = (d); \
+ (ix0) = ew_u.xbits.expsign; \
+ (ix1) = ew_u.xbits.man; \
+} while (0)
+
+/*
+ * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
+ * long double.
+ */
+
+#define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \
+do { \
+ union IEEEl2bits ew_u; \
+ ew_u.e = (d); \
+ (ix0) = ew_u.xbits.expsign; \
+ (ix1) = ew_u.xbits.manh; \
+ (ix2) = ew_u.xbits.manl; \
+} while (0)
+
/* Get expsign as a 16 bit int from a long double. */
#define GET_LDBL_EXPSIGN(i,d) \
@@ -197,6 +224,33 @@ do { \
(i) = ge_u.xbits.expsign; \
} while (0)
+/*
+ * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
+ * mantissa.
+ */
+
+#define INSERT_LDBL80_WORDS(d,ix0,ix1) \
+do { \
+ union IEEEl2bits iw_u; \
+ iw_u.xbits.expsign = (ix0); \
+ iw_u.xbits.man = (ix1); \
+ (d) = iw_u.e; \
+} while (0)
+
+/*
+ * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
+ * comprising the mantissa.
+ */
+
+#define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \
+do { \
+ union IEEEl2bits iw_u; \
+ iw_u.xbits.expsign = (ix0); \
+ iw_u.xbits.manh = (ix1); \
+ iw_u.xbits.manl = (ix2); \
+ (d) = iw_u.e; \
+} while (0)
+
/* Set expsign of a long double from a 16 bit int. */
#define SET_LDBL_EXPSIGN(d,v) \
@@ -260,6 +314,110 @@ do { \
/* Default return statement if hack*_t() is not used. */
#define RETURNF(v) return (v)
+/*
+ * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
+ * a == 0, but is slower.
+ */
+#define _2sum(a, b) do { \
+ __typeof(a) __s, __w; \
+ \
+ __w = (a) + (b); \
+ __s = __w - (a); \
+ (b) = ((a) - (__w - __s)) + ((b) - __s); \
+ (a) = __w; \
+} while (0)
+
+/*
+ * 2sumF algorithm.
+ *
+ * "Normalize" the terms in the infinite-precision expression a + b for
+ * the sum of 2 floating point values so that b is as small as possible
+ * relative to 'a'. (The resulting 'a' is the value of the expression in
+ * the same precision as 'a' and the resulting b is the rounding error.)
+ * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
+ * exponent overflow or underflow must not occur. This uses a Theorem of
+ * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum"
+ * is apparently due to Skewchuk (1997).
+ *
+ * For this to always work, assignment of a + b to 'a' must not retain any
+ * extra precision in a + b. This is required by C standards but broken
+ * in many compilers. The brokenness cannot be worked around using
+ * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
+ * algorithm would be destroyed by non-null strict assignments. (The
+ * compilers are correct to be broken -- the efficiency of all floating
+ * point code calculations would be destroyed similarly if they forced the
+ * conversions.)
+ *
+ * Fortunately, a case that works well can usually be arranged by building
+ * any extra precision into the type of 'a' -- 'a' should have type float_t,
+ * double_t or long double. b's type should be no larger than 'a's type.
+ * Callers should use these types with scopes as large as possible, to
+ * reduce their own extra-precision and efficiciency problems. In
+ * particular, they shouldn't convert back and forth just to call here.
+ */
+#ifdef DEBUG
+#define _2sumF(a, b) do { \
+ __typeof(a) __w; \
+ volatile __typeof(a) __ia, __ib, __r, __vw; \
+ \
+ __ia = (a); \
+ __ib = (b); \
+ assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \
+ \
+ __w = (a) + (b); \
+ (b) = ((a) - __w) + (b); \
+ (a) = __w; \
+ \
+ /* The next 2 assertions are weak if (a) is already long double. */ \
+ assert((long double)__ia + __ib == (long double)(a) + (b)); \
+ __vw = __ia + __ib; \
+ __r = __ia - __vw; \
+ __r += __ib; \
+ assert(__vw == (a) && __r == (b)); \
+} while (0)
+#else /* !DEBUG */
+#define _2sumF(a, b) do { \
+ __typeof(a) __w; \
+ \
+ __w = (a) + (b); \
+ (b) = ((a) - __w) + (b); \
+ (a) = __w; \
+} while (0)
+#endif /* DEBUG */
+
+/*
+ * Set x += c, where x is represented in extra precision as a + b.
+ * x must be sufficiently normalized and sufficiently larger than c,
+ * and the result is then sufficiently normalized.
+ *
+ * The details of ordering are that |a| must be >= |c| (so that (a, c)
+ * can be normalized without extra work to swap 'a' with c). The details of
+ * the normalization are that b must be small relative to the normalized 'a'.
+ * Normalization of (a, c) makes the normalized c tiny relative to the
+ * normalized a, so b remains small relative to 'a' in the result. However,
+ * b need not ever be tiny relative to 'a'. For example, b might be about
+ * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
+ * That is usually enough, and adding c (which by normalization is about
+ * 2**53 times smaller than a) cannot change b significantly. However,
+ * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
+ * significantly relative to b. The caller must ensure that significant
+ * cancellation doesn't occur, either by having c of the same sign as 'a',
+ * or by having |c| a few percent smaller than |a|. Pre-normalization of
+ * (a, b) may help.
+ *
+ * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
+ * exercise 19). We gain considerable efficiency by requiring the terms to
+ * be sufficiently normalized and sufficiently increasing.
+ */
+#define _3sumF(a, b, c) do { \
+ __typeof(a) __tmp; \
+ \
+ __tmp = (c); \
+ _2sumF(__tmp, (a)); \
+ (b) += (a); \
+ (a) = __tmp; \
+} while (0)
+
/*
* Common routine to process the arguments to nan(), nanf(), and nanl().
*/
@@ -370,6 +528,140 @@ irintl(long double x)
#endif /* __GNUCLIKE_ASM */
+#ifdef DEBUG
+#if defined(__amd64__) || defined(__i386__)
+#define breakpoint() asm("int $3")
+#else
+#include <signal.h>
+
+#define breakpoint() raise(SIGTRAP)
+#endif
+#endif
+
+/* Write a pari script to test things externally. */
+#ifdef DOPRINT
+#include <stdio.h>
+
+#ifndef DOPRINT_SWIZZLE
+#define DOPRINT_SWIZZLE 0
+#endif
+
+#ifdef DOPRINT_LD80
+
+#define DOPRINT_START(xp) do { \
+ uint64_t __lx; \
+ uint16_t __hx; \
+ \
+ /* Hack to give more-problematic args. */ \
+ EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
+ __lx ^= DOPRINT_SWIZZLE; \
+ INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
+ printf("x = %.21Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_D64)
+
+#define DOPRINT_START(xp) do { \
+ uint32_t __hx, __lx; \
+ \
+ EXTRACT_WORDS(__hx, __lx, *xp); \
+ __lx ^= DOPRINT_SWIZZLE; \
+ INSERT_WORDS(*xp, __hx, __lx); \
+ printf("x = %.21Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_F32)
+
+#define DOPRINT_START(xp) do { \
+ uint32_t __hx; \
+ \
+ GET_FLOAT_WORD(__hx, *xp); \
+ __hx ^= DOPRINT_SWIZZLE; \
+ SET_FLOAT_WORD(*xp, __hx); \
+ printf("x = %.21Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
+
+#ifndef DOPRINT_SWIZZLE_HIGH
+#define DOPRINT_SWIZZLE_HIGH 0
+#endif
+
+#define DOPRINT_START(xp) do { \
+ uint64_t __lx, __llx; \
+ uint16_t __hx; \
+ \
+ EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
+ __llx ^= DOPRINT_SWIZZLE; \
+ __lx ^= DOPRINT_SWIZZLE_HIGH; \
+ INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
+ printf("x = %.36Lg; ", (long double)*xp); \
+} while (0)
+#define DOPRINT_END1(v) \
+ printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define DOPRINT_END2(hi, lo) \
+ printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
+ (long double)(hi), (long double)(lo))
+
+#endif /* DOPRINT_LD80 */
+
+#else /* !DOPRINT */
+#define DOPRINT_START(xp)
+#define DOPRINT_END1(v)
+#define DOPRINT_END2(hi, lo)
+#endif /* DOPRINT */
+
+#define RETURNP(x) do { \
+ DOPRINT_END1(x); \
+ RETURNF(x); \
+} while (0)
+#define RETURNPI(x) do { \
+ DOPRINT_END1(x); \
+ RETURNI(x); \
+} while (0)
+#define RETURN2P(x, y) do { \
+ DOPRINT_END2((x), (y)); \
+ RETURNF((x) + (y)); \
+} while (0)
+#define RETURN2PI(x, y) do { \
+ DOPRINT_END2((x), (y)); \
+ RETURNI((x) + (y)); \
+} while (0)
+#ifdef STRUCT_RETURN
+#define RETURNSP(rp) do { \
+ if (!(rp)->lo_set) \
+ RETURNP((rp)->hi); \
+ RETURN2P((rp)->hi, (rp)->lo); \
+} while (0)
+#define RETURNSPI(rp) do { \
+ if (!(rp)->lo_set) \
+ RETURNPI((rp)->hi); \
+ RETURN2PI((rp)->hi, (rp)->lo); \
+} while (0)
+#endif
+#define SUM2P(x, y) ({ \
+ const __typeof (x) __x = (x); \
+ const __typeof (y) __y = (y); \
+ \
+ DOPRINT_END2(__x, __y); \
+ __x + __y; \
+})
+
/*
* ieee style elementary functions
*
diff --git a/libm/upstream-freebsd/lib/msun/src/s_asinh.c b/libm/upstream-freebsd/lib/msun/src/s_asinh.c
index f3fdf741a..cbb3d463a 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_asinh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_asinh.c
@@ -24,6 +24,8 @@ __FBSDID("$FreeBSD$");
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
+#include <float.h>
+
#include "math.h"
#include "math_private.h"
@@ -54,3 +56,7 @@ asinh(double x)
}
if(hx>0) return w; else return -w;
}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(asinh, asinhl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_exp2.c b/libm/upstream-freebsd/lib/msun/src/s_exp2.c
index 485b4e3c3..fde11c27d 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_exp2.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_exp2.c
@@ -36,7 +36,6 @@ __FBSDID("$FreeBSD$");
#define TBLSIZE (1 << TBLBITS)
static const double
- huge = 0x1p1000,
redux = 0x1.8p52 / TBLSIZE,
P1 = 0x1.62e42fefa39efp-1,
P2 = 0x1.ebfbdff82c575p-3,
@@ -44,7 +43,9 @@ static const double
P4 = 0x1.3b2ab88f70400p-7,
P5 = 0x1.5d88003875c74p-10;
-static volatile double twom1000 = 0x1p-1000;
+static volatile double
+ huge = 0x1p1000,
+ twom1000 = 0x1p-1000;
static const double tbl[TBLSIZE * 2] = {
/* exp2(z + eps) eps */
diff --git a/libm/upstream-freebsd/lib/msun/src/s_exp2f.c b/libm/upstream-freebsd/lib/msun/src/s_exp2f.c
index 0a97bf64d..9ac7c1f6c 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_exp2f.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_exp2f.c
@@ -36,14 +36,15 @@ __FBSDID("$FreeBSD$");
#define TBLSIZE (1 << TBLBITS)
static const float
- huge = 0x1p100f,
redux = 0x1.8p23f / TBLSIZE,
P1 = 0x1.62e430p-1f,
P2 = 0x1.ebfbe0p-3f,
P3 = 0x1.c6b348p-5f,
P4 = 0x1.3b2c9cp-7f;
-static volatile float twom100 = 0x1p-100f;
+static volatile float
+ huge = 0x1p100f,
+ twom100 = 0x1p-100f;
static const double exp2ft[TBLSIZE] = {
0x1.6a09e667f3bcdp-1,
diff --git a/libm/upstream-freebsd/lib/msun/src/s_expm1.c b/libm/upstream-freebsd/lib/msun/src/s_expm1.c
index 5aa1917b5..37998a387 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_expm1.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_expm1.c
@@ -115,7 +115,6 @@ __FBSDID("$FreeBSD$");
static const double
one = 1.0,
-huge = 1.0e+300,
tiny = 1.0e-300,
o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
@@ -128,6 +127,8 @@ Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
+static volatile double huge = 1.0e+300;
+
double
expm1(double x)
{
@@ -215,3 +216,7 @@ expm1(double x)
}
return y;
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(expm1, expm1l);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_expm1f.c b/libm/upstream-freebsd/lib/msun/src/s_expm1f.c
index fb3749430..c0a39340f 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_expm1f.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_expm1f.c
@@ -23,7 +23,6 @@ __FBSDID("$FreeBSD$");
static const float
one = 1.0,
-huge = 1.0e+30,
tiny = 1.0e-30,
o_threshold = 8.8721679688e+01,/* 0x42b17180 */
ln2_hi = 6.9313812256e-01,/* 0x3f317180 */
@@ -37,6 +36,8 @@ invln2 = 1.4426950216e+00,/* 0x3fb8aa3b */
Q1 = -3.3333212137e-2, /* -0x888868.0p-28 */
Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */
+static volatile float huge = 1.0e+30;
+
float
expm1f(float x)
{
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fma.c b/libm/upstream-freebsd/lib/msun/src/s_fma.c
index dfbd13c46..452bece37 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fma.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fma.c
@@ -238,6 +238,8 @@ fma(double x, double y, double z)
zs = copysign(DBL_MIN, zs);
fesetround(FE_TONEAREST);
+ /* work around clang bug 8100 */
+ volatile double vxs = xs;
/*
* Basic approach for round-to-nearest:
@@ -247,7 +249,7 @@ fma(double x, double y, double z)
* adj = xy.lo + r.lo (inexact; low bit is sticky)
* result = r.hi + adj (correctly rounded)
*/
- xy = dd_mul(xs, ys);
+ xy = dd_mul(vxs, ys);
r = dd_add(xy.hi, zs);
spread = ex + ey;
@@ -268,7 +270,9 @@ fma(double x, double y, double z)
* rounding modes.
*/
fesetround(oround);
- adj = r.lo + xy.lo;
+ /* work around clang bug 8100 */
+ volatile double vrlo = r.lo;
+ adj = vrlo + xy.lo;
return (ldexp(r.hi + adj, spread));
}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmal.c b/libm/upstream-freebsd/lib/msun/src/s_fmal.c
index c2a691350..92719010e 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmal.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmal.c
@@ -226,6 +226,8 @@ fmal(long double x, long double y, long double z)
zs = copysignl(LDBL_MIN, zs);
fesetround(FE_TONEAREST);
+ /* work around clang bug 8100 */
+ volatile long double vxs = xs;
/*
* Basic approach for round-to-nearest:
@@ -235,7 +237,7 @@ fmal(long double x, long double y, long double z)
* adj = xy.lo + r.lo (inexact; low bit is sticky)
* result = r.hi + adj (correctly rounded)
*/
- xy = dd_mul(xs, ys);
+ xy = dd_mul(vxs, ys);
r = dd_add(xy.hi, zs);
spread = ex + ey;
@@ -256,7 +258,9 @@ fmal(long double x, long double y, long double z)
* rounding modes.
*/
fesetround(oround);
- adj = r.lo + xy.lo;
+ /* work around clang bug 8100 */
+ volatile long double vrlo = r.lo;
+ adj = vrlo + xy.lo;
return (ldexpl(r.hi + adj, spread));
}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_log1p.c b/libm/upstream-freebsd/lib/msun/src/s_log1p.c
index b062a8ad9..3cc77bda2 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_log1p.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_log1p.c
@@ -96,6 +96,7 @@ Lp6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static const double zero = 0.0;
+static volatile double vzero = 0.0;
double
log1p(double x)
@@ -109,7 +110,7 @@ log1p(double x)
k = 1;
if (hx < 0x3FDA827A) { /* 1+x < sqrt(2)+ */
if(ax>=0x3ff00000) { /* x <= -1.0 */
- if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */
+ if(x==-1.0) return -two54/vzero; /* log1p(-1)=+inf */
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if(ax<0x3e200000) { /* |x| < 2**-29 */
@@ -173,3 +174,7 @@ log1p(double x)
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
}
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log1p, log1pl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_log1pf.c b/libm/upstream-freebsd/lib/msun/src/s_log1pf.c
index 01d3457d3..df04c6756 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_log1pf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_log1pf.c
@@ -34,6 +34,7 @@ Lp6 = 1.5313838422e-01, /* 3E1CD04F */
Lp7 = 1.4798198640e-01; /* 3E178897 */
static const float zero = 0.0;
+static volatile float vzero = 0.0;
float
log1pf(float x)
@@ -47,7 +48,7 @@ log1pf(float x)
k = 1;
if (hx < 0x3ed413d0) { /* 1+x < sqrt(2)+ */
if(ax>=0x3f800000) { /* x <= -1.0 */
- if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
+ if(x==(float)-1.0) return -two25/vzero; /* log1p(-1)=+inf */
else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
}
if(ax<0x38000000) { /* |x| < 2**-15 */
diff --git a/libm/upstream-freebsd/lib/msun/src/s_nearbyint.c b/libm/upstream-freebsd/lib/msun/src/s_nearbyint.c
index 12493d292..063f8d736 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_nearbyint.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_nearbyint.c
@@ -36,12 +36,16 @@ __FBSDID("$FreeBSD$");
* instead of feclearexcept()/feupdateenv() to restore the environment
* because the only exception defined for rint() is overflow, and
* rounding can't overflow as long as emax >= p.
+ *
+ * The volatile keyword is needed below because clang incorrectly assumes
+ * that rint won't raise any floating-point exceptions. Declaring ret volatile
+ * is sufficient to trick the compiler into doing the right thing.
*/
#define DECL(type, fn, rint) \
type \
fn(type x) \
{ \
- type ret; \
+ volatile type ret; \
fenv_t env; \
\
fegetenv(&env); \