ieee754[sd]p_neg workaround

ieee754[sd]p_neg workaround

[Atsushi Nemoto]

I have a long standing patch for FPU emulator to fix a segmentation
fault in pow() library function.

Here is a test program to reproduce it.

main()
{
	union {
		double d;
		struct {
#ifdef __MIPSEB
			unsigned int high, low;
#else
			unsigned int low, high;
#endif
		} i;
	} x, y, z;
        x.i.low = 0x00000000;
        x.i.high = 0xfff00001;
        y.i.low = 0x80000000;
        y.i.high = 0xcff00000;
        z.d = pow(x.d, y.d);
        printf("%x %x\n", z.i.high, z.i.low);
        return 0;
}


If you run this program, you will get segmentation fault (unless your
FPU does not raise Unimplemented exception for NaN operands).  The
segmentation fault is caused by endless recursion in __ieee754_pow().

It looks glibc's pow() assume unary '-' operation for any number
(including NaN) always invert its sign bit.

Here is a revised (just added a few comments) patch.  Please review.
Thank you.

diff -u linux-mips/arch/mips/math-emu/dp_simple.c linux/arch/mips/math-emu/dp_simple.c
--- linux-mips/arch/mips/math-emu/dp_simple.c	2005-01-31 11:05:18.000000000 +0900
+++ linux/arch/mips/math-emu/dp_simple.c	2005-04-20 17:02:02.613112541 +0900
@@ -48,16 +48,22 @@
 	CLEARCX;
 	FLUSHXDP;
 
+	/*
+	 * Invert the sign ALWAYS to prevent an endless recursion on
+	 * pow() in libc.
+	 */
+	/* quick fix up */
+	DPSIGN(x) ^= 1;
+
 	if (xc == IEEE754_CLASS_SNAN) {
+		ieee754dp y = ieee754dp_indef();
 		SETCX(IEEE754_INVALID_OPERATION);
-		return ieee754dp_nanxcpt(ieee754dp_indef(), "neg");
+		DPSIGN(y) = DPSIGN(x);
+		return ieee754dp_nanxcpt(y, "neg");
 	}
 
 	if (ieee754dp_isnan(x))	/* but not infinity */
 		return ieee754dp_nanxcpt(x, "neg", x);
-
-	/* quick fix up */
-	DPSIGN(x) ^= 1;
 	return x;
 }
 
diff -u linux-mips/arch/mips/math-emu/sp_simple.c linux/arch/mips/math-emu/sp_simple.c
--- linux-mips/arch/mips/math-emu/sp_simple.c	2005-01-31 11:05:18.000000000 +0900
+++ linux/arch/mips/math-emu/sp_simple.c	2005-04-20 17:02:13.678391113 +0900
@@ -48,16 +48,22 @@
 	CLEARCX;
 	FLUSHXSP;
 
+	/*
+	 * Invert the sign ALWAYS to prevent an endless recursion on
+	 * pow() in libc.
+	 */
+	/* quick fix up */
+	SPSIGN(x) ^= 1;
+
 	if (xc == IEEE754_CLASS_SNAN) {
+		ieee754sp y = ieee754sp_indef();
 		SETCX(IEEE754_INVALID_OPERATION);
-		return ieee754sp_nanxcpt(ieee754sp_indef(), "neg");
+		SPSIGN(y) = SPSIGN(x);
+		return ieee754sp_nanxcpt(y, "neg");
 	}
 
 	if (ieee754sp_isnan(x))	/* but not infinity */
 		return ieee754sp_nanxcpt(x, "neg", x);
-
-	/* quick fix up */
-	SPSIGN(x) ^= 1;
 	return x;
 }
 

---
Atsushi Nemoto

Re: ieee754[sd]p_neg workaround

[Maciej W. Rozycki]

On Wed, 20 Apr 2005, Atsushi Nemoto wrote:

