[PATCH 2/2] Simplify powl computation for small integral y [BZ #33411]

[Siddhesh Poyarekar]

The powl implementation for x86_64 ends up multiplying X once more than
necessary and then throwing away that result.  This results in an
overflow flag being set in cases where there is no overflow.

Simplify the relevant portion by special casing the -3 to 3 range and
simply multiplying repetitively.

Resolves: BZ #33411
Signed-off-by: Siddhesh Poyarekar <siddhesh@sourceware.org>
---
 math/test-powl.c            | 22 +++++++++++--
 sysdeps/x86_64/fpu/e_powl.S | 61 +++++++++++++++++++++----------------
 2 files changed, 54 insertions(+), 29 deletions(-)

diff --git a/math/test-powl.c b/math/test-powl.c
index 2c716043b8..25253d51d7 100644
--- a/math/test-powl.c
+++ b/math/test-powl.c
@@ -16,10 +16,11 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
-#include <stdio.h>
-#include <math.h>
+#include <fenv.h>
 #include <float.h>
 #include <ieee754.h>
+#include <math.h>
+#include <stdio.h>
 #include <support/check.h>
 
 /* Smoke test for bz12775, which discovered a typo in the x86_64 powl
@@ -46,10 +47,27 @@ do_test_bz12775 (void)
 #endif
 }
 
+/* Make sure that a long double ^ 1 does not set the overflow flag.  For
+   context see BZ #33411.  */
+static void
+do_test_overflow (void)
+{
+    long double x = 1.090748135619415929463e+2466L;
+
+    feclearexcept( FE_ALL_EXCEPT );
+
+    long double res = powl (x, 1.L);
+    int flags = fetestexcept( FE_ALL_EXCEPT );
+
+    TEST_VERIFY (res == x);
+    TEST_COMPARE (flags & FE_OVERFLOW, 0);
+}
+
 static int
 do_test (void)
 {
   do_test_bz12775 ();
+  do_test_overflow ();
 
   return 0;
 }
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 620ef765a7..5695e1a9a6 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -144,40 +144,47 @@ ENTRY(__ieee754_powl)
 	fcomip	%st(1), %st	// 4 : y : x
 	fstp	%st(0)		// y : x
 	jnc	3f
-	mov	-8(%rsp),%eax
-	mov	-4(%rsp),%edx
-	orl	$0, %edx
+
+	/* Here onwards, it's just integral y in range [-3, 3].  */
+	movq	-8(%rsp),%rax
+	orq	$0, %rax
 	fstp	%st(0)		// x
 	jns	4f		// y >= 0, jump
 	fdivrl	MO(one)		// 1/x		(now referred to as x)
-	negl	%eax
-	adcl	$0, %edx
-	negl	%edx
+	negq	%rax
 4:	fldl	MO(one)		// 1 : x
-	fxch
 
-	/* If y is even, take the absolute value of x.  Otherwise,
-	   ensure all intermediate values that might overflow have the
-	   sign of x.  */
+	/* y range is further reduced to [0, 3].  Simply walk through the
+	   options.  */
+	test	%eax, %eax
+	jnz	5f
+	fstp	%st(1)		// 1
+	jmp	7f
+
+	/* y == 1.  */
+5:	movl	%eax, %ecx
+	xorl	$1, %ecx
+	fstp	%st(0)		// x
+	jz	7f
+
+	/* y == 2.  */
 	testb	$1, %al
 	jnz	6f
-	fabs
-
-6:	shrdl	$1, %edx, %eax
-	jnc	5f
-	fxch
-	fabs
-	fmul	%st(1)		// x : ST*x
-	fxch
-5:	fld	%st		// x : x : ST*x
-	fabs			// |x| : x : ST*x
-	fmulp			// |x|*x : ST*x
-	shrl	$1, %edx
-	movl	%eax, %ecx
-	orl	%edx, %ecx
-	jnz	6b
-	fstp	%st(0)		// ST*x
-	LDBL_CHECK_FORCE_UFLOW_NONNAN
+	fabs			// |x|
+	fld	%st		// |x| : |x|
+	fmulp	%st(1)		// |x| * |x|
+	jmp	7f
+
+	/* y == 3.  Special cased because x * x and |x| * |x| decay towards
+	   infinity faster compared to x * |x|.  */
+6:	fld	%st		// x : x
+	fabs			// |x| : x
+	fxch			// x : |x|
+	fld	%st(1)		// |x| : x : |x|
+	fmulp			// |x| * x : |x|
+	fmulp			// |x| * x * |x|
+
+7:	LDBL_CHECK_FORCE_UFLOW_NONNAN
 	ret
 
 	/* y is ±NAN */
-- 
2.51.0