[glibc] Vectorise special cases for SVE inverse hyperbolics

https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=f7eea1b2585b7b8c1f3ee9646627a9cf03fc33f2

commit f7eea1b2585b7b8c1f3ee9646627a9cf03fc33f2
Author: Thomas Daubney <[email protected]>
Date:   Fri May 8 15:58:11 2026 +0000

    Vectorise special cases for SVE inverse hyperbolics
    
    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.
    Reviewed-by: Adhemerval Zanella  <[email protected]>

Diff:
---
 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