diff mbox

Fix log2 (1) in round-downward mode (bug 17042)

Message ID Pine.LNX.4.64.1406092254280.4133@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers June 9, 2014, 10:55 p.m. UTC 
As with other issues of this kind, bug 17042 is log2 (1) wrongly
returning -0 instead of +0 in round-downward mode because of
implementations effectively in terms of log1p (x - 1).  This patch
fixes the issue in the same way used for log and log10.

Tested x86_64 and x86 and ulps updated accordingly.  Also tested for
mips64 to confirm a fix was needed for ldbl-128 and to validate that
fix (also applied to ldbl-128ibm since that version of log2l is
essentially the same as the ldbl-128 one).

(Other functions remaining to be converted to all-rounding-mode
testing with all inputs are pow
<https://sourceware.org/ml/libc-alpha/2014-06/msg00077.html>, cexp,
cpow, jn, tgamma, yn.)

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

	[BZ #17042]
	* sysdeps/i386/fpu/e_log2.S (__ieee754_log2): Take absolete value
	when x - 1 is zero.
	* sysdeps/i386/fpu/e_log2f.S (__ieee754_log2f): Likewise.
	* sysdeps/i386/fpu/e_log2l.S (__ieee754_log2l): Likewise.
	* sysdeps/ieee754/ldbl-128/e_log2l.c (__ieee754_log2l): Return
	0.0L for an argument of 1.0L.
	* sysdeps/ieee754/ldbl-128ibm/e_log2l.c (__ieee754_log2l):
	Likewise.
	* sysdeps/x86_64/fpu/e_log2l.S (__ieee754_log2l): Take absolute
	value when x - 1 is zero.
	* math/libm-test.inc (log2_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

Comments

Carlos O'Donell June 10, 2014, 3:14 a.m. UTC | #1 
On 06/09/2014 06:55 PM, Joseph S. Myers wrote:
> As with other issues of this kind, bug 17042 is log2 (1) wrongly
> returning -0 instead of +0 in round-downward mode because of
> implementations effectively in terms of log1p (x - 1).  This patch
> fixes the issue in the same way used for log and log10.
> 
> Tested x86_64 and x86 and ulps updated accordingly.  Also tested for
> mips64 to confirm a fix was needed for ldbl-128 and to validate that
> fix (also applied to ldbl-128ibm since that version of log2l is
> essentially the same as the ldbl-128 one).
> 
> (Other functions remaining to be converted to all-rounding-mode
> testing with all inputs are pow
> <https://sourceware.org/ml/libc-alpha/2014-06/msg00077.html>, cexp,
> cpow, jn, tgamma, yn.)
> 
> 2014-06-09  Joseph Myers  <joseph@codesourcery.com>
> 
> 	[BZ #17042]
> 	* sysdeps/i386/fpu/e_log2.S (__ieee754_log2): Take absolete value
> 	when x - 1 is zero.
> 	* sysdeps/i386/fpu/e_log2f.S (__ieee754_log2f): Likewise.
> 	* sysdeps/i386/fpu/e_log2l.S (__ieee754_log2l): Likewise.
> 	* sysdeps/ieee754/ldbl-128/e_log2l.c (__ieee754_log2l): Return
> 	0.0L for an argument of 1.0L.
> 	* sysdeps/ieee754/ldbl-128ibm/e_log2l.c (__ieee754_log2l):
> 	Likewise.
> 	* sysdeps/x86_64/fpu/e_log2l.S (__ieee754_log2l): Take absolute
> 	value when x - 1 is zero.
> 	* math/libm-test.inc (log2_test): Use ALL_RM_TEST.
> 	* sysdeps/i386/fpu/libm-test-ulps: Update.
> 	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

OK for me. 

My preference is to avoid the magic numbers so I don't have to
keep looking up the constants for the status word or provide a
comment.

Any idea how the assembly variants perform against the C-variants?

> diff --git a/math/libm-test.inc b/math/libm-test.inc
> index 0d467a2..fa8e238 100644
> --- a/math/libm-test.inc
> +++ b/math/libm-test.inc
> @@ -7840,9 +7840,7 @@ static const struct test_f_f_data log2_test_data[] =
>  static void
>  log2_test (void)
>  {
> -  START (log2, 0);
> -  RUN_TEST_LOOP_f_f (log2, log2_test_data, );
> -  END;
> +  ALL_RM_TEST (log2, 0, log2_test_data, RUN_TEST_LOOP_f_f, END);

OK.

>  }
>  
>  
> diff --git a/sysdeps/i386/fpu/e_log2.S b/sysdeps/i386/fpu/e_log2.S
> index a202bc7..73ff0ff 100644
> --- a/sysdeps/i386/fpu/e_log2.S
> +++ b/sysdeps/i386/fpu/e_log2.S
> @@ -47,7 +47,13 @@ ENTRY(__ieee754_log2)
>  	fnstsw			// x-1 : x : 1
>  	andb	$0x45, %ah
>  	jz	2f
> -	fstp	%st(1)		// x-1 : 1
> +	fxam
> +	fnstsw
> +	andb	$0x45, %ah
> +	cmpb	$0x40, %ah
> +	jne	5f
> +	fabs			// log2(1) is +0 in all rounding modes.
> +5:	fstp	%st(1)		// x-1 : 1

OK.

>  	fyl2xp1			// log(x)
>  	ret
>  
> diff --git a/sysdeps/i386/fpu/e_log2f.S b/sysdeps/i386/fpu/e_log2f.S
> index f4f9a8c..344eeb4 100644
> --- a/sysdeps/i386/fpu/e_log2f.S
> +++ b/sysdeps/i386/fpu/e_log2f.S
> @@ -47,7 +47,13 @@ ENTRY(__ieee754_log2f)
>  	fnstsw			// x-1 : x : 1
>  	andb	$0x45, %ah
>  	jz	2f
> -	fstp	%st(1)		// x-1 : 1
> +	fxam
> +	fnstsw
> +	andb	$0x45, %ah
> +	cmpb	$0x40, %ah
> +	jne	5f
> +	fabs			// log2(1) is +0 in all rounding modes.
> +5:	fstp	%st(1)		// x-1 : 1

OK.

>  	fyl2xp1			// log(x)
>  	ret
>  
> diff --git a/sysdeps/i386/fpu/e_log2l.S b/sysdeps/i386/fpu/e_log2l.S
> index bd51b56..0f5f7e5 100644
> --- a/sysdeps/i386/fpu/e_log2l.S
> +++ b/sysdeps/i386/fpu/e_log2l.S
> @@ -47,7 +47,13 @@ ENTRY(__ieee754_log2l)
>  	fnstsw			// x-1 : x : 1
>  	andb	$0x45, %ah
>  	jz	2f
> -	fstp	%st(1)		// x-1 : 1
> +	fxam
> +	fnstsw
> +	andb	$0x45, %ah
> +	cmpb	$0x40, %ah
> +	jne	5f
> +	fabs			// log2(1) is +0 in all rounding modes.
> +5:	fstp	%st(1)		// x-1 : 1

OK.

>  	fyl2xp1			// log(x)
>  	ret
>  
> diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
> index 1e89284..d7424a6 100644
> --- a/sysdeps/i386/fpu/libm-test-ulps
> +++ b/sysdeps/i386/fpu/libm-test-ulps
> @@ -1588,6 +1588,22 @@ ifloat: 1
>  ildouble: 1
>  ldouble: 1
>  
> +Function: "log2_towardzero":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +
> +Function: "log2_upward":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +
>  Function: "log_downward":
>  ildouble: 1
>  ldouble: 1
> diff --git a/sysdeps/ieee754/ldbl-128/e_log2l.c b/sysdeps/ieee754/ldbl-128/e_log2l.c
> index 6c7da0e..991a3b7 100644
> --- a/sysdeps/ieee754/ldbl-128/e_log2l.c
> +++ b/sysdeps/ieee754/ldbl-128/e_log2l.c
> @@ -188,6 +188,9 @@ __ieee754_log2l (x)
>    if (hx >= 0x7fff000000000000LL)
>      return (x + x);
>  
> +  if (x == 1.0L)
> +    return 0.0L;

OK.

> +
>  /* separate mantissa from exponent */
>  
>  /* Note, frexp is used so that denormal numbers
> diff --git a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
> index 323ded0..442ad97 100644
> --- a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
> +++ b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
> @@ -190,6 +190,9 @@ __ieee754_log2l (x)
>    if (hx >= 0x7ff0000000000000LL)
>      return (x + x);
>  
> +  if (x == 1.0L)
> +    return 0.0L;

OK.

> +
>  /* separate mantissa from exponent */
>  
>  /* Note, frexp is used so that denormal numbers
> diff --git a/sysdeps/x86_64/fpu/e_log2l.S b/sysdeps/x86_64/fpu/e_log2l.S
> index 956489f..c12906d 100644
> --- a/sysdeps/x86_64/fpu/e_log2l.S
> +++ b/sysdeps/x86_64/fpu/e_log2l.S
> @@ -45,7 +45,13 @@ ENTRY(__ieee754_log2l)
>  	fnstsw			// x-1 : x : 1
>  	andb	$0x45, %ah
>  	jz	2f
> -	fstp	%st(1)		// x-1 : 1
> +	fxam
> +	fnstsw
> +	andb	$0x45, %ah
> +	cmpb	$0x40, %ah
> +	jne	5f
> +	fabs			// log2(1) is +0 in all rounding modes.
> +5:	fstp	%st(1)		// x-1 : 1

OK.

>  	fyl2xp1			// log(x)
>  	ret
>  
> diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
> index bb549d2..92fef5a 100644
> --- a/sysdeps/x86_64/fpu/libm-test-ulps
> +++ b/sysdeps/x86_64/fpu/libm-test-ulps
> @@ -1665,6 +1665,28 @@ ifloat: 1
>  ildouble: 1
>  ldouble: 1
>  
> +Function: "log2_downward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +
> +Function: "log2_towardzero":
> +double: 1
> +float: 1
> +idouble: 1
> +ifloat: 1
> +ildouble: 1
> +ldouble: 1
> +
> +Function: "log2_upward":
> +double: 2
> +float: 2
> +idouble: 2
> +ifloat: 2
> +ildouble: 1
> +ldouble: 1
> +
>  Function: "log_downward":
>  float: 1
>  ifloat: 1
> 

Cheers,
Carlos.
Joseph Myers June 10, 2014, 12:01 p.m. UTC | #2 
On Mon, 9 Jun 2014, Carlos O'Donell wrote:

> Any idea how the assembly variants perform against the C-variants?

No idea.  We don't have C variants for log2l for ldbl-96, just separate 
assembly variants for x86 / x86_64 / ia64 / m68k (the last via inline 
asm).
diff mbox

Patch

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0d467a2..fa8e238 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -7840,9 +7840,7 @@  static const struct test_f_f_data log2_test_data[] =
 static void
 log2_test (void)
 {
-  START (log2, 0);
-  RUN_TEST_LOOP_f_f (log2, log2_test_data, );
-  END;
+  ALL_RM_TEST (log2, 0, log2_test_data, RUN_TEST_LOOP_f_f, END);
 }
 
 
diff --git a/sysdeps/i386/fpu/e_log2.S b/sysdeps/i386/fpu/e_log2.S
index a202bc7..73ff0ff 100644
--- a/sysdeps/i386/fpu/e_log2.S
+++ b/sysdeps/i386/fpu/e_log2.S
@@ -47,7 +47,13 @@  ENTRY(__ieee754_log2)
 	fnstsw			// x-1 : x : 1
 	andb	$0x45, %ah
 	jz	2f
-	fstp	%st(1)		// x-1 : 1
+	fxam
+	fnstsw
+	andb	$0x45, %ah
+	cmpb	$0x40, %ah
+	jne	5f
+	fabs			// log2(1) is +0 in all rounding modes.
+5:	fstp	%st(1)		// x-1 : 1
 	fyl2xp1			// log(x)
 	ret
 
diff --git a/sysdeps/i386/fpu/e_log2f.S b/sysdeps/i386/fpu/e_log2f.S
index f4f9a8c..344eeb4 100644
--- a/sysdeps/i386/fpu/e_log2f.S
+++ b/sysdeps/i386/fpu/e_log2f.S
@@ -47,7 +47,13 @@  ENTRY(__ieee754_log2f)
 	fnstsw			// x-1 : x : 1
 	andb	$0x45, %ah
 	jz	2f
-	fstp	%st(1)		// x-1 : 1
+	fxam
+	fnstsw
+	andb	$0x45, %ah
+	cmpb	$0x40, %ah
+	jne	5f
+	fabs			// log2(1) is +0 in all rounding modes.
+5:	fstp	%st(1)		// x-1 : 1
 	fyl2xp1			// log(x)
 	ret
 
diff --git a/sysdeps/i386/fpu/e_log2l.S b/sysdeps/i386/fpu/e_log2l.S
index bd51b56..0f5f7e5 100644
--- a/sysdeps/i386/fpu/e_log2l.S
+++ b/sysdeps/i386/fpu/e_log2l.S
@@ -47,7 +47,13 @@  ENTRY(__ieee754_log2l)
 	fnstsw			// x-1 : x : 1
 	andb	$0x45, %ah
 	jz	2f
-	fstp	%st(1)		// x-1 : 1
+	fxam
+	fnstsw
+	andb	$0x45, %ah
+	cmpb	$0x40, %ah
+	jne	5f
+	fabs			// log2(1) is +0 in all rounding modes.
+5:	fstp	%st(1)		// x-1 : 1
 	fyl2xp1			// log(x)
 	ret
 
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 1e89284..d7424a6 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1588,6 +1588,22 @@  ifloat: 1
 ildouble: 1
 ldouble: 1
 
+Function: "log2_towardzero":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "log2_upward":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
 Function: "log_downward":
 ildouble: 1
 ldouble: 1
diff --git a/sysdeps/ieee754/ldbl-128/e_log2l.c b/sysdeps/ieee754/ldbl-128/e_log2l.c
index 6c7da0e..991a3b7 100644
--- a/sysdeps/ieee754/ldbl-128/e_log2l.c
+++ b/sysdeps/ieee754/ldbl-128/e_log2l.c
@@ -188,6 +188,9 @@  __ieee754_log2l (x)
   if (hx >= 0x7fff000000000000LL)
     return (x + x);
 
+  if (x == 1.0L)
+    return 0.0L;
+
 /* separate mantissa from exponent */
 
 /* Note, frexp is used so that denormal numbers
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
index 323ded0..442ad97 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_log2l.c
@@ -190,6 +190,9 @@  __ieee754_log2l (x)
   if (hx >= 0x7ff0000000000000LL)
     return (x + x);
 
+  if (x == 1.0L)
+    return 0.0L;
+
 /* separate mantissa from exponent */
 
 /* Note, frexp is used so that denormal numbers
diff --git a/sysdeps/x86_64/fpu/e_log2l.S b/sysdeps/x86_64/fpu/e_log2l.S
index 956489f..c12906d 100644
--- a/sysdeps/x86_64/fpu/e_log2l.S
+++ b/sysdeps/x86_64/fpu/e_log2l.S
@@ -45,7 +45,13 @@  ENTRY(__ieee754_log2l)
 	fnstsw			// x-1 : x : 1
 	andb	$0x45, %ah
 	jz	2f
-	fstp	%st(1)		// x-1 : 1
+	fxam
+	fnstsw
+	andb	$0x45, %ah
+	cmpb	$0x40, %ah
+	jne	5f
+	fabs			// log2(1) is +0 in all rounding modes.
+5:	fstp	%st(1)		// x-1 : 1
 	fyl2xp1			// log(x)
 	ret
 
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index bb549d2..92fef5a 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1665,6 +1665,28 @@  ifloat: 1
 ildouble: 1
 ldouble: 1
 
+Function: "log2_downward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+
+Function: "log2_towardzero":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "log2_upward":
+double: 2
+float: 2
+idouble: 2
+ifloat: 2
+ildouble: 1
+ldouble: 1
+
 Function: "log_downward":
 float: 1
 ifloat: 1