[PATCH] Fix i386 frexpl

[Jakub Jelinek]

Hi!

i386 assembly frexpl.S has a bug in test for denormals (where C code tests
        if (ix==0x0000) {               /* subnormal */
            x *= two65;
            GET_LDOUBLE_EXP(se,x);
            ix = se&0x7fff;
            *eptr = -65;
        }
assembly has
        cmpl    $0, %eax
	je      2f
	special denormal handling
2:
(%eax holds ix)), so frexpl was needlessly slow for most normalized numbers
and was broken for denormals and values with large exponents (the former
because the needed special care was not given to it, the latter because it
overflowed). In December last year frexpl got a workaround for the large
exponent case, but denormals are still broken.
The following patch backs out 2000-12-03 changes (leaves fwait in) and fixes
the branch, plus adds testcase for this which passed glibc make check.

2001-06-02  Jakub Jelinek  <jakub@redhat.com>

	* sysdeps/i386/fpu/s_frexpl.S (__frexpl): Mostly revert 2000-12-03
	changes, do the special handling for denormal numbers, not for
	normalized numbers (patch by <trevin@xmission.com>).

	* math/test-misc.c (main): Test frexpl with denormal arguments.

--- libc/math/test-misc.c.jj	Wed Jan 31 16:35:03 2001
+++ libc/math/test-misc.c	Sat Jun  2 20:36:01 2001
@@ -79,6 +79,28 @@ main (void)
 	  }
 	puts ("ok");
       }
+
+    for (i = LDBL_MIN_EXP, x = LDBL_MIN; i >= LDBL_MIN_EXP - LDBL_MANT_DIG + 1;
+	 --i, x /= 2.0L)
+      {
+        printf ("2^%d: ", i);
+
+        r = frexpl (x, &e);
+        if (r != 0.5L)
+          {
+            printf ("mantissa incorrect: %.20La\n", r);
+            result = 1;
+            continue;
+          }
+        if (e != i)
+          {
+            printf ("exponent wrong %d (%.20Lg)\n", e, x);
+            result = 1;
+            continue;
+          }
+        puts ("ok");
+      }
+
   }
 # endif
 
--- libc/sysdeps/i386/fpu/s_frexpl.S.jj	Mon Dec  4 15:07:46 2000
+++ libc/sysdeps/i386/fpu/s_frexpl.S	Sat Jun  2 19:56:15 2001
@@ -1,5 +1,5 @@
 /* ix87 specific frexp implementation for long double.
-   Copyright (C) 1997, 2000 Free Software Foundation, Inc.
+   Copyright (C) 1997, 2000, 2001 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -32,12 +32,6 @@
 	ASM_TYPE_DIRECTIVE(two64,@object)
 two64:	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
 	ASM_SIZE_DIRECTIVE(two64)
-	/* The following is LDBL_MAX / ldexp (1.0, 64), the largest
-	   number we can handle the normal way.  */
-	ASM_TYPE_DIRECTIVE(largest,@object)
-largest:
-	.byte 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xbe, 0x7f, 0, 0
-	ASM_SIZE_DIRECTIVE(largest)
 
 #ifdef PIC
 #define MO(op) op##@GOTOFF(%edx)
@@ -67,10 +61,7 @@ ENTRY (BP_SYM (__frexpl))
 	je	1f
 
 	cmpl	$0, %eax
-	je	2f
-
-	cmpl	$0x7fbe, %eax
-	ja	4f
+	jne	2f
 
 	fldt	VAL0(%esp)
 #ifdef	PIC
@@ -102,15 +93,5 @@ ENTRY (BP_SYM (__frexpl))
 
 	LEAVE
 	ret
-
-4:	movl	VAL2(%esp), %ecx
-	movl	%ecx, %edx
-	andl	$0x7fff, %ecx
-
-	andl	$0x8000, %edx
-	subl	$16382, %ecx
-	orl	$0x3ffe, %edx
-	movl	%edx, VAL2(%esp)
-	jmp	1b
 END (BP_SYM (__frexpl))
 weak_alias (BP_SYM (__frexpl), BP_SYM (frexpl))

	Jakub