Re: for information, gcc-4.2.3 miscompiles musl math

[Szabolcs Nagy <nsz@port70.net>]

* Rich Felker <dalias@libc.org> [2015-11-21 14:51:44 -0500]:
> On Sat, Nov 21, 2015 at 08:41:32PM +0100, Szabolcs Nagy wrote:
> > * Rich Felker <dalias@libc.org> [2015-11-21 14:25:48 -0500]:
> > 
> > > On Sat, Nov 21, 2015 at 06:24:18PM +0100, u-uy74@aetey.se wrote:
> > > > Good to be aware of:
> > > > gcc-4.2.3 miscompiles musl math since at least 1.1.6,
> > > > tested while targeting i486,
> > > > 1.0.x seems to have been alright.
> > > > 
> > > > The symptom is that sin(larger-than-2*pi) yields large values
> > > > like "sin(8.000000) = 21.709544".
> > > > Looks like the argument reduction logic has changed in a way
> > > > which is not compatible with gcc-4.2.3.
> > > 
> > > Are you using configure or a hand-written config.mak? configure sets
> > > up a big hammer, -ffloat-store, when -fexcess-precision=standard is
> > > not supported (i.e. on old gcc), which hopefully suffices to make this
> > > code work, but it's possible it doesn't always do the job.
> > > 
> > 
> > i think this change might be it:
> > http://git.musl-libc.org/cgit/musl/commit/?id=0ce946cf808274c2d6e5419b139e130c8ad4bd30
> > 
> > the new code avoids an extra store,
> > but then i rely on the evaluation
> > being in long double.
> > 
> > with -ffloat-store this breaks,
> > adding extra store rounds at the
> > wrong place.
> 
> Hmm, which places does it add the stores around? Could you fix it with
> an explicit conversion to double_t? That might be nice to harden
> against broken compilers without penalizing correct ones.
> 

yeah that might work

i dont have gcc-4.2, can you try:

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