Re: for information, gcc-4.2.3 miscompiles musl math

[Rich Felker <dalias@libc.org>]

On Sat, Nov 21, 2015 at 09:03:42PM +0100, Szabolcs Nagy wrote:
> * 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 */

I just tested with 4.2.1 binaries from Aboriginal Linux and (1) was
able to reproduce the bug, and (2) made it go away with your patch.
Not only did the bogus out-of-range result go away; the new result is
a bit-exact match for modern gcc.

Conceptually, it seems to me that for code that explicitly uses
float_t and double_t in expressions and only converts down to float or
double with stores, -ffloat-store should be just as good (albeit
overly pessimizing) as -fexcess-precision=standard. So IMO this looks
like a really good solution, and might come in handy as hardening
against mistakes in new compilers too. IIRC firm used to do this wrong
and I would not be at all surprised if pcc gets it wrong.

Rich