math: explicitly promote expressions to excess-precision types
authorRich Felker <dalias@aerifal.cx>
Sat, 21 Nov 2015 21:23:30 +0000 (21:23 +0000)
committerRich Felker <dalias@aerifal.cx>
Sun, 22 Nov 2015 02:41:42 +0000 (21:41 -0500)
commit8eead3ef18ea71a64ef3cbff8c09bac3b82f1242
tree6494875b77a28a8cc62b3c7606d8801a32497030
parent19caa25d0a8e587bb89b79c3f629085548709dd4
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
src/math/__rem_pio2f.c
src/math/hypot.c