[PATCH] Vectorise special cases for SVE log1p(f)

[Thomas Daubney]

This patch adds vectorised special cases for the SVE log functions
log1p and log1pf.

When built with GCC-15 and executed on a Neoverse V2 platform, the
following benchmarking throughput uplifts were measured:

log1pf -> 285% speed-up (4.85 ns/element to 1.26 ns/element)
log1p  -> 117% speed-up (8.25 ns/element to 3.80 ns/element)

Note that the numbers here are for the special case path only and that
the fast path performance has been maintained. These changes have also
maintained the same level of accuracy as before.
---
 sysdeps/aarch64/fpu/log1p_sve.c        | 12 ++++++++----
 sysdeps/aarch64/fpu/log1pf_sve.c       | 23 ++++++++++++++++-------
 sysdeps/aarch64/fpu/sv_log1pf_inline.h |  5 ++---
 3 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/sysdeps/aarch64/fpu/log1p_sve.c b/sysdeps/aarch64/fpu/log1p_sve.c
index bde0e8a7d1..331cfb34d3 100644
--- a/sysdeps/aarch64/fpu/log1p_sve.c
+++ b/sysdeps/aarch64/fpu/log1p_sve.c
@@ -63,9 +63,13 @@ static const struct data
 #define BottomMask 0xffffffff
 
 static svfloat64_t NOINLINE
-special_case (svfloat64_t x, svfloat64_t y, svbool_t special)
+special_case (svfloat64_t x, svfloat64_t y, svbool_t special, svbool_t pg)
 {
-  return sv_call_f64 (log1p, x, y, special);
+  y = svsel (special, sv_f64 (NAN), y);
+  svbool_t ret_pinf = svcmpeq (pg, x, INFINITY);
+  svbool_t ret_minf = svcmpeq (pg, x, -1.0);
+  y = svsel (ret_pinf, sv_f64 (INFINITY), y);
+  return svsel (ret_minf, sv_f64 (-INFINITY), y);
 }
 
 /* Vector approximation for log1p using polynomial on reduced interval. Maximum
@@ -159,8 +163,8 @@ svfloat64_t SV_NAME_D1 (log1p) (svfloat64_t x, svbool_t pg)
   if (__glibc_unlikely (svptest_any (pg, special)))
     return special_case (
 	x, svmla_x (svptrue_b64 (), svadd_x (svptrue_b64 (), ylo, yhi), f2, p),
-	special);
+	special, pg);
   return svmla_x (svptrue_b64 (), svadd_x (svptrue_b64 (), ylo, yhi), f2, p);
 }
 
-strong_alias (SV_NAME_D1 (log1p), SV_NAME_D1 (logp1))
+strong_alias (SV_NAME_D1 (log1p), SV_NAME_D1 (logp1))
\ No newline at end of file
diff --git a/sysdeps/aarch64/fpu/log1pf_sve.c b/sysdeps/aarch64/fpu/log1pf_sve.c
index f4d409eb50..297f5374ca 100644
--- a/sysdeps/aarch64/fpu/log1pf_sve.c
+++ b/sysdeps/aarch64/fpu/log1pf_sve.c
@@ -17,14 +17,21 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
+#include "../../ieee754/flt-32/math_config.h"
 #include "sv_math.h"
 #include "sv_log1pf_inline.h"
 
 static svfloat32_t NOINLINE
-special_case (svfloat32_t x, svbool_t special)
+special_case (svfloat32_t x, svfloat32_t y, svbool_t pg, svbool_t special,
+	      const struct sv_log1pf_data *d)
 {
-  return sv_call_f32 (log1pf, x, sv_log1pf_inline (x, svptrue_b32 ()),
-		      special);
+  y = svsel_f32 (special, svreinterpret_f32 (sv_u32 (d->nan)), y);
+
+  svbool_t ret_pinf = svcmpeq (pg, x, asfloat (d->inf));
+  svbool_t ret_minf = svcmpeq (pg, x, -1.0f);
+
+  y = svsel_f32 (ret_pinf, svreinterpret_f32 (sv_u32 (d->inf)), y);
+  return svsel_f32 (ret_minf, sv_f32 (-d->inf), y);
 }
 
 /* Vector log1pf approximation using polynomial on reduced interval. Worst-case
@@ -33,14 +40,16 @@ special_case (svfloat32_t x, svbool_t special)
 				 want 0x1.9f323ep-2.  */
 svfloat32_t SV_NAME_F1 (log1p) (svfloat32_t x, svbool_t pg)
 {
+
+  const struct sv_log1pf_data *d = ptr_barrier (&sv_log1pf_data);
   /* x < -1, Inf/Nan.  */
-  svbool_t special = svcmpeq (pg, svreinterpret_u32 (x), 0x7f800000);
-  special = svorn_z (pg, special, svcmpge (pg, x, -1));
+  svbool_t special = svcmpeq (pg, svreinterpret_u32 (x), d->inf);
+  special = svorn_z (pg, special, svcmpge (pg, x, -1.0f));
 
   if (__glibc_unlikely (svptest_any (pg, special)))
-    return special_case (x, special);
+    return special_case (x, sv_log1pf_inline (x, pg), pg, special, d);
 
   return sv_log1pf_inline (x, pg);
 }
 
-strong_alias (SV_NAME_F1 (log1p), SV_NAME_F1 (logp1))
+strong_alias (SV_NAME_F1 (log1p), SV_NAME_F1 (logp1))
\ No newline at end of file
diff --git a/sysdeps/aarch64/fpu/sv_log1pf_inline.h b/sysdeps/aarch64/fpu/sv_log1pf_inline.h
index 494fa279f8..3677578c59 100644
--- a/sysdeps/aarch64/fpu/sv_log1pf_inline.h
+++ b/sysdeps/aarch64/fpu/sv_log1pf_inline.h
@@ -21,7 +21,6 @@
 #define AARCH64_FPU_SV_LOG1PF_INLINE_H
 
 #include "sv_math.h"
-#include "vecmath_config.h"
 
 #define SignExponentMask 0xff800000
 
@@ -31,6 +30,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 +39,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 +93,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
-- 
2.43.0