[glibc] AArch64: Implement AdvSIMD and SVE powr(f) routines
Adhemerval Zanella via Glibc-cvs <[email protected]>
| Newsgroups | gmane.comp.lib.glibc.cvs |
|---|---|
| Message-ID | <[email protected]> |
https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=2ff2565df78160b392b120793d6d1aac2d716701 commit 2ff2565df78160b392b120793d6d1aac2d716701 Author: Pierre Blanchard <[email protected]> Date: Wed Apr 15 08:32:44 2026 +0000 AArch64: Implement AdvSIMD and SVE powr(f) routines Vector variants of the new C23 powr routines. These provide same maximum error error as pow by virtue of relying on shared approximation techniques and sources. Note: Benchmark inputs for powr(f) are identical to pow(f). Performance gain over pow on V1 with GCC@15: - SVE powr: 10-12% on subnormal x, 12-13% on x < 0. - SVE powrf: 15% on all x < 0. - AdvSIMD powr: for x < 0, 40% if x subnormal, 60% otherwise. - AdvSIMD powrf: 4% on x subnormals or x < 0. Diff: --- bits/libm-simd-decl-stubs.h | 11 ++ math/bits/mathcalls.h | 1 + sysdeps/aarch64/fpu/Makefile | 1 + sysdeps/aarch64/fpu/Versions | 7 + sysdeps/aarch64/fpu/advsimd_f32_protos.h | 1 + sysdeps/aarch64/fpu/bits/math-vector.h | 8 ++ sysdeps/aarch64/fpu/finclude/math-vector-fortran.h | 2 + sysdeps/aarch64/fpu/powr_advsimd.c | 148 +++++++++++++++++++++ sysdeps/aarch64/fpu/powr_sve.c | 123 +++++++++++++++++ sysdeps/aarch64/fpu/powrf_advsimd.c | 135 +++++++++++++++++++ sysdeps/aarch64/fpu/powrf_sve.c | 135 +++++++++++++++++++ sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c | 1 + sysdeps/aarch64/fpu/test-double-sve-wrappers.c | 1 + sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c | 1 + sysdeps/aarch64/fpu/test-float-sve-wrappers.c | 1 + sysdeps/aarch64/fpu/v_powrf_inline.h | 4 +- sysdeps/unix/sysv/linux/aarch64/libmvec.abilist | 5 + 17 files changed, 584 insertions(+), 1 deletion(-) diff --git a/bits/libm-simd-decl-stubs.h b/bits/libm-simd-decl-stubs.h index 2b19901d3e..5cb0e4a245 100644 --- a/bits/libm-simd-decl-stubs.h +++ b/bits/libm-simd-decl-stubs.h @@ -99,6 +99,17 @@ #define __DECL_SIMD_powf64x #define __DECL_SIMD_powf128x +#define __DECL_SIMD_powr +#define __DECL_SIMD_powrf +#define __DECL_SIMD_powrl +#define __DECL_SIMD_powrf16 +#define __DECL_SIMD_powrf32 +#define __DECL_SIMD_powrf64 +#define __DECL_SIMD_powrf128 +#define __DECL_SIMD_powrf32x +#define __DECL_SIMD_powrf64x +#define __DECL_SIMD_powrf128x + #define __DECL_SIMD_acos #define __DECL_SIMD_acosf #define __DECL_SIMD_acosl diff --git a/math/bits/mathcalls.h b/math/bits/mathcalls.h index a4c994d0a7..4e983f5c7b 100644 --- a/math/bits/mathcalls.h +++ b/math/bits/mathcalls.h @@ -197,6 +197,7 @@ __MATHCALL (compoundn,, (_Mdouble_ __x, long long int __y)); __MATHCALL (pown,, (_Mdouble_ __x, long long int __y)); /* Return X to the Y power. */ +__MATHCALL_VEC (powr,, (_Mdouble_ __x, _Mdouble_ __y)); __MATHCALL (powr,, (_Mdouble_ __x, _Mdouble_ __y)); /* Return the Yth root of X. */ diff --git a/sysdeps/aarch64/fpu/Makefile b/sysdeps/aarch64/fpu/Makefile index 998fc08d43..6c8cacf21d 100644 --- a/sysdeps/aarch64/fpu/Makefile +++ b/sysdeps/aarch64/fpu/Makefile @@ -29,6 +29,7 @@ libmvec-supported-funcs = acos \ log2 \ log2p1 \ pow \ + powr \ rsqrt \ sin \ sinh \ diff --git a/sysdeps/aarch64/fpu/Versions b/sysdeps/aarch64/fpu/Versions index d68510a20e..2337f9d331 100644 --- a/sysdeps/aarch64/fpu/Versions +++ b/sysdeps/aarch64/fpu/Versions @@ -206,4 +206,11 @@ libmvec { _ZGVsMxv_rsqrt; _ZGVsMxv_rsqrtf; } + GLIBC_2.44 { + _ZGVnN2vv_powr; + _ZGVnN2vv_powrf; + _ZGVnN4vv_powrf; + _ZGVsMxvv_powr; + _ZGVsMxvv_powrf; + } } diff --git a/sysdeps/aarch64/fpu/advsimd_f32_protos.h b/sysdeps/aarch64/fpu/advsimd_f32_protos.h index 81de7351f1..59210b11ad 100644 --- a/sysdeps/aarch64/fpu/advsimd_f32_protos.h +++ b/sysdeps/aarch64/fpu/advsimd_f32_protos.h @@ -47,6 +47,7 @@ libmvec_hidden_proto (V_NAME_F1(log2p1)); libmvec_hidden_proto (V_NAME_F1(logp1)); libmvec_hidden_proto (V_NAME_F1(log)); libmvec_hidden_proto (V_NAME_F2(pow)); +libmvec_hidden_proto (V_NAME_F2(powr)); libmvec_hidden_proto (V_NAME_F1(rsqrt)); libmvec_hidden_proto (V_NAME_F1(sin)); libmvec_hidden_proto (V_NAME_F1(sinh)); diff --git a/sysdeps/aarch64/fpu/bits/math-vector.h b/sysdeps/aarch64/fpu/bits/math-vector.h index 442cd2a02c..db218eedf9 100644 --- a/sysdeps/aarch64/fpu/bits/math-vector.h +++ b/sysdeps/aarch64/fpu/bits/math-vector.h @@ -157,6 +157,10 @@ # define __DECL_SIMD_pow __DECL_SIMD_aarch64 # undef __DECL_SIMD_powf # define __DECL_SIMD_powf __DECL_SIMD_aarch64 +# undef __DECL_SIMD_powr +# define __DECL_SIMD_powr __DECL_SIMD_aarch64 +# undef __DECL_SIMD_powrf +# define __DECL_SIMD_powrf __DECL_SIMD_aarch64 # undef __DECL_SIMD_rsqrt # define __DECL_SIMD_rsqrt __DECL_SIMD_aarch64 # undef __DECL_SIMD_rsqrtf @@ -243,6 +247,7 @@ __vpcs __f32x4_t _ZGVnN4v_log2f (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_log2p1f (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_logp1f (__f32x4_t); __vpcs __f32x4_t _ZGVnN4vv_powf (__f32x4_t, __f32x4_t); +__vpcs __f32x4_t _ZGVnN4vv_powrf (__f32x4_t, __f32x4_t); __vpcs __f32x4_t _ZGVnN4v_rsqrtf (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_sinf (__f32x4_t); __vpcs __f32x4_t _ZGVnN4v_sinhf (__f32x4_t); @@ -283,6 +288,7 @@ __vpcs __f64x2_t _ZGVnN2v_log2 (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_log2p1 (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_logp1 (__f64x2_t); __vpcs __f64x2_t _ZGVnN2vv_pow (__f64x2_t, __f64x2_t); +__vpcs __f64x2_t _ZGVnN2vv_powr (__f64x2_t, __f64x2_t); __vpcs __f64x2_t _ZGVnN2v_rsqrt (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_sin (__f64x2_t); __vpcs __f64x2_t _ZGVnN2v_sinh (__f64x2_t); @@ -328,6 +334,7 @@ __sv_f32_t _ZGVsMxv_log2f (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_log2p1f (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_logp1f (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxvv_powf (__sv_f32_t, __sv_f32_t, __sv_bool_t); +__sv_f32_t _ZGVsMxvv_powrf (__sv_f32_t, __sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_rsqrtf (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_sinf (__sv_f32_t, __sv_bool_t); __sv_f32_t _ZGVsMxv_sinhf (__sv_f32_t, __sv_bool_t); @@ -368,6 +375,7 @@ __sv_f64_t _ZGVsMxv_log2 (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_log2p1 (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_logp1 (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxvv_pow (__sv_f64_t, __sv_f64_t, __sv_bool_t); +__sv_f64_t _ZGVsMxvv_powr (__sv_f64_t, __sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_rsqrt (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_sin (__sv_f64_t, __sv_bool_t); __sv_f64_t _ZGVsMxv_sinh (__sv_f64_t, __sv_bool_t); diff --git a/sysdeps/aarch64/fpu/finclude/math-vector-fortran.h b/sysdeps/aarch64/fpu/finclude/math-vector-fortran.h index 46fc8a627c..71ee5d6a0e 100644 --- a/sysdeps/aarch64/fpu/finclude/math-vector-fortran.h +++ b/sysdeps/aarch64/fpu/finclude/math-vector-fortran.h @@ -80,6 +80,8 @@ !GCC$ builtin (logp1f) attributes simd (notinbranch) if('fastmath') !GCC$ builtin (pow) attributes simd (notinbranch) if('fastmath') !GCC$ builtin (powf) attributes simd (notinbranch) if('fastmath') +!GCC$ builtin (powr) attributes simd (notinbranch) if('fastmath') +!GCC$ builtin (powrf) attributes simd (notinbranch) if('fastmath') !GCC$ builtin (rsqrt) attributes simd (notinbranch) if('fastmath') !GCC$ builtin (rsqrtf) attributes simd (notinbranch) if('fastmath') !GCC$ builtin (sin) attributes simd (notinbranch) if('fastmath') diff --git a/sysdeps/aarch64/fpu/powr_advsimd.c b/sysdeps/aarch64/fpu/powr_advsimd.c new file mode 100644 index 0000000000..8163ae87ad --- /dev/null +++ b/sysdeps/aarch64/fpu/powr_advsimd.c @@ -0,0 +1,148 @@ +/* Double-precision vector (AdvSIMD) powr function + + Copyright (C) 2026 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#include "pow_common.h" +#include "v_math.h" + +#include "v_pow_inline.h" + +static double NOINLINE +powr_scalar_special_case (double x, double y) +{ + /* Negative x returns NaN (+0/-0 and NaN x not handled here). */ + if (x < 0) + return __builtin_nan (""); + + uint64_t ix = asuint64 (x); + uint64_t iy = asuint64 (y); + uint32_t topx = top12 (x); + uint32_t topy = top12 (y); + + /* Special cases: (x < 0x1p-126 or inf or nan) or + (|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */ + if (__glibc_unlikely (topx - SmallPowX >= ThresPowX + || (topy & 0x7ff) - SmallPowY >= ThresPowY)) + { + /* |y| is 0, Inf or NaN. */ + if (__glibc_unlikely (zeroinfnan (iy))) + { + if (2 * ix > 2 * asuint64 (INFINITY) + || 2 * iy > 2 * asuint64 (INFINITY)) + return __builtin_nan (""); + if (2 * iy == 0) + { + /* |x| = 0 or inf. */ + if ((2 * ix == 0) || (2 * ix == 2 * asuint64 (INFINITY))) + return __builtin_nan (""); + /* x is finite. */ + return 1.0; + } + /* |y| = Inf and x = 1.0. */ + if (ix == asuint64 (1.0)) + return __builtin_nan (""); + /* |x| < 1 and y = Inf or |x| > 1 and y = -Inf. */ + if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63)) + return 0.0; + /* |y| = Inf and previous conditions not met. */ + return y * y; + } + /* |x| is 0, Inf or NaN. */ + if (__glibc_unlikely (zeroinfnan (ix))) + { + double x2 = x * x; + return iy >> 63 ? 1 / x2 : x2; + } + /* Here x and y are non-zero finite. */ + /* Note: if |y| > 1075 * ln2 * 2^53 ~= 0x1.749p62 then powr(x,y) = inf/0 + and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then powr(x,y) = +-1. */ + if ((topy & 0x7ff) - SmallPowY >= ThresPowY) + { + if (ix == asuint64 (1.0)) + return 1.0; + /* |y| < 2^-65, x^y ~= 1 + y*log(x). */ + if ((topy & 0x7ff) < SmallPowY) + return 1.0; + return (ix > asuint64 (1.0)) == (topy < 0x800) ? INFINITY : 0; + } + if (topx == 0) + { + /* Normalize subnormal x so exponent becomes negative. */ + ix = asuint64 (x * 0x1p52); + ix -= 52ULL << 52; + } + } + + /* Core computation of exp (y * log (x)). */ + double lo; + double hi = log_inline (ix, &lo); + double ehi = y * hi; + double elo = y * lo + fma (y, hi, -ehi); + return exp_inline (ehi, elo, 0); +} + +static float64x2_t VPCS_ATTR NOINLINE +scalar_fallback (float64x2_t x, float64x2_t y) +{ + return (float64x2_t){ powr_scalar_special_case (x[0], y[0]), + powr_scalar_special_case (x[1], y[1]) }; +} + +/* Implementation of AdvSIMD powr. + Maximum measured error is 1.04 ULPs: + _ZGVnN2vv_powr(0x1.024a3e56b3c3p-136, 0x1.87910248b58acp-13) + got 0x1.f71162f473251p-1 + want 0x1.f71162f473252p-1. */ +float64x2_t VPCS_ATTR V_NAME_D2 (powr) (float64x2_t x, float64x2_t y) +{ + const struct data *d = ptr_barrier (&data); + + /* Case of x <= 0 is too complicated to be vectorised efficiently here, + fallback to scalar pow for all lanes if any x < 0 detected. */ + if (v_any_u64 (vclezq_s64 (vreinterpretq_s64_f64 (x)))) + return scalar_fallback (x, y); + + uint64x2_t vix = vreinterpretq_u64_f64 (x); + uint64x2_t viy = vreinterpretq_u64_f64 (y); + + /* Special cases of x or y. + The case y==0 does not trigger a special case, since in this case it is + necessary to fix the result only if x is a signalling nan, which already + triggers a special case. We test y==0 directly in the scalar fallback. */ + uint64x2_t x_is_inf_or_nan = vcgeq_u64 (vandq_u64 (vix, d->inf), d->inf); + uint64x2_t y_is_inf_or_nan = vcgeq_u64 (vandq_u64 (viy, d->inf), d->inf); + uint64x2_t special = vorrq_u64 (x_is_inf_or_nan, y_is_inf_or_nan); + + /* Fallback to scalar on all lanes if any lane is inf or nan. */ + if (__glibc_unlikely (v_any_u64 (special))) + return scalar_fallback (x, y); + + /* Cases of subnormal x: |x| < 0x1p-1022. */ + uint64x2_t x_is_subnormal = vcaltq_f64 (x, d->subnormal_bound); + if (__glibc_unlikely (v_any_u64 (x_is_subnormal))) + { + /* Normalize subnormal x so exponent becomes negative. */ + uint64x2_t vix_norm + = vreinterpretq_u64_f64 (vmulq_f64 (x, d->subnormal_scale)); + vix_norm = vsubq_u64 (vix_norm, d->subnormal_bias); + x = vbslq_f64 (x_is_subnormal, vreinterpretq_f64_u64 (vix_norm), x); + } + + /* Core computation of exp (y * log (x)). */ + return v_pow_inline (x, y, d); +} diff --git a/sysdeps/aarch64/fpu/powr_sve.c b/sysdeps/aarch64/fpu/powr_sve.c new file mode 100644 index 0000000000..ae599a037d --- /dev/null +++ b/sysdeps/aarch64/fpu/powr_sve.c @@ -0,0 +1,123 @@ +/* Double-precision vector (SVE) powr function + + Copyright (C) 2026 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#include "math_config.h" +#include "pow_common.h" +#include "sv_math.h" + +#define WANT_SV_POW_SIGN_BIAS 0 +#include "sv_pow_inline.h" + +/* A scalar subroutine used to fix main powr special cases. */ +static inline double +powr_specialcase (double x, double y) +{ + uint64_t ix = asuint64 (x); + uint64_t iy = asuint64 (y); + /* |y| is 0, Inf or NaN. */ + if (__glibc_unlikely (zeroinfnan (iy))) + { + /* |x| or |y| is NaN. */ + if (2 * ix > 2 * asuint64 (INFINITY) || 2 * iy > 2 * asuint64 (INFINITY)) + return __builtin_nan (""); + /* |y| is 0.0. */ + if (2 * iy == 0) + { + /* |x| = 0 or Inf. */ + if ((2 * ix == 0) || (2 * ix == 2 * asuint64 (INFINITY))) + return __builtin_nan (""); + /* x is finite. */ + return 1.0; + } + /* x is 1.0. */ + if (ix == asuint64 (1.0)) + return __builtin_nan (""); + /* |x| < 1 and y = Inf or |x| > 1 and y = -Inf. */ + if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63)) + return 0.0; + /* |y| = Inf and previous conditions not met. */ + return y * y; + } + /* x is 0, Inf or NaN. Negative x are handled in the core. */ + if (__glibc_unlikely (zeroinfnan (ix))) + { + double x2 = x * x; + return (iy >> 63) ? 1 / x2 : x2; + } + /* Return x for convenience, but make sure result is never used. */ + return x; +} + +/* Scalar fallback for special case routines with custom signature. */ +static svfloat64_t NOINLINE +sv_powr_specialcase (svfloat64_t x1, svfloat64_t x2, svfloat64_t y, + svbool_t cmp) +{ + return sv_call2_f64 (powr_specialcase, x1, x2, y, cmp); +} + +/* Implementation of SVE powr. + + Provides the same accuracy as AdvSIMD pow and powr, since it relies on the + same algorithm. + + Maximum measured error is 1.04 ULPs: + SV_NAME_D2 (powr) (0x1.3d2d45bc848acp+63, -0x1.a48a38b40cd43p-12) + got 0x1.f7116284221fcp-1 + want 0x1.f7116284221fdp-1. */ +svfloat64_t SV_NAME_D2 (powr) (svfloat64_t x, svfloat64_t y, const svbool_t pg) +{ + const struct data *d = ptr_barrier (&data); + + svuint64_t vix = svreinterpret_u64 (x); + svuint64_t viy = svreinterpret_u64 (y); + + svbool_t xpos = svcmpge (pg, x, sv_f64 (0.0)); + + /* Special cases of x or y: zero, inf and nan. */ + svbool_t xspecial = sv_zeroinfnan (xpos, vix); + svbool_t yspecial = sv_zeroinfnan (xpos, viy); + svbool_t cmp = svorr_z (xpos, xspecial, yspecial); + + /* Cases of positive subnormal x: 0 < x < 0x1p-1022. */ + svbool_t x_is_subnormal = svaclt (xpos, x, 0x1p-1022); + if (__glibc_unlikely (svptest_any (xpos, x_is_subnormal))) + { + /* Normalize subnormal x so exponent becomes negative. */ + svuint64_t vix_norm + = svreinterpret_u64 (svmul_m (x_is_subnormal, x, 0x1p52)); + vix = svsub_m (x_is_subnormal, vix_norm, 52ULL << 52); + } + + svfloat64_t vlo; + svfloat64_t vhi = sv_log_inline (xpos, vix, &vlo, d); + + svfloat64_t vehi = svmul_x (svptrue_b64 (), y, vhi); + svfloat64_t vemi = svmls_x (xpos, vehi, y, vhi); + svfloat64_t velo = svnmls_x (xpos, vemi, y, vlo); + svfloat64_t vz = sv_exp_inline (xpos, vehi, velo, sv_u64 (0), d); + + /* Cases of negative x. */ + vz = svsel (xpos, vz, sv_f64 (__builtin_nan (""))); + + if (__glibc_unlikely (svptest_any (cmp, cmp))) + return sv_powr_specialcase (x, y, vz, cmp); + + return vz; +} diff --git a/sysdeps/aarch64/fpu/powrf_advsimd.c b/sysdeps/aarch64/fpu/powrf_advsimd.c new file mode 100644 index 0000000000..c0d5de3d07 --- /dev/null +++ b/sysdeps/aarch64/fpu/powrf_advsimd.c @@ -0,0 +1,135 @@ +/* Single-precision vector (AdvSIMD) powr function + + Copyright (C) 2026 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#include "flt-32/math_config.h" +#include "v_math.h" +#include "v_powrf_inline.h" + +/* A scalar subroutine used to fix main powrf special cases. */ +static inline float +powrf_specialcase (float x, float y) +{ + /* Negative x returns NaN (+0/-0 and NaN x not handled here). */ + if (x < 0) + return __builtin_nanf (""); + + uint32_t ix = asuint (x); + uint32_t iy = asuint (y); + /* y is 0, Inf or NaN. */ + if (__glibc_unlikely (zeroinfnan (iy))) + { + /* |x| or |y| is NaN. */ + if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000) + return __builtin_nanf (""); + /* |y| = 0. */ + if (2 * iy == 0) + { + /* |x| = 0 or inf. */ + if ((2 * ix == 0) || (2 * ix == 2u * 0x7f800000)) + return __builtin_nanf (""); + /* x is finite. */ + return 1.0f; + } + /* |y| = Inf and x = 1.0. */ + if (ix == 0x3f800000) + return __builtin_nanf (""); + /* |x| < 1 and y = Inf or |x| > 1 and y = -Inf. */ + if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) + return 0.0f; + /* |y| = Inf and previous conditions not met. */ + return y * y; + } + /* x is 0, Inf or NaN. Negative x are handled in the core. */ + if (__glibc_unlikely (zeroinfnan (ix))) + { + float x2 = x * x; + return iy & 0x80000000 ? 1 / x2 : x2; + } + + /* Return x for convenience, but make sure result is never used. */ + return x; +} + +/* Special case function wrapper. */ +static float32x4_t VPCS_ATTR NOINLINE +special_case (float32x4_t x, float32x4_t y, float32x4_t ret, uint32x4_t cmp) +{ + return v_call2_f32 (powrf_specialcase, x, y, ret, cmp); +} + +/* Power implementation for x containing negative or subnormal lanes. */ +static inline float32x4_t +v_powrf_x_is_neg_or_sub (float32x4_t x, float32x4_t y, const struct data *d) +{ + uint32x4_t xsmall = vcaltq_f32 (x, v_f32 (0x1p-126f)); + + /* Normalize subnormals. */ + float32x4_t a = vabsq_f32 (x); + uint32x4_t ia_norm = vreinterpretq_u32_f32 (vmulq_f32 (a, d->norm)); + ia_norm = vsubq_u32 (ia_norm, d->subnormal_bias); + a = vbslq_f32 (xsmall, vreinterpretq_f32_u32 (ia_norm), a); + + /* Evaluate exp (y * log(x)) using |x| and sign bias correction. */ + float32x4_t ret = v_powrf_core (a, y, d); + + /* Cases of finite y and finite negative x. */ + uint32x4_t xisneg = vcltzq_f32 (x); + return vbslq_f32 (xisneg, d->nan, ret); +} + +/* Implementation of AdvSIMD powrf. + + powr(x,y) := exp(y * log (x)) + + This means powr(x,y) core computation matches that of pow(x,y) + but powr returns NaN for negative x even if y is an integer. + + Maximum measured error is 2.57 ULPs: + V_NAME_F2 (powr) (0x1.031706p+0, 0x1.ce2ec2p+12) + got 0x1.fff868p+127 + want 0x1.fff862p+127. */ +float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (powr) (float32x4_t x, float32x4_t y) +{ + const struct data *d = ptr_barrier (&data); + + /* Special cases of x or y: zero, inf and nan. */ + uint32x4_t ix = vreinterpretq_u32_f32 (x); + uint32x4_t iy = vreinterpretq_u32_f32 (y); + uint32x4_t xspecial = v_zeroinfnan (d, ix); + uint32x4_t yspecial = v_zeroinfnan (d, iy); + uint32x4_t cmp = vorrq_u32 (xspecial, yspecial); + + /* Evaluate pow(x, y) for x containing negative or subnormal lanes. */ + uint32x4_t x_is_neg_or_sub = vcltq_f32 (x, v_f32 (0x1p-126f)); + if (__glibc_unlikely (v_any_u32 (x_is_neg_or_sub))) + { + float32x4_t ret = v_powrf_x_is_neg_or_sub (x, y, d); + if (__glibc_unlikely (v_any_u32 (cmp))) + return special_case (x, y, ret, cmp); + return ret; + } + + /* Else evaluate pow(x, y) for normal and positive x only. */ + if (__glibc_unlikely (v_any_u32 (cmp))) + return special_case (x, y, v_powrf_core (x, y, d), cmp); + return v_powrf_core (x, y, d); +} + +libmvec_hidden_def (V_NAME_F2 (powr)) +HALF_WIDTH_ALIAS_F2 (powr) diff --git a/sysdeps/aarch64/fpu/powrf_sve.c b/sysdeps/aarch64/fpu/powrf_sve.c new file mode 100644 index 0000000000..32891be0a0 --- /dev/null +++ b/sysdeps/aarch64/fpu/powrf_sve.c @@ -0,0 +1,135 @@ +/* Single-precision vector (SVE) powr function + + Copyright (C) 2026 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <https://www.gnu.org/licenses/>. */ + +#include "flt-32/math_config.h" +#include "sv_math.h" + +#define WANT_SV_POWF_SIGN_BIAS 0 +#include "sv_powf_inline.h" + +/* A scalar subroutine used to fix main powrf special cases. */ +static inline float +powrf_specialcase (float x, float y) +{ + uint32_t ix = asuint (x); + uint32_t iy = asuint (y); + /* |y| is 0, Inf or NaN. */ + if (__glibc_unlikely (zeroinfnan (iy))) + { + /* |x| or |y| is NaN. */ + if (2 * ix > 2u * 0x7f800000 || 2 * iy > 2u * 0x7f800000) + return __builtin_nanf (""); + /* |y| = 0. */ + if (2 * iy == 0) + { + /* |x| = 0 or Inf. */ + if ((2 * ix == 0) || (2 * ix == 2u * 0x7f800000)) + return __builtin_nanf (""); + /* x is finite. */ + return 1.0f; + } + /* |y| = Inf and x = 1.0. */ + if (ix == 0x3f800000) + return __builtin_nanf (""); + /* |x| < 1 and y = Inf or |x| > 1 and y = -Inf. */ + if ((2 * ix < 2 * 0x3f800000) == !(iy & 0x80000000)) + return 0.0f; + /* |y| = Inf and previous conditions not met. */ + return y * y; + } + /* x is 0, Inf or NaN. Negative x are handled in the core. */ + if (__glibc_unlikely (zeroinfnan (ix))) + { + float x2 = x * x; + return iy & 0x80000000 ? 1 / x2 : x2; + } + /* Return x for convenience, but make sure result is never used. */ + return x; +} + +/* Scalar fallback for special case routines with custom signature. */ +static svfloat32_t NOINLINE +sv_call_powrf_sc (svfloat32_t x1, svfloat32_t x2, svfloat32_t y, svbool_t cmp) +{ + return sv_call2_f32 (powrf_specialcase, x1, x2, y, cmp); +} + +/* Implementation of SVE powrf. + + Provides the same accuracy as AdvSIMD powf and powrf, since it relies on the + same algorithm. + + Maximum measured error is 2.57 ULPs: + SV_NAME_F2 (powr) (0x1.031706p+0, 0x1.ce2ec2p+12) + got 0x1.fff868p+127 + want 0x1.fff862p+127. */ +svfloat32_t SV_NAME_F2 (powr) (svfloat32_t x, svfloat32_t y, const svbool_t pg) +{ + const struct data *d = ptr_barrier (&data); + + svuint32_t vix = svreinterpret_u32 (x); + svuint32_t viy = svreinterpret_u32 (y); + + svbool_t xpos = svcmpge (pg, x, sv_f32 (0.0f)); + + /* Special cases of x or y: zero, inf and nan. */ + svbool_t xspecial = sv_zeroinfnan (xpos, vix); + svbool_t yspecial = sv_zeroinfnan (xpos, viy); + svbool_t cmp = svorr_z (xpos, xspecial, yspecial); + + /* Cases of subnormal x: |x| < 0x1p-126. */ + svbool_t x_is_subnormal = svaclt (xpos, x, d->small_bound); + if (__glibc_unlikely (svptest_any (xpos, x_is_subnormal))) + { + /* Normalize subnormal x so exponent becomes negative. */ + vix = svreinterpret_u32 (svmul_m (x_is_subnormal, x, 0x1p23f)); + vix = svsub_m (x_is_subnormal, vix, d->subnormal_bias); + } + + /* Part of core computation carried in working precision. */ + svuint32_t tmp = svsub_x (xpos, vix, d->off); + svuint32_t i + = svand_x (xpos, svlsr_x (xpos, tmp, (23 - V_POWF_LOG2_TABLE_BITS)), + V_POWF_LOG2_N - 1); + svuint32_t top = svand_x (xpos, tmp, 0xff800000); + svuint32_t iz = svsub_x (xpos, vix, top); + svint32_t k + = svasr_x (xpos, svreinterpret_s32 (top), (23 - V_POWF_EXP2_TABLE_BITS)); + + /* Compute core in extended precision and return intermediate ylogx results + to handle cases of underflow and underflow in exp. */ + svfloat32_t ylogx; + /* Pass a dummy sign_bias so we can re-use powf core. + The core is simplified by setting WANT_SV_POWF_SIGN_BIAS = 0. */ + svfloat32_t ret = sv_powf_core (xpos, i, iz, k, y, sv_u32 (0), &ylogx, d); + + /* Handle exp special cases of underflow and overflow. */ + svbool_t no_uflow = svcmpgt (xpos, ylogx, d->uflow_bound); + svbool_t oflow = svcmpgt (xpos, ylogx, d->oflow_bound); + svfloat32_t ret_flow = svdup_n_f32_z (no_uflow, INFINITY); + ret = svsel (svorn_z (xpos, oflow, no_uflow), ret_flow, ret); + + /* Cases of negative x. */ + ret = svsel (xpos, ret, sv_f32 (__builtin_nanf (""))); + + if (__glibc_unlikely (svptest_any (cmp, cmp))) + return sv_call_powrf_sc (x, y, ret, cmp); + + return ret; +} diff --git a/sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c b/sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c index 19adf79fde..74980a7f6f 100644 --- a/sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c +++ b/sysdeps/aarch64/fpu/test-double-advsimd-wrappers.c @@ -54,6 +54,7 @@ VPCS_VECTOR_WRAPPER (log1p_advsimd, _ZGVnN2v_log1p) VPCS_VECTOR_WRAPPER (log2_advsimd, _ZGVnN2v_log2) VPCS_VECTOR_WRAPPER (log2p1_advsimd, _ZGVnN2v_log2p1) VPCS_VECTOR_WRAPPER_ff (pow_advsimd, _ZGVnN2vv_pow) +VPCS_VECTOR_WRAPPER_ff (powr_advsimd, _ZGVnN2vv_powr) VPCS_VECTOR_WRAPPER (rsqrt_advsimd, _ZGVnN2v_rsqrt) VPCS_VECTOR_WRAPPER (sin_advsimd, _ZGVnN2v_sin) VPCS_VECTOR_WRAPPER (sinh_advsimd, _ZGVnN2v_sinh) diff --git a/sysdeps/aarch64/fpu/test-double-sve-wrappers.c b/sysdeps/aarch64/fpu/test-double-sve-wrappers.c index 86e73756a2..e6e3d652c9 100644 --- a/sysdeps/aarch64/fpu/test-double-sve-wrappers.c +++ b/sysdeps/aarch64/fpu/test-double-sve-wrappers.c @@ -73,6 +73,7 @@ SVE_VECTOR_WRAPPER (log1p_sve, _ZGVsMxv_log1p) SVE_VECTOR_WRAPPER (log2_sve, _ZGVsMxv_log2) SVE_VECTOR_WRAPPER (log2p1_sve, _ZGVsMxv_log2p1) SVE_VECTOR_WRAPPER_ff (pow_sve, _ZGVsMxvv_pow) +SVE_VECTOR_WRAPPER_ff (powr_sve, _ZGVsMxvv_powr) SVE_VECTOR_WRAPPER (rsqrt_sve, _ZGVsMxv_rsqrt) SVE_VECTOR_WRAPPER (sin_sve, _ZGVsMxv_sin) SVE_VECTOR_WRAPPER (sinh_sve, _ZGVsMxv_sinh) diff --git a/sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c b/sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c index 3bd3f5c950..223e491007 100644 --- a/sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c +++ b/sysdeps/aarch64/fpu/test-float-advsimd-wrappers.c @@ -54,6 +54,7 @@ VPCS_VECTOR_WRAPPER (log1pf_advsimd, _ZGVnN4v_log1pf) VPCS_VECTOR_WRAPPER (log2f_advsimd, _ZGVnN4v_log2f) VPCS_VECTOR_WRAPPER (log2p1f_advsimd, _ZGVnN4v_log2p1f) VPCS_VECTOR_WRAPPER_ff (powf_advsimd, _ZGVnN4vv_powf) +VPCS_VECTOR_WRAPPER_ff (powrf_advsimd, _ZGVnN4vv_powrf) VPCS_VECTOR_WRAPPER (rsqrtf_advsimd, _ZGVnN4v_rsqrtf) VPCS_VECTOR_WRAPPER (sinf_advsimd, _ZGVnN4v_sinf) VPCS_VECTOR_WRAPPER (sinhf_advsimd, _ZGVnN4v_sinhf) diff --git a/sysdeps/aarch64/fpu/test-float-sve-wrappers.c b/sysdeps/aarch64/fpu/test-float-sve-wrappers.c index 0d9a7e5b93..2a01f93b5a 100644 --- a/sysdeps/aarch64/fpu/test-float-sve-wrappers.c +++ b/sysdeps/aarch64/fpu/test-float-sve-wrappers.c @@ -73,6 +73,7 @@ SVE_VECTOR_WRAPPER (log1pf_sve, _ZGVsMxv_log1pf) SVE_VECTOR_WRAPPER (log2f_sve, _ZGVsMxv_log2f) SVE_VECTOR_WRAPPER (log2p1f_sve, _ZGVsMxv_log2p1f) SVE_VECTOR_WRAPPER_ff (powf_sve, _ZGVsMxvv_powf) +SVE_VECTOR_WRAPPER_ff (powrf_sve, _ZGVsMxvv_powrf) SVE_VECTOR_WRAPPER (rsqrtf_sve, _ZGVsMxv_rsqrtf) SVE_VECTOR_WRAPPER (sinf_sve, _ZGVsMxv_sinf) SVE_VECTOR_WRAPPER (sinhf_sve, _ZGVsMxv_sinhf) diff --git a/sysdeps/aarch64/fpu/v_powrf_inline.h b/sysdeps/aarch64/fpu/v_powrf_inline.h index 0168c32338..ea87a2b85b 100644 --- a/sysdeps/aarch64/fpu/v_powrf_inline.h +++ b/sysdeps/aarch64/fpu/v_powrf_inline.h @@ -1,6 +1,6 @@ /* Helper for AdvSIMD single-precision powr - Copyright (C) 2025 Free Software Foundation, Inc. + Copyright (C) 2025-2026 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -17,6 +17,8 @@ License along with the GNU C Library; if not, see <https://www.gnu.org/licenses/>. */ +#include "powf_common.h" + #define Log2IdxMask (V_POWF_LOG2_N - 1) #define Exp2IdxMask (V_POWF_EXP2_N - 1) #define Scale ((double) V_POWF_EXP2_N) diff --git a/sysdeps/unix/sysv/linux/aarch64/libmvec.abilist b/sysdeps/unix/sysv/linux/aarch64/libmvec.abilist index 6d13d53613..638e34e500 100644 --- a/sysdeps/unix/sysv/linux/aarch64/libmvec.abilist +++ b/sysdeps/unix/sysv/linux/aarch64/libmvec.abilist @@ -193,3 +193,8 @@ GLIBC_2.43 _ZGVsMxv_log2p1 F GLIBC_2.43 _ZGVsMxv_log2p1f F GLIBC_2.43 _ZGVsMxv_rsqrt F GLIBC_2.43 _ZGVsMxv_rsqrtf F +GLIBC_2.44 _ZGVnN2vv_powr F +GLIBC_2.44 _ZGVnN2vv_powrf F +GLIBC_2.44 _ZGVnN4vv_powrf F +GLIBC_2.44 _ZGVsMxvv_powr F +GLIBC_2.44 _ZGVsMxvv_powrf F