endless recursive loop in pow()

[Atsushi Nemoto <nemoto@toshiba-tops.co.jp>]

I found that pow() of glibc 2.2.4 causes stack overflow on some
testcase with NaN operand.  Try this simple (and dirty) test program
(for little endian).

main()
{
        double x, y, z;
        *(int*)(&x) = 0x00000000;
        *((int*)(&x) + 1) = 0xfff00001;
        *(int*)(&y) = 0x80000000;
        *((int*)(&y) + 1) = 0xcff00000;
        z = pow(x, y);
        printf("%x %x\n", *((int*)&z + 1), *(int*)&z);
        return 0;
}

This program fall into endless recursive call at line 128 in
sysdeps/ieee754/dbl-64/e_pow.c.

128:     return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */


The problem occur if sign-bit of a NaN and negation of the NaN is
same.  I saw this problem only if the neg.d instruction was
interpreted by kernel fp emulation.

The patch below is a fix for kernel fp emulation.  Is this a right
fix?  Or glibc should be fixed?

---
Atsushi Nemoto
diff -u linux-sgi-cvs/arch/mips/math-emu/dp_simple.c linux.new/arch/mips/math-emu/dp_simple.c
--- linux-sgi-cvs/arch/mips/math-emu/dp_simple.c	Mon Oct 22 10:29:56 2001
+++ linux.new/arch/mips/math-emu/dp_simple.c	Tue Feb  5 17:36:29 2002
@@ -48,16 +48,18 @@
 	CLEARCX;
 	FLUSHXDP;
 
+	/* 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-sgi-cvs/arch/mips/math-emu/sp_simple.c linux.new/arch/mips/math-emu/sp_simple.c
--- linux-sgi-cvs/arch/mips/math-emu/sp_simple.c	Mon Oct 22 10:29:56 2001
+++ linux.new/arch/mips/math-emu/sp_simple.c	Tue Feb  5 17:36:29 2002
@@ -48,16 +48,18 @@
 	CLEARCX;
 	FLUSHXSP;
 
+	/* 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;
 }