From 8eead3ef18ea71a64ef3cbff8c09bac3b82f1242 Mon Sep 17 00:00:00 2001 From: Rich Felker Date: Sat, 21 Nov 2015 21:23:30 +0000 Subject: [PATCH] math: explicitly promote expressions to excess-precision types a conforming compiler for an arch with excess precision floating point (FLT_EVAL_METHOD!=0; presently i386 is the only such arch supported) computes all intermediate results in the types float_t and double_t rather than the nominal type of the expression. some incorrect compilers, however, only keep excess precision in registers, and convert down to the nominal type when spilling intermediate results to memory, yielding unpredictable results that depend on the compiler's choices of what/when to spill. in particular, this happens on old gcc versions with -ffloat-store, which we need in order to work around bugs where the compiler wrongly keeps explicitly-dropped excess precision. by explicitly converting to double_t where expressions are expected be be evaluated in double_t precision, we can avoid depending on the compiler to get types correct when spilling; the nominal and intermediate precision now match. this commit should not change the code generated by correct compilers, or by old ones on non-i386 archs where double_t is defined as double. this fixes a serious bug in argument reduction observed on i386 with gcc 4.2: for values of x outside the unit circle, sin(x) was producing results outside the interval [-1,1]. changes made in commit 0ce946cf808274c2d6e5419b139e130c8ad4bd30 were likely responsible for breaking compatibility with this and other old gcc versions. patch by Szabolcs Nagy. --- src/math/__rem_pio2.c | 2 +- src/math/__rem_pio2f.c | 2 +- src/math/hypot.c | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index a40db9f..d403f81 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -118,7 +118,7 @@ int __rem_pio2(double x, double *y) if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ medium: /* rint(x/(pi/2)), Assume round-to-nearest. */ - fn = x*invpio2 + toint - toint; + fn = (double_t)x*invpio2 + toint - toint; n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c index f516385..4473c1c 100644 --- a/src/math/__rem_pio2f.c +++ b/src/math/__rem_pio2f.c @@ -51,7 +51,7 @@ int __rem_pio2f(float x, double *y) /* 25+53 bit pi is good enough for medium size */ if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - fn = x*invpio2 + toint - toint; + fn = (double_t)x*invpio2 + toint - toint; n = (int32_t)fn; *y = x - fn*pio2_1 - fn*pio2_1t; return n; diff --git a/src/math/hypot.c b/src/math/hypot.c index 29ec6a4..6071bf1 100644 --- a/src/math/hypot.c +++ b/src/math/hypot.c @@ -12,10 +12,10 @@ static void sq(double_t *hi, double_t *lo, double x) { double_t xh, xl, xc; - xc = x*SPLIT; + xc = (double_t)x*SPLIT; xh = x - xc + xc; xl = x - xh; - *hi = x*x; + *hi = (double_t)x*x; *lo = xh*xh - *hi + 2*xh*xl + xl*xl; } -- 2.7.0