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

[H.J. Lu]

On Tue, Oct 21, 2025 at 8:22 PM Siddhesh Poyarekar
<siddhesh@sourceware.org> wrote:
>
> 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 | 56 +++++++++++++++++++------------------
>  2 files changed, 49 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..39f77480e8 100644
> --- a/sysdeps/x86_64/fpu/e_powl.S
> +++ b/sysdeps/x86_64/fpu/e_powl.S
> @@ -144,39 +144,41 @@ 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.  First up, 0 and 1.  */
> +       test    %eax, %eax
> +       jz      6f
> +       fxch                    // x : 1
> +       subl    $1, %eax
> +       jz      6f
> +
> +       /* Finally, y == 2 and 3.  For y == 3 we do |x| * x * |x| because x * x
> +          and |x| * |x| decay faster towards infinity compared to x * |x|.  */
> +       fld     %st             // x : x : 1
> +       fabs                    // |x| : x : 1
> +       fxch                    // x : |x| : 1
> +       fld     %st(1)          // |x| : x : |x| : 1
>         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
> +       jz      7f
> +       fmulp   %st(2)          // x : |x| * |x| : 1
> +       fstp    %st(0)          // |x| * |x| : 1
> +       jmp     6f
> +7:     fmulp                   // |x| * x : |x| : 1
> +       fmulp                   // |x| * x * |x| : 1
> +
> +       /* We come here with the stack as RES : <something>, so pop off
> +          <something>.  */
> +6:     fstp    %st(1)
>         LDBL_CHECK_FORCE_UFLOW_NONNAN
>         ret
>
> --
> 2.51.0
>

This caused:

FAIL: math/test-float64x-pow
original exit status 1
testing _Float64x (without inline functions)
Failure: pow (0x1p+8192, 0x1p+0): Exception "Overflow" set
Failure: pow_downward (0x1p+8192, 0x1p+0): Exception "Overflow" set
Failure: pow_towardzero (0x1p+8192, 0x1p+0): Exception "Overflow" set
Failure: pow_upward (0x1p+8192, 0x1p+0): Exception "Overflow" set

on i686 with -O2 -march=i686.

-- 
H.J.