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
0ce946cf80 were likely responsible for
breaking compatibility with this and other old gcc versions.

patch by Szabolcs Nagy.

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 */

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;

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;
 }