[PATCH] AArch64: Optimise SVE scalar callbacks

[Joe Ramsay]

Instead of using SVE instructions to marshall special results into the
correct lane, just write the entire vector (and the predicate) to
memory, then use cheaper scalar operations.

-fno-stack-protector is added to the CFLAGS for SVE routines, as the
fallbacks incur a small stack allocation. sv_call_* allocates VL + VL
/ 8 bytes, and sv_call2_* allocates 2 * VL + VL / 8, where VL is the
SVE vector length in bytes. The maximum VL is 256, so the largest
possible stack allocation from the new handlers is 544 bytes.

Though scalar fallbacks are supposed to be unlikely, they are still
relevant at Ofast for certain routines where they are used for
'difficult' intervals. Below is the effect on throughput (ns/elem) on
a subset of libmvec, with GCC 14 on Neoverse V1, in the relevant
intervals.

Func   Low    High      After/Before
====================================
log    0      0x1p-1022 0.92x
logf   0      0x1p-126  0.81x
sin    0x1p23 inf       1.00x
sinf   0x1p20 inf       0.91x
cos    0x1p23 inf       0.95x
cosf   0x1p20 inf       0.66x
atan2  inf    nan       0.97x
atan2f inf    nan       0.80x
---
OK for master? If so please commit for me as I don't have commit rights.
Thanks,
Joe
 sysdeps/aarch64/fpu/Makefile  |   3 +
 sysdeps/aarch64/fpu/sv_math.h | 100 ++++++++++++++++++++++------------
 2 files changed, 69 insertions(+), 34 deletions(-)

diff --git a/sysdeps/aarch64/fpu/Makefile b/sysdeps/aarch64/fpu/Makefile
index 1ba04590cf..b2beee706e 100644
--- a/sysdeps/aarch64/fpu/Makefile
+++ b/sysdeps/aarch64/fpu/Makefile
@@ -64,6 +64,9 @@ endif
 
 sve-cflags = -mcpu=generic+sve -march=armv8-a+sve -mtune=neoverse-v1
 
+# Add -fno-stack-protector for optimal code in the SVE scalar fallbacks
+sve-cflags += -fno-stack-protector
+
 ifeq ($(build-mathvec),yes)
 bench-libmvec = $(addprefix float-advsimd-,$(float-advsimd-funcs)) \
                 $(addprefix double-advsimd-,$(double-advsimd-funcs)) \
diff --git a/sysdeps/aarch64/fpu/sv_math.h b/sysdeps/aarch64/fpu/sv_math.h
index 3d576df4cc..d5e7e1765b 100644
--- a/sysdeps/aarch64/fpu/sv_math.h
+++ b/sysdeps/aarch64/fpu/sv_math.h
@@ -29,6 +29,13 @@
 #define SV_NAME_F2(fun) _ZGVsMxvv_##fun##f
 #define SV_NAME_D2(fun) _ZGVsMxvv_##fun
 
+static inline void
+svstr_p (uint8_t *dst, svbool_t p)
+{
+  /* Predicate STR does not currently have an intrinsic.  */
+  __asm__("str %0, [%x1]\n" : : "Upa"(p), "r"(dst) : "memory");
+}
+
 /* Double precision.  */
 static inline svint64_t
 sv_s64 (int64_t x)
