Re: for information, gcc-4.2.3 miscompiles musl math

[Szabolcs Nagy <nsz@port70.net>]

[-- Attachment #1: Type: text/plain, Size: 973 bytes --]

* Rich Felker <dalias@libc.org> [2015-11-21 15:15:51 -0500]:
> 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.
> 

i did the changes where i think it was necessary
(hypot has asm implementation on i386 though).

(in most cases the code already uses float_t, double_t
where this would matter.)

[-- Attachment #2: 0001-math-make-ffloat-store-work-on-i386.patch --]
[-- Type: text/x-diff, Size: 2257 bytes --]

From c4cf94d37e2c556a2f77b1debe2fe9ec269fdb3e Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Sat, 21 Nov 2015 21:23:30 +0000
Subject: [PATCH] math: make -ffloat-store work on i386

old i386 toolchains don't support -fexcess-precision=standard
and then musl configure defaults to -ffloat-store to avoid
problems because of excess precision.

however there are cases when excess precision must not be rounded
away so -ffloat-store breaks the code.

this patch adds (double_t) casts to double precision arithmetics
when the excess precision is required for correct results.
(without -ffloat-store or on non-i386 targets the generated code
is not changed.)
---
 src/math/__rem_pio2.c  | 2 +-
 src/math/__rem_pio2f.c | 2 +-
 src/math/hypot.c       | 4 ++--
 3 files changed, 4 insertions(+), 4 deletions(-)

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 */
diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c
index f516385..4473c1c 100644
--- 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 29ec6a4..6071bf1 100644
--- 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;
 }
 
-- 
2.4.1