[PATCH] Vectorise special cases for SVE inverse hyperbolics
Adhemerval Zanella Netto
adhemerval.zanella@linaro.org
Thu May 21 12:06:39 GMT 2026
On 08/05/26 12:58, Thomas Daubney wrote:
> This patch adds vectorised special cases for the SVE inverse hyperbolic
> functions atanh, acosh and asinh for single precision floats. It also
> moves the commonly used inf and nan bit values into the sv_log1pf_inline
> data struct for resuse.
>
> When built with GCC-15 and executed on a Neoverse V2 platform, the
> following benchmarking throughput uplifts were measured:
>
> atanh -> 215% speed-up (5.51 ns/element to 1.75 ns/element)
> acosh -> 152% speed-up (4.63 ns/element to 1.84 ns/element)
> asinh -> 51% speed-up (5.00 ns/element to 3.31 ns/element)
>
> Note that the numbers here are for the special case path only and that
> the fast path performance has been maintained.
LGTM, thanks.
Reviewed-by: Adhemerval Zanella <adhemerval.zanella@linaro.org>
> ---
> sysdeps/aarch64/fpu/acoshf_sve.c | 43 +++++++++++++++++++-------
> sysdeps/aarch64/fpu/asinhf_sve.c | 36 ++++++++++++++++-----
> sysdeps/aarch64/fpu/atanhf_sve.c | 35 +++++++++++----------
> sysdeps/aarch64/fpu/sv_log1pf_inline.h | 4 +--
> 4 files changed, 80 insertions(+), 38 deletions(-)
>
> diff --git a/sysdeps/aarch64/fpu/acoshf_sve.c b/sysdeps/aarch64/fpu/acoshf_sve.c
> index e133db5490..e3a42e8fe4 100644
> --- a/sysdeps/aarch64/fpu/acoshf_sve.c
> +++ b/sysdeps/aarch64/fpu/acoshf_sve.c
> @@ -20,17 +20,36 @@
> #include "sv_math.h"
> #include "sv_log1pf_inline.h"
>
> -#define One 0x3f800000
> -#define Thres 0x20000000 /* asuint(0x1p64) - One. */
> +#define One (0x3f800000U)
> +#define Thres (0x20000000U) /* asuint(0x1p64) - One. */
>
> -static svfloat32_t NOINLINE
> -special_case (svfloat32_t xm1, svfloat32_t tmp, svbool_t special)
> +/* Acosh is defined on [1, inf). Its formula can be re-written knowing that 1
> + becomes negligible when x is a very large number. So for special numbers,
> + where x >= 2^64, acosh ~= ln(2x). But, ln(2x) = ln(2) + ln(x) and below we
> + calculate ln(x) and then add ln(2) to the result.
> +
> + Right before returning we check if x is infinity or if x is lower than 1,
> + in which case we return infinity or NaN. */
> +static inline svfloat32_t special_case (svfloat32_t x, svfloat32_t xm1,
> + svfloat32_t y, svbool_t pg,
> + svbool_t special,
> + const struct sv_log1pf_data *d)
> {
> - svfloat32_t x = svadd_x (svptrue_b32 (), xm1, 1.0f);
> - svfloat32_t y = sv_log1pf_inline (tmp, svptrue_b32 ());
> - return sv_call_f32 (acoshf, x, y, special);
> -}
> + svfloat32_t logv = sv_log1pf_inline (svsel (special, xm1, y), pg);
> + svfloat32_t result = svadd_m (special, logv, sv_f32 (d->ln2));
> +
> + /* Catch x<1, and x==inf.
> + Also catch x==nan using negation of fp comparison. */
> + svbool_t is_x_ge1 = svcmpge (special, xm1, 0.0f);
> + svbool_t is_x_pinf
> + = svcmpeq (special, x, svreinterpret_f32 (sv_u32 (d->inf)));
> +
> + svbool_t res_is_inf_nan = svorn_b_z (special, is_x_pinf, is_x_ge1);
>
> + svuint32_t inf_or_nan = svsel (is_x_pinf, sv_u32 (d->inf), sv_u32 (d->nan));
> +
> + return svsel (res_is_inf_nan, svreinterpret_f32 (inf_or_nan), result);
> +}
> /* Single-precision SVE acosh(x) routine. Implements the same algorithm as
> vector acoshf and log1p.
>
> @@ -39,14 +58,16 @@ special_case (svfloat32_t xm1, svfloat32_t tmp, svbool_t special)
> want 0x1.e435a2p-4. */
> svfloat32_t SV_NAME_F1 (acosh) (svfloat32_t x, const svbool_t pg)
> {
> + const struct sv_log1pf_data *d = ptr_barrier (&sv_log1pf_data);
> +
> svuint32_t ix = svreinterpret_u32 (x);
> svbool_t special = svcmpge (pg, svsub_x (pg, ix, One), Thres);
>
> svfloat32_t xm1 = svsub_x (pg, x, 1.0f);
> svfloat32_t u = svmul_x (pg, xm1, svadd_x (pg, x, 1.0f));
> - svfloat32_t tmp = svadd_x (pg, xm1, svsqrt_x (pg, u));
> + svfloat32_t y = svadd_x (pg, xm1, svsqrt_x (pg, u));
>
> if (__glibc_unlikely (svptest_any (pg, special)))
> - return special_case (xm1, tmp, special);
> - return sv_log1pf_inline (tmp, pg);
> + return special_case (x, xm1, y, pg, special, d);
> + return sv_log1pf_inline (y, pg);
> }
> diff --git a/sysdeps/aarch64/fpu/asinhf_sve.c b/sysdeps/aarch64/fpu/asinhf_sve.c
> index b41ed93d71..9f87c2d737 100644
> --- a/sysdeps/aarch64/fpu/asinhf_sve.c
> +++ b/sysdeps/aarch64/fpu/asinhf_sve.c
> @@ -23,12 +23,30 @@
> #define BigBound 0x5f800000 /* asuint(0x1p64). */
>
> static svfloat32_t NOINLINE
> -special_case (svuint32_t iax, svuint32_t sign, svfloat32_t y, svbool_t special)
> +special_case (svfloat32_t ax, svfloat32_t y, svuint32_t sign, svbool_t pg,
> + svbool_t special, const struct sv_log1pf_data *d)
> {
> - svfloat32_t x = svreinterpret_f32 (sveor_x (svptrue_b32 (), iax, sign));
> - y = svreinterpret_f32 (
> - svorr_x (svptrue_b32 (), sign, svreinterpret_u32 (y)));
> - return sv_call_f32 (asinhf, x, y, special);
> + /* For very large inputs (x > 2^64), asinh(x) ≈ ln(2x).
> + In this range the +sqrt(x^2+1) term is negligible, so we compute
> + asinh(x) as ln(x) + ln(2) later in this function. */
> + svfloat32_t log_ax = sv_log1pf_inline (ax, special);
> +
> + /* The only special cases that need considering are infinity and NaNs since
> + 0 will be handled by other calculations. */
> + svfloat32_t inf = svreinterpret_f32 (sv_u32 (d->inf));
> + svbool_t is_inf = svcmpeq (special, ax, inf);
> + svbool_t is_nan = svcmpne (special, ax, ax);
> + svfloat32_t inf_ln2 = svsel (is_inf, inf, sv_f32 (d->ln2));
> + svfloat32_t inf_nan_ln2
> + = svsel (is_nan, svreinterpret_f32 (sv_u32 (d->nan)), inf_ln2);
> + svfloat32_t asinh_x_res = svadd_x (special, log_ax, inf_nan_ln2);
> +
> + /* Now select (based on special) between x and y to change the type and,
> + return either the positive or negative value, considering the input and
> + its sign. */
> + svfloat32_t result = svsel (special, asinh_x_res, y);
> + svuint32_t result_uint = svreinterpret_u32 (result);
> + return svreinterpret_f32 (sveor_m (pg, result_uint, sign));
> }
>
> /* Single-precision SVE asinh(x) routine. Implements the same algorithm as
> @@ -39,6 +57,8 @@ special_case (svuint32_t iax, svuint32_t sign, svfloat32_t y, svbool_t special)
> want -0x1.fd0bc8p-2. */
> svfloat32_t SV_NAME_F1 (asinh) (svfloat32_t x, const svbool_t pg)
> {
> + const struct sv_log1pf_data *d = ptr_barrier (&sv_log1pf_data);
> +
> svfloat32_t ax = svabs_x (pg, x);
> svuint32_t iax = svreinterpret_u32 (ax);
> svuint32_t sign = sveor_x (pg, svreinterpret_u32 (x), iax);
> @@ -47,11 +67,11 @@ svfloat32_t SV_NAME_F1 (asinh) (svfloat32_t x, const svbool_t pg)
> /* asinh(x) = log(x + sqrt(x * x + 1)).
> For positive x, asinh(x) = log1p(x + x * x / (1 + sqrt(x * x + 1))). */
> svfloat32_t ax2 = svmul_x (pg, ax, ax);
> - svfloat32_t d = svadd_x (pg, svsqrt_x (pg, svadd_x (pg, ax2, 1.0f)), 1.0f);
> + svfloat32_t dx = svadd_x (pg, svsqrt_x (pg, svadd_x (pg, ax2, 1.0f)), 1.0f);
> svfloat32_t y
> - = sv_log1pf_inline (svadd_x (pg, ax, svdiv_x (pg, ax2, d)), pg);
> + = sv_log1pf_inline (svadd_x (pg, ax, svdiv_x (pg, ax2, dx)), pg);
>
> if (__glibc_unlikely (svptest_any (pg, special)))
> - return special_case (iax, sign, y, special);
> + return special_case (ax, y, sign, pg, special, d);
> return svreinterpret_f32 (svorr_x (pg, sign, svreinterpret_u32 (y)));
> }
> diff --git a/sysdeps/aarch64/fpu/atanhf_sve.c b/sysdeps/aarch64/fpu/atanhf_sve.c
> index 8e07fe5a82..6c245c6d36 100644
> --- a/sysdeps/aarch64/fpu/atanhf_sve.c
> +++ b/sysdeps/aarch64/fpu/atanhf_sve.c
> @@ -20,16 +20,14 @@
> #include "sv_math.h"
> #include "sv_log1pf_inline.h"
>
> -#define One (0x3f800000)
> -#define Half (0x3f000000)
> -
> -static svfloat32_t NOINLINE
> -special_case (svuint32_t iax, svuint32_t sign, svfloat32_t halfsign,
> - svfloat32_t y, svbool_t special)
> +static inline svfloat32_t special_case (svfloat32_t ax, svfloat32_t y,
> + svbool_t special, svfloat32_t halfsign,
> + const struct sv_log1pf_data *d)
> {
> - svfloat32_t x = svreinterpret_f32 (sveor_x (svptrue_b32 (), iax, sign));
> - y = svmul_x (svptrue_b32 (), halfsign, y);
> - return sv_call_f32 (atanhf, x, y, special);
> + svfloat32_t res = svsel (special, svreinterpret_f32 (sv_u32 (d->nan)), y);
> + res = svsel (svcmpeq (special, ax, sv_f32 (1.0)),
> + svreinterpret_f32 (sv_u32 (d->inf)), res);
> + return svmul_x (svptrue_b32 (), res, halfsign);
> }
>
> /* Approximation for vector single-precision atanh(x) using modified log1p.
> @@ -38,21 +36,24 @@ special_case (svuint32_t iax, svuint32_t sign, svfloat32_t halfsign,
> want 0x1.f1f4f6p-5. */
> svfloat32_t SV_NAME_F1 (atanh) (svfloat32_t x, const svbool_t pg)
> {
> + const struct sv_log1pf_data *d = ptr_barrier (&sv_log1pf_data);
> +
> svfloat32_t ax = svabs_x (pg, x);
> svuint32_t iax = svreinterpret_u32 (ax);
> - svuint32_t sign = sveor_x (pg, svreinterpret_u32 (x), iax);
> - svfloat32_t halfsign = svreinterpret_f32 (svorr_x (pg, sign, Half));
> - svbool_t special = svcmpge (pg, iax, One);
> + svuint32_t sign = sveor_x (svptrue_b32 (), svreinterpret_u32 (x), iax);
> + svfloat32_t halfsign
> + = svreinterpret_f32 (svorr_x (svptrue_b32 (), sign, 0x3f000000));
> + svbool_t special = svcmpge (pg, iax, svreinterpret_u32 (sv_f32 (1)));
>
> /* Computation is performed based on the following sequence of equality:
> - * (1+x)/(1-x) = 1 + 2x/(1-x). */
> - svfloat32_t y = svadd_x (pg, ax, ax);
> - y = svdiv_x (pg, y, svsub_x (pg, sv_f32 (1), ax));
> + (1+x)/(1-x) = 1 + 2x/(1-x). */
> + svfloat32_t y = svadd_x (svptrue_b32 (), ax, ax);
> + y = svdiv_x (pg, y, svsubr_x (pg, ax, 1.0f));
> /* ln((1+x)/(1-x)) = ln(1+2x/(1-x)) = ln(1 + y). */
> y = sv_log1pf_inline (y, pg);
>
> if (__glibc_unlikely (svptest_any (pg, special)))
> - return special_case (iax, sign, halfsign, y, special);
> + return special_case (ax, y, special, halfsign, d);
>
> - return svmul_x (pg, halfsign, y);
> + return svmul_x (svptrue_b32 (), halfsign, y);
> }
> diff --git a/sysdeps/aarch64/fpu/sv_log1pf_inline.h b/sysdeps/aarch64/fpu/sv_log1pf_inline.h
> index 494fa279f8..0043df952c 100644
> --- a/sysdeps/aarch64/fpu/sv_log1pf_inline.h
> +++ b/sysdeps/aarch64/fpu/sv_log1pf_inline.h
> @@ -31,6 +31,7 @@ static const struct sv_log1pf_data
> float c1, c3, c5, c7;
> float ln2, exp_bias, quarter;
> uint32_t four, three_quarters;
> + uint32_t inf, nan;
> } sv_log1pf_data = {
> /* Do not store first term of polynomial, which is -0.5, as
> this can be fmov-ed directly instead of including it in
> @@ -39,7 +40,7 @@ static const struct sv_log1pf_data
> .c3 = -0x1.54ef78p-3f, .c4 = 0x1.28a1f4p-3f, .c5 = -0x1.0da91p-3f,
> .c6 = 0x1.abcb6p-4f, .c7 = -0x1.6f0d5ep-5f, .ln2 = 0x1.62e43p-1f,
> .exp_bias = 0x1p-23f, .quarter = 0x1p-2f, .four = 0x40800000,
> - .three_quarters = 0x3f400000,
> + .three_quarters = 0x3f400000, .inf = 0x7f800000, .nan = 0x7fc00000,
> };
>
> static inline svfloat32_t
> @@ -93,5 +94,4 @@ sv_log1pf_inline (svfloat32_t x, svbool_t pg)
> svfloat32_t scale_back = svmul_lane_f32 (svcvt_f32_x (pg, k), fconst, 1);
> return svmla_lane_f32 (p, scale_back, fconst, 0);
> }
> -
> #endif
More information about the Libc-alpha
mailing list