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.

diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c
index a40db9fc05fdd028d90f11d8f9ce8f8db6992fcf..d403f81c792217a195d9680d9356eec17f734f41 100644 (file)
--- 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 f5163856512880fdf82597b169e0b75d6ae0257c..4473c1c420e922167814e2602eb889be62d2b283 100644 (file)
--- 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 29ec6a476d1bf6787b393c0f2a498ea338b583a6..6071bf1e284f376f417dbab5bc306964a7491f3f 100644 (file)
--- 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;
 }