Fix x86/x86_64 expm1l spurious underflow exceptions (bug 16539)

Message ID Pine.LNX.4.64.1406242046500.21652@digraph.polyomino.org.uk
State Committed
Headers

Commit Message

Joseph Myers June 24, 2014, 8:47 p.m. UTC
  This patch fixes bug 16539, spurious underflow exceptions from x86 /
x86-64 expm1l.  The problem is that the computation of a base-2
exponent with extra precision involves spurious underflows for
arguments that are small but not subnormal, so a check is added to
just return the argument in those cases.  (If the argument *is*
subnormal, underflowing is correct and the existing code will always
underflow, so it suffices to keep using the existing code in that
case; some expm1 implementations have a bug (bug 16353) with missing
underflow exceptions, but I don't think there's such a bug in this
particular version.)

Tested x86_64 and x86; no ulps updates needed.

(auto-libm-test-out diffs omitted below.)

2014-06-24  Joseph Myers  <joseph@codesourcery.com>

	[BZ #16539]
	* sysdeps/i386/fpu/e_expl.S (IEEE754_EXPL) [USE_AS_EXPM1L]: Just
	return the argument for normal arguments with exponent below -64.
	* sysdeps/x86_64/fpu/e_expl.S (IEEE754_EXPL) [USE_AS_EXPM1L]:
	Likewise.
	* math/auto-libm-test-in: Add another test of expm1.
	* math/auto-libm-test-out: Regenerated.
  

Comments

Andreas Jaeger June 24, 2014, 8:50 p.m. UTC | #1
On 06/24/2014 10:47 PM, Joseph S. Myers wrote:
> This patch fixes bug 16539, spurious underflow exceptions from x86 /
> x86-64 expm1l.  The problem is that the computation of a base-2
> exponent with extra precision involves spurious underflows for
> arguments that are small but not subnormal, so a check is added to
> just return the argument in those cases.  (If the argument *is*
> subnormal, underflowing is correct and the existing code will always
> underflow, so it suffices to keep using the existing code in that
> case; some expm1 implementations have a bug (bug 16353) with missing
> underflow exceptions, but I don't think there's such a bug in this
> particular version.)

Thanks,
Andreas
  
Stefan Liebler June 26, 2014, 1:10 p.m. UTC | #2
Hi Joseph,

on S390 test-ldouble fails after this patch:
Failure: expm1 (0x4.0000000000000028p-16384): Exception "Underflow" set
Failure: expm1_downward (0x4.0000000000000028p-16384): Exception 
"Underflow" set
Failure: expm1_towardzero (0x4.0000000000000028p-16384): Exception 
"Underflow" set
Failure: expm1_upward (0x4.0000000000000028p-16384): Exception 
"Underflow" set

But I havenĀ“t debugged it yet.

Bye

On 06/24/2014 10:47 PM, Joseph S. Myers wrote:
> This patch fixes bug 16539, spurious underflow exceptions from x86 /
> x86-64 expm1l.  The problem is that the computation of a base-2
> exponent with extra precision involves spurious underflows for
> arguments that are small but not subnormal, so a check is added to
> just return the argument in those cases.  (If the argument *is*
> subnormal, underflowing is correct and the existing code will always
> underflow, so it suffices to keep using the existing code in that
> case; some expm1 implementations have a bug (bug 16353) with missing
> underflow exceptions, but I don't think there's such a bug in this
> particular version.)
>
> Tested x86_64 and x86; no ulps updates needed.
>
> (auto-libm-test-out diffs omitted below.)
>
> 2014-06-24  Joseph Myers  <joseph@codesourcery.com>
>
> 	[BZ #16539]
> 	* sysdeps/i386/fpu/e_expl.S (IEEE754_EXPL) [USE_AS_EXPM1L]: Just
> 	return the argument for normal arguments with exponent below -64.
> 	* sysdeps/x86_64/fpu/e_expl.S (IEEE754_EXPL) [USE_AS_EXPM1L]:
> 	Likewise.
> 	* math/auto-libm-test-in: Add another test of expm1.
> 	* math/auto-libm-test-out: Regenerated.
>
> diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
> index 4eaa013..0826825 100644
> --- a/math/auto-libm-test-in
> +++ b/math/auto-libm-test-in
> @@ -960,6 +960,8 @@ expm1 0x1p-64
>   expm1 -0x1p-64
>   expm1 0x1p-100
>   expm1 -0x1p-100
> +# Bug 16353: underflow exception may be missing
> +expm1 0x4.0000000000000028p-16384 missing-underflow
>
>   fma 1.0 2.0 3.0
>   fma 1.25 0.75 0.0625
> diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
> index eb4086b..c7e4373 100644
> --- a/sysdeps/i386/fpu/e_expl.S
> +++ b/sysdeps/i386/fpu/e_expl.S
> @@ -108,6 +108,16 @@ ENTRY(IEEE754_EXPL)
>   	andb	%ah, %dh
>   	cmpb	$0x40, %dh
>   	je	2f
> +
> +	/* Test for arguments that are small but not subnormal.  */
> +	movzwl	4+8(%esp), %eax
> +	andl	$0x7fff, %eax
> +	cmpl	$0x3fbf, %eax
> +	jge	3f
> +	/* Argument's exponent below -64; avoid spurious underflow if
> +	   normal.  */
> +	cmpl	$0x0001, %eax
> +	jge	2f
>   #else
>   	movzwl	4+8(%esp), %eax
>   	andl	$0x7fff, %eax
> diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S
> index 338527b..0ebe388 100644
> --- a/sysdeps/x86_64/fpu/e_expl.S
> +++ b/sysdeps/x86_64/fpu/e_expl.S
> @@ -105,6 +105,16 @@ ENTRY(IEEE754_EXPL)
>   	andb	%ah, %dh
>   	cmpb	$0x40, %dh
>   	je	2f
> +
> +	/* Test for arguments that are small but not subnormal.  */
> +	movzwl	8+8(%rsp), %eax
> +	andl	$0x7fff, %eax
> +	cmpl	$0x3fbf, %eax
> +	jge	3f
> +	/* Argument's exponent below -64; avoid spurious underflow if
> +	   normal.  */
> +	cmpl	$0x0001, %eax
> +	jge	2f
>   #else
>   	movzwl	8+8(%rsp), %eax
>   	andl	$0x7fff, %eax
>
  

Patch

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4eaa013..0826825 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -960,6 +960,8 @@  expm1 0x1p-64
 expm1 -0x1p-64
 expm1 0x1p-100
 expm1 -0x1p-100
+# Bug 16353: underflow exception may be missing
+expm1 0x4.0000000000000028p-16384 missing-underflow
 
 fma 1.0 2.0 3.0
 fma 1.25 0.75 0.0625
diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index eb4086b..c7e4373 100644
--- a/sysdeps/i386/fpu/e_expl.S
+++ b/sysdeps/i386/fpu/e_expl.S
@@ -108,6 +108,16 @@  ENTRY(IEEE754_EXPL)
 	andb	%ah, %dh
 	cmpb	$0x40, %dh
 	je	2f
+
+	/* Test for arguments that are small but not subnormal.  */
+	movzwl	4+8(%esp), %eax
+	andl	$0x7fff, %eax
+	cmpl	$0x3fbf, %eax
+	jge	3f
+	/* Argument's exponent below -64; avoid spurious underflow if
+	   normal.  */
+	cmpl	$0x0001, %eax
+	jge	2f
 #else
 	movzwl	4+8(%esp), %eax
 	andl	$0x7fff, %eax
diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S
index 338527b..0ebe388 100644
--- a/sysdeps/x86_64/fpu/e_expl.S
+++ b/sysdeps/x86_64/fpu/e_expl.S
@@ -105,6 +105,16 @@  ENTRY(IEEE754_EXPL)
 	andb	%ah, %dh
 	cmpb	$0x40, %dh
 	je	2f
+
+	/* Test for arguments that are small but not subnormal.  */
+	movzwl	8+8(%rsp), %eax
+	andl	$0x7fff, %eax
+	cmpl	$0x3fbf, %eax
+	jge	3f
+	/* Argument's exponent below -64; avoid spurious underflow if
+	   normal.  */
+	cmpl	$0x0001, %eax
+	jge	2f
 #else
 	movzwl	8+8(%rsp), %eax
 	andl	$0x7fff, %eax