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

[H.J. Lu]

On Thu, Oct 30, 2025 at 6:53 AM H.J. Lu <hjl.tools@gmail.com> wrote:
>
> 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.

I opened:

https://sourceware.org/bugzilla/show_bug.cgi?id=33586

-- 
H.J.