@@ -51,33 +58,39 @@ sv_f64 (double x)
 static inline svfloat64_t
 sv_call_f64 (double (*f) (double), svfloat64_t x, svfloat64_t y, svbool_t cmp)
 {
-  svbool_t p = svpfirst (cmp, svpfalse ());
-  while (svptest_any (cmp, p))
+  double tmp[svcntd ()];
+  uint8_t pg_bits[svcntd ()];
+  svstr_p (pg_bits, cmp);
+  svst1 (svptrue_b64 (), tmp, svsel (cmp, x, y));
+
+  for (int i = 0; i < svcntd (); i++)
     {
-      double elem = svclastb_n_f64 (p, 0, x);
-      elem = (*f) (elem);
-      svfloat64_t y2 = svdup_n_f64 (elem);
-      y = svsel_f64 (p, y2, y);
-      p = svpnext_b64 (cmp, p);
+      if (pg_bits[i])
+	{
+	  tmp[i] = f (tmp[i]);
+	}
     }
-  return y;
+  return svld1 (svptrue_b64 (), tmp);
 }
 
 static inline svfloat64_t
 sv_call2_f64 (double (*f) (double, double), svfloat64_t x1, svfloat64_t x2,
 	      svfloat64_t y, svbool_t cmp)
 {
-  svbool_t p = svpfirst (cmp, svpfalse ());
-  while (svptest_any (cmp, p))
+  double tmp1[svcntd ()], tmp2[svcntd ()];
+  uint8_t pg_bits[svcntd ()];
+  svstr_p (pg_bits, cmp);
+  svst1 (svptrue_b64 (), tmp1, svsel (cmp, x1, y));
+  svst1 (cmp, tmp2, x2);
+
+  for (int i = 0; i < svcntd (); i++)
     {
-      double elem1 = svclastb_n_f64 (p, 0, x1);
-      double elem2 = svclastb_n_f64 (p, 0, x2);
-      double ret = (*f) (elem1, elem2);
-      svfloat64_t y2 = svdup_n_f64 (ret);
-      y = svsel_f64 (p, y2, y);
-      p = svpnext_b64 (cmp, p);
+      if (pg_bits[i])
+	{
+	  tmp1[i] = f (tmp1[i], tmp2[i]);
+	}
     }
-  return y;
+  return svld1 (svptrue_b64 (), tmp1);
 }
 
 static inline svuint64_t
@@ -109,33 +122,52 @@ sv_f32 (float x)
 static inline svfloat32_t
 sv_call_f32 (float (*f) (float), svfloat32_t x, svfloat32_t y, svbool_t cmp)
 {
-  svbool_t p = svpfirst (cmp, svpfalse ());
-  while (svptest_any (cmp, p))
+  float tmp[svcntw ()];
+  /* svcntd, not svcntw, is correct for pg_bits because each bit of pg_bits
+     maps to 1 byte of the vector, so a uint8_t indicates predication of two
+     floats.  */
+  uint8_t pg_bits[svcntd ()];
+  svstr_p (pg_bits, cmp);
+  svst1 (svptrue_b32 (), tmp, svsel (cmp, x, y));
+
+  for (int i = 0; i < svcntd (); i++)
     {
-      float elem = svclastb_n_f32 (p, 0, x);
-      elem = f (elem);
-      svfloat32_t y2 = svdup_n_f32 (elem);
-      y = svsel_f32 (p, y2, y);
-      p = svpnext_b32 (cmp, p);
+      uint8_t p = pg_bits[i];
+      if (p & 1)
+	{
+	  tmp[i * 2] = f (tmp[i * 2]);
+	}
+      if (p & (1 << 4))
+	{
+	  tmp[i * 2 + 1] = f (tmp[i * 2 + 1]);
+	}
     }
-  return y;
+  return svld1 (svptrue_b32 (), tmp);
 }
 
 static inline svfloat32_t
 sv_call2_f32 (float (*f) (float, float), svfloat32_t x1, svfloat32_t x2,
 	      svfloat32_t y, svbool_t cmp)
 {
-  svbool_t p = svpfirst (cmp, svpfalse ());
-  while (svptest_any (cmp, p))
+  float tmp1[svcntw ()], tmp2[svcntw ()];
+  uint8_t pg_bits[svcntd ()];
+  svstr_p (pg_bits, cmp);
+  svst1 (svptrue_b32 (), tmp1, svsel (cmp, x1, y));
+  svst1 (cmp, tmp2, x2);
+
+  for (int i = 0; i < svcntd (); i++)
     {
-      float elem1 = svclastb_n_f32 (p, 0, x1);
-      float elem2 = svclastb_n_f32 (p, 0, x2);
-      float ret = f (elem1, elem2);
-      svfloat32_t y2 = svdup_n_f32 (ret);
-      y = svsel_f32 (p, y2, y);
-      p = svpnext_b32 (cmp, p);
+      uint8_t p = pg_bits[i];
+      if (p & 1)
+	{
+	  tmp1[i * 2] = f (tmp1[i * 2], tmp2[i * 2]);
+	}
+      if (p & (1 << 4))
+	{
+	  tmp1[i * 2 + 1] = f (tmp1[i * 2 + 1], tmp2[i * 2 + 1]);
+	}
     }
-  return y;
+  return svld1 (svptrue_b32 (), tmp1);
 }
 
 #endif
-- 
2.34.1