[RFC] x86_64: 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.

This does not introduce any regressions other than the following:

testing long double (without inline functions)
Failure: Test: pown_upward (-0x1p+8192, 3LL)
Result:
 is:         -inf  -inf
 should be:  -1.18973149535723176502e+4932  -0xf.fffffffffffffff0p+16380
Failure: Test: pown_upward (-0x1p-8192, -3LL)
Result:
 is:         -inf  -inf
 should be:  -1.18973149535723176502e+4932  -0xf.fffffffffffffff0p+16380
Failure: Test: pown_upward (-0x2p-16384, -3LL)
Result:
 is:         -inf  -inf
 should be:  -1.18973149535723176502e+4932  -0xf.fffffffffffffff0p+16380
Failure: Test: pown_upward (-0x4p-16384, -3LL)
Result:
 is:         -inf  -inf
 should be:  -1.18973149535723176502e+4932  -0xf.fffffffffffffff0p+16380
Failure: Test: pown_upward (-0xf.fffffffffffffffp+16380, 3LL)
Result:
 is:         -inf  -inf
 should be:  -1.18973149535723176502e+4932  -0xf.fffffffffffffff0p+16380

The results seem OK to me and I'm inclined to update the results file,
but I'm not entirely sure, so I need someone to tell me I'm right, or
of course, wrong :)

Resolves: BZ #33411
Signed-off-by: Siddhesh Poyarekar <siddhesh@sourceware.org>
---

PS: I just noticed that powl does not actually have any test inputs and it's
likely only being tested via powr/pown (I assume? since they get affected by
changes to __ieee754_powl?).  That's why this doesn't have a test yet, I'll
post a proper patch with a test if this approach looks fine.

 sysdeps/x86_64/fpu/e_powl.S | 49 +++++++++++++++----------------------
 1 file changed, 20 insertions(+), 29 deletions(-)

diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 620ef765a7..35c50f37f3 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -144,39 +144,30 @@ 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.  */
-	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
+
+	/* y range is further reduced to [0, 3], so just run through the
+	   combination.  */
+	test	%eax, %eax
+	jz	6f
+
+5:	fmul	%st(1)		// RES*x : x
+	subl	$1, %eax
+	test	%eax, %eax
+	jz	6f
+	jmp	5b
+
+	/* Result is in %st(1), so pop the stack once.  */
+6:	fxch
+	fstp	%st(0)		// RES
 	LDBL_CHECK_FORCE_UFLOW_NONNAN
 	ret
 
-- 
2.51.0