> I have a long standing patch for FPU emulator to fix a segmentation
> fault in pow() library function.
> 
> Here is a test program to reproduce it.
> 
> main()
> {
> 	union {
> 		double d;
> 		struct {
> #ifdef __MIPSEB
> 			unsigned int high, low;
> #else
> 			unsigned int low, high;
> #endif
> 		} i;
> 	} x, y, z;
>         x.i.low = 0x00000000;
>         x.i.high = 0xfff00001;
>         y.i.low = 0x80000000;
>         y.i.high = 0xcff00000;
>         z.d = pow(x.d, y.d);
>         printf("%x %x\n", z.i.high, z.i.low);
>         return 0;
> }
> 
> 
> If you run this program, you will get segmentation fault (unless your
> FPU does not raise Unimplemented exception for NaN operands).  The
> segmentation fault is caused by endless recursion in __ieee754_pow().
> 
> It looks glibc's pow() assume unary '-' operation for any number
> (including NaN) always invert its sign bit.

 AFAICS, the IEEE 754 standard explicitly leaves interpretation of the 
sign bit for NaNs as unspecified.  Therefore our implementation is correct 
and its glibc that should be fixed instead.  Please file a bug report 
against glibc.

  Maciej

Re: ieee754[sd]p_neg workaround

[Dominic Sweetman]

Maciej W. Rozycki (macro@linux-mips.org) writes:

> > I have a long standing patch for FPU emulator to fix a segmentation
> > fault in pow() library function.
> > 
> > Here is a test program to reproduce it.
> > 
> > main()
> > {
> > 	union {
> > 		double d;
> > 		struct {
> > #ifdef __MIPSEB
> > 			unsigned int high, low;
> > #else
> > 			unsigned int low, high;
> > #endif
> > 		} i;
> > 	} x, y, z;
> >         x.i.low = 0x00000000;
> >         x.i.high = 0xfff00001;
> >         y.i.low = 0x80000000;
> >         y.i.high = 0xcff00000;
> >         z.d = pow(x.d, y.d);
> >         printf("%x %x\n", z.i.high, z.i.low);
> >         return 0;
> > }
> > 
> > 
> > If you run this program, you will get segmentation fault (unless your
> > FPU does not raise Unimplemented exception for NaN operands).  The
> > segmentation fault is caused by endless recursion in __ieee754_pow().
> > 
> > It looks glibc's pow() assume unary '-' operation for any number
> > (including NaN) always invert its sign bit.
> 
>  AFAICS, the IEEE 754 standard explicitly leaves interpretation of the 
> sign bit for NaNs as unspecified.  Therefore our implementation is correct 
> and its glibc that should be fixed instead.  Please file a bug report 
> against glibc.

You are both right, in a sense.

"Annex A Recommended Functions and Predicates" of IEEE754 (the
annexe is not part of the formal spec) includes a recommendation:

   "­x is x copied with its sign reversed, not 0­x; the distinction is
    germane when x is ±0 or NaN."

And in fact that's how the MIPS neg.d operation works: it never
generates an exception, just blindly flips the sign bit.

In the body of IEEE754 it says:

   "This standard does not interpret the sign of an NaN."

So there you are.  According to the book, the library function should
not depend on the sign of a NaN, but it would also be better if the
compiler/emulator/whatever ensures that "-x" is always and only
implemented by reversing the sign bit.

So file a bug against glibc, but we should fix the emulator so it
correctly imitates the MIPS instruction set...

--
Dominic Sweetman
MIPS Technologies.

Re: ieee754[sd]p_neg workaround

[Ralf Baechle]

On Wed, Apr 20, 2005 at 01:57:41PM +0100, Dominic Sweetman wrote:

> So file a bug against glibc, but we should fix the emulator so it
> correctly imitates the MIPS instruction set...

As a matter of defensive design I think we should try to follow the
establish behaviour if nothing more specific is defined anywhere.

  Ralf

Re: ieee754[sd]p_neg workaround

[Atsushi Nemoto]

>>>>> On Wed, 20 Apr 2005 14:13:04 +0100, Ralf Baechle <ralf@linux-mips.org> said:
>> So file a bug against glibc, but we should fix the emulator so it
>> correctly imitates the MIPS instruction set...

ralf> As a matter of defensive design I think we should try to follow
ralf> the establish behaviour if nothing more specific is defined
ralf> anywhere.

OK, I sent a bug reoport to glibc bugzilla. (Bug# 864)

---
Atsushi Nemoto

Re: ieee754[sd]p_neg workaround

[Atsushi Nemoto]

>>>>> On Thu, 21 Apr 2005 16:19:45 +0900 (JST), Atsushi Nemoto <anemo@mba.ocn.ne.jp> said:

>>> So file a bug against glibc, but we should fix the emulator so it
>>> correctly imitates the MIPS instruction set...

ralf> As a matter of defensive design I think we should try to follow
ralf> the establish behaviour if nothing more specific is defined
ralf> anywhere.

anemo> OK, I sent a bug reoport to glibc bugzilla. (Bug# 864)

The bug was resolved ... marked as WONTFIX.

glibc developers said math-emu should be fixed.  Please look at:

http://sourceware.org/bugzilla/show_bug.cgi?id=864

for their comments.

So I send a patch for math-emu again.

Description:

It looks glibc's pow() assume an unary '-' operation for any number
(including NaNs) always invert its sign bit (though IEEE754 does not
specify the sign bit for NaNs).  This patch make the kernel math-emu
emulates real MIPS neg.[ds] instruction.

Signed-off-by: Atsushi Nemoto <anemo@mba.ocn.ne.jp>

diff --git a/arch/mips/math-emu/dp_simple.c b/arch/mips/math-emu/dp_simple.c
index 495c1ac..1c555e6 100644
--- a/arch/mips/math-emu/dp_simple.c
+++ b/arch/mips/math-emu/dp_simple.c
@@ -48,16 +48,22 @@ ieee754dp ieee754dp_neg(ieee754dp x)
 	CLEARCX;
 	FLUSHXDP;
 
+	/*
+	 * Invert the sign ALWAYS to prevent an endless recursion on
+	 * pow() in libc.
+	 */
+	/* quick fix up */
+	DPSIGN(x) ^= 1;
+
 	if (xc == IEEE754_CLASS_SNAN) {
+		ieee754dp y = ieee754dp_indef();
 		SETCX(IEEE754_INVALID_OPERATION);
-		return ieee754dp_nanxcpt(ieee754dp_indef(), "neg");
+		DPSIGN(y) = DPSIGN(x);
+		return ieee754dp_nanxcpt(y, "neg");
 	}
 
 	if (ieee754dp_isnan(x))	/* but not infinity */
 		return ieee754dp_nanxcpt(x, "neg", x);
-
-	/* quick fix up */
-	DPSIGN(x) ^= 1;
 	return x;
 }
 
diff --git a/arch/mips/math-emu/sp_simple.c b/arch/mips/math-emu/sp_simple.c
index c809830..770f0f4 100644
--- a/arch/mips/math-emu/sp_simple.c
+++ b/arch/mips/math-emu/sp_simple.c
@@ -48,16 +48,22 @@ ieee754sp ieee754sp_neg(ieee754sp x)
 	CLEARCX;
 	FLUSHXSP;
 
+	/*
+	 * Invert the sign ALWAYS to prevent an endless recursion on
+	 * pow() in libc.
+	 */
+	/* quick fix up */
+	SPSIGN(x) ^= 1;
+
 	if (xc == IEEE754_CLASS_SNAN) {
+		ieee754sp y = ieee754sp_indef();
 		SETCX(IEEE754_INVALID_OPERATION);
-		return ieee754sp_nanxcpt(ieee754sp_indef(), "neg");
+		SPSIGN(y) = SPSIGN(x);
+		return ieee754sp_nanxcpt(y, "neg");
 	}
 
 	if (ieee754sp_isnan(x))	/* but not infinity */
 		return ieee754sp_nanxcpt(x, "neg", x);
-
-	/* quick fix up */
-	SPSIGN(x) ^= 1;
 	return x;
 }