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.

3 files changed, 4 insertions, 4 deletions

diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c
index a40db9fc..d403f81c 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 f5163856..4473c1c4 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 29ec6a47..6071bf1e 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;
 }