[glibc] AArch64: Vectorise SVE fp32 sin/cos special case

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

commit e585d88982e79fbc117c1b1e7bec541c70981182
Author: Dylan Fleming <[email protected]>
Date:   Thu Jun 18 14:55:15 2026 +0000

    AArch64: Vectorise SVE fp32 sin/cos special case
    
    Implemented single precision vector fallback for SVE sinf/cosf.
    
    Includes a rewrite of SVE sinf to use the trig instructions.
    
    Reviewed-by: Wilco Dijkstra  <[email protected]>

Diff:
---
 sysdeps/aarch64/fpu/cosf_sve.c          |  75 ++++++++-----
 sysdeps/aarch64/fpu/sinf_sve.c          | 137 +++++++++++++-----------
 sysdeps/aarch64/fpu/sv_trigf_fallback.h | 184 ++++++++++++++++++++++++++++++++
 3 files changed, 309 insertions(+), 87 deletions(-)

diff --git a/sysdeps/aarch64/fpu/cosf_sve.c b/sysdeps/aarch64/fpu/cosf_sve.c
index 680d89a8fa..9aedcb14fa 100644
--- a/sysdeps/aarch64/fpu/cosf_sve.c
+++ b/sysdeps/aarch64/fpu/cosf_sve.c
@@ -18,10 +18,11 @@
    <https://www.gnu.org/licenses/>.  */
 
 #include "sv_math.h"
+#include "sv_trigf_fallback.h"
 
 static const struct data
 {
-  float neg_pio2_1, neg_pio2_2, neg_pio2_3, inv_pio2, shift;
+  float neg_pio2_1, neg_pio2_2, neg_pio2_3, inv_pio2, shift, range_val;
 } data = {
   /* Polynomial coefficients are hard-wired in FTMAD instructions.  */
   .neg_pio2_1 = -0x1.921fb6p+0f,
@@ -31,46 +32,71 @@ static const struct data
   /* Original shift used in AdvSIMD cosf,
      plus a contribution to set the bit #0 of q
      as expected by trigonometric instructions.  */
-  .shift = 0x1.800002p+23f
+  .shift = 0x1.800002p+23f,
+  .range_val = 0x1p20f,
 };
 
-#define RangeVal 0x49800000 /* asuint32(0x1p20f).  */
-
 static svfloat32_t NOINLINE
-special_case (svfloat32_t x, svfloat32_t y, svbool_t oob)
+special_case (svfloat32_t x, svfloat32_t y, svbool_t special)
 {
-  return sv_call_f32 (cosf, x, y, oob);
+  special = svaclt (special, x, sv_f32 (INFINITY));
+
+  svfloat32x2_t reduction = large_range_reduction (svptrue_b32 (), x);
+
+  /* Unpack the quadrant from the return struct.  */
+  svuint32_t quadrant = svreinterpret_u32 (svget2 (reduction, 1));
+  svfloat32_t r = svget2 (reduction, 0);
+
+  /* Adjust quadrant to select cosine polynomial.  */
+  quadrant = svadd_x (svptrue_b32 (), quadrant, 1);
+
+  svfloat32_t f = svtssel (r, quadrant);
+  svfloat32_t r2 = svtsmul (r, quadrant);
+  svfloat32_t cos = sv_f32 (0.0f);
+  cos = svtmad (cos, r2, 4);
+  cos = svtmad (cos, r2, 3);
+  cos = svtmad (cos, r2, 2);
+  cos = svtmad (cos, r2, 1);
+  cos = svtmad (cos, r2, 0);
+  cos = svmul_x (svptrue_b32 (), f, cos);
+
+  return svsel (special, cos, y);
 }
 
-/* A fast SVE implementation of cosf based on trigonometric
-   instructions (FTMAD, FTSSEL, FTSMUL).
-   Maximum measured error: 2.06 ULPs.
-   SV_NAME_F1 (cos)(0x1.dea2f2p+19) got 0x1.fffe7ap-6
-				   want 0x1.fffe76p-6.  */
+/* Vector version of cosf.
+   The maximum observed error is 1.56 + 0.5 ULP if |x| < 0x1p20.
+   _ZGVsMxv_cosf (0x1.dea2f2p+19)
+    got 0x1.fffe7ap-6
+   want 0x1.fffe76p-6
+   The special domain has a higher maximum error than the fast path:
+   Maximum observed error is 2.65 + 0.5ULP
+   _ZGVsMxv_cosf (0x1.ff3afcp+53)
+    got -0x1.ffe74p-3
+   want -0x1.ffe73ap-3.  */
 svfloat32_t SV_NAME_F1 (cos) (svfloat32_t x, const svbool_t pg)
 {
   const struct data *d = ptr_barrier (&data);
-
-  svfloat32_t r = svabs_x (pg, x);
-  svbool_t oob = svcmpge (pg, svreinterpret_u32 (r), RangeVal);
+  svbool_t ptrue = svptrue_b32 ();
 
   /* Load some constants in quad-word chunks to minimise memory access.  */
-  svfloat32_t negpio2_and_invpio2 = svld1rq (svptrue_b32 (), &d->neg_pio2_1);
+  svfloat32_t negpio2_and_invpio2 = svld1rq (ptrue, &d->neg_pio2_1);
 
-  /* n = rint(|x|/(pi/2)).  */
-  svfloat32_t q = svmla_lane (sv_f32 (d->shift), r, negpio2_and_invpio2, 3);
-  svfloat32_t n = svsub_x (pg, q, d->shift);
+  /* n = rint(x/(pi/2)).  */
+  svfloat32_t q = svmla_lane (sv_f32 (d->shift), x, negpio2_and_invpio2, 3);
+  svfloat32_t n = svsub_x (ptrue, q, d->shift);
 
-  /* r = |x| - n*(pi/2)  (range reduction into -pi/4 .. pi/4).  */
+  /* r = x - n*(pi/2)  (range reduction into -pi/4 .. pi/4).  */
+  svfloat32_t r = x;
   r = svmla_lane (r, n, negpio2_and_invpio2, 0);
   r = svmla_lane (r, n, negpio2_and_invpio2, 1);
   r = svmla_lane (r, n, negpio2_and_invpio2, 2);
 
   /* Final multiplicative factor: 1.0 or x depending on bit #0 of q.  */
-  svfloat32_t f = svtssel (r, svreinterpret_u32 (q));
+  svuint32_t q_u = svreinterpret_u32 (q);
+  svfloat32_t f = svtssel (r, q_u);
 
   /* cos(r) poly approx.  */
-  svfloat32_t r2 = svtsmul (r, svreinterpret_u32 (q));
+  svfloat32_t r2 = svtsmul (r, q_u);
   svfloat32_t y = sv_f32 (0.0f);
   y = svtmad (y, r2, 4);
   y = svtmad (y, r2, 3);
@@ -78,8 +104,9 @@ svfloat32_t SV_NAME_F1 (cos) (svfloat32_t x, const svbool_t pg)
   y = svtmad (y, r2, 1);
   y = svtmad (y, r2, 0);
 
-  if (__glibc_unlikely (svptest_any (pg, oob)))
-    return special_case (x, svmul_x (svnot_z (pg, oob), f, y), oob);
+  svbool_t cmp = svacge (pg, x, sv_f32 (d->range_val));
+  if (__glibc_unlikely (svptest_any (pg, cmp)))
+    return special_case (x, svmul_x (ptrue, f, y), cmp);
   /* Apply factor.  */
-  return svmul_x (pg, f, y);
+  return svmul_x (ptrue, f, y);
 }
diff --git a/sysdeps/aarch64/fpu/sinf_sve.c b/sysdeps/aarch64/fpu/sinf_sve.c
index 5f5aa60924..d4668bfafd 100644
--- a/sysdeps/aarch64/fpu/sinf_sve.c
+++ b/sysdeps/aarch64/fpu/sinf_sve.c
@@ -18,81 +18,92 @@
    <https://www.gnu.org/licenses/>.  */
 
 #include "sv_math.h"
+#include "sv_trigf_fallback.h"
 
 static const struct data
 {
-  float poly[4];
-  /* Pi-related values to be loaded as one quad-word and used with
-     svmla_lane.  */
-  float negpi1, negpi2, negpi3, invpi;
-  float shift;
+  float neg_pio2_1, neg_pio2_2, neg_pio2_3, inv_pio2, shift, range_val;
 } data = {
-  .poly = {
-    /* Non-zero coefficients from the degree 9 Taylor series expansion of
-       sin.  */
-    -0x1.555548p-3f, 0x1.110df4p-7f, -0x1.9f42eap-13f, 0x1.5b2e76p-19f
-  },
-  .negpi1 = -0x1.921fb6p+1f,
-  .negpi2 = 0x1.777a5cp-24f,
-  .negpi3 = 0x1.ee59dap-49f,
-  .invpi = 0x1.45f306p-2f,
-  .shift = 0x1.8p+23f
+  /* Polynomial coefficients are hard-wired in FTMAD instructions.  */
+  .neg_pio2_1 = -0x1.921fb6p+0f,
+  .neg_pio2_2 = 0x1.777a5cp-25f,
+  .neg_pio2_3 = 0x1.ee59dap-50f,
+  .inv_pio2 = 0x1.45f306p-1f,
+  /* Original shift used in AdvSIMD cosf,
+     plus a contribution to set the bit #0 of q
+     as expected by trigonometric instructions.  */
+  .shift = 0x1.8p+23f,
+  .range_val = 0x1p20f,
 };
 
-#define RangeVal 0x49800000 /* asuint32 (0x1p20f).  */
-#define C(i) sv_f32 (d->poly[i])
-
 static svfloat32_t NOINLINE
-special_case (svfloat32_t x, svfloat32_t y, svbool_t cmp)
+special_case (svfloat32_t x, svfloat32_t y, svbool_t special)
 {
-  return sv_call_f32 (sinf, x, y, cmp);
+  special = svaclt (special, x, sv_f32 (INFINITY));
+
+  svfloat32x2_t reduction = large_range_reduction (svptrue_b32 (), x);
+
+  /* Unpack the quadrant from the return struct.  */
+  svuint32_t quadrant = svreinterpret_u32 (svget2 (reduction, 1));
+  svfloat32_t r = svget2 (reduction, 0);
+
+  svfloat32_t f = svtssel (r, quadrant);
+  svfloat32_t r2 = svtsmul (r, quadrant);
+  svfloat32_t sin = sv_f32 (0.0f);
+  sin = svtmad (sin, r2, 4);
+  sin = svtmad (sin, r2, 3);
+  sin = svtmad (sin, r2, 2);
+  sin = svtmad (sin, r2, 1);
+  sin = svtmad (sin, r2, 0);
+  sin = svmul_x (svptrue_b32 (), f, sin);
+
+  return svsel (special, sin, y);
 }
 
-/* A fast SVE implementation of sinf.
-   Maximum error: 1.89 ULPs.
-   This maximum error is achieved at multiple values in [-2^18, 2^18]
-   but one example is:
-   SV_NAME_F1 (sin)(0x1.9247a4p+0) got 0x1.fffff6p-1 want 0x1.fffffap-1.  */
+/* Vector version of sinf.
+   The maximum observed error is 1.44 + 0.5 ULP when |x| < 0x1p20.
+   _ZGVsMxv_sinf(0x1.4b0d9cp+13)
+    got 0x1.fc28cep-3
+   want 0x1.fc28d2p-3.
+   The special domain has a higher maximum error than the fast path:
+   The maximum observed error is 2.69 + 0.5 ULP when |x| >= 0x1p20.
+   _ZGVsMxv_sinf (0x1.be07aap+77)
+    got 0x1.ffe05ep-5
+   want 0x1.ffe058p-5.  */
 svfloat32_t SV_NAME_F1 (sin) (svfloat32_t x, const svbool_t pg)
 {
   const struct data *d = ptr_barrier (&data);
-
-  svfloat32_t ax = svabs_x (pg, x);
-  svuint32_t sign
-      = sveor_x (pg, svreinterpret_u32 (x), svreinterpret_u32 (ax));
-  svbool_t cmp = svcmpge (pg, svreinterpret_u32 (ax), RangeVal);
-
-  /* pi_vals are a quad-word of helper values - the first 3 elements contain
-     -pi in extended precision, the last contains 1 / pi.  */
-  svfloat32_t pi_vals = svld1rq (svptrue_b32 (), &d->negpi1);
-
-  /* n = rint(|x|/pi).  */
-  svfloat32_t n = svmla_lane (sv_f32 (d->shift), ax, pi_vals, 3);
-  svuint32_t odd = svlsl_x (pg, svreinterpret_u32 (n), 31);
-  n = svsub_x (pg, n, d->shift);
-
-  /* r = |x| - n*pi  (range reduction into -pi/2 .. pi/2).  */
-  svfloat32_t r;
-  r = svmla_lane (ax, n, pi_vals, 0);
-  r = svmla_lane (r, n, pi_vals, 1);
-  r = svmla_lane (r, n, pi_vals, 2);
-
-  /* sin(r) approx using a degree 9 polynomial from the Taylor series
-     expansion. Note that only the odd terms of this are non-zero.  */
-  svfloat32_t r2 = svmul_x (pg, r, r);
-  svfloat32_t y;
-  y = svmla_x (pg, C (2), r2, C (3));
-  y = svmla_x (pg, C (1), r2, y);
-  y = svmla_x (pg, C (0), r2, y);
-  y = svmla_x (pg, r, r, svmul_x (pg, y, r2));
-
-  /* sign = y^sign^odd.  */
-  sign = sveor_x (pg, sign, odd);
-
+  svbool_t ptrue = svptrue_b32 ();
+
+  /* Load some constants in quad-word chunks to minimise memory access.  */
+  svfloat32_t negpio2_and_invpio2 = svld1rq (ptrue, &d->neg_pio2_1);
+
+  /* n = rint(x/(pi/2)).  */
+  svfloat32_t q = svmla_lane (sv_f32 (d->shift), x, negpio2_and_invpio2, 3);
+  svfloat32_t n = svsub_x (ptrue, q, d->shift);
+
+  /* r = x - n*(pi/2)  (range reduction into -pi/4 .. pi/4).  */
+  svfloat32_t r = x;
+  r = svmla_lane (r, n, negpio2_and_invpio2, 0);
+  r = svmla_lane (r, n, negpio2_and_invpio2, 1);
+  r = svmla_lane (r, n, negpio2_and_invpio2, 2);
+
+  /* Final multiplicative factor: 1.0 or x depending on bit #0 of q.  */
+  svuint32_t q_u = svreinterpret_u32 (q);
+  svfloat32_t f = svtssel (r, q_u);
+
+  /* sin(r) poly approx.  */
+  svfloat32_t r2 = svtsmul (r, q_u);
+  svfloat32_t y = sv_f32 (0.0f);
+  y = svtmad (y, r2, 4);
+  y = svtmad (y, r2, 3);
+  y = svtmad (y, r2, 2);
+  y = svtmad (y, r2, 1);
+  y = svtmad (y, r2, 0);
+
+  svbool_t cmp = svacge (pg, x, sv_f32 (d->range_val));
   if (__glibc_unlikely (svptest_any (pg, cmp)))
-    return special_case (x,
-			 svreinterpret_f32 (sveor_x (
-			     svnot_z (pg, cmp), svreinterpret_u32 (y), sign)),
-			 cmp);
-  return svreinterpret_f32 (sveor_x (pg, svreinterpret_u32 (y), sign));
+    return special_case (x, svmul_x (ptrue, f, y), cmp);
+  /* Apply factor.  */
+  return svmul_x (ptrue, f, y);
 }
diff --git a/sysdeps/aarch64/fpu/sv_trigf_fallback.h b/sysdeps/aarch64/fpu/sv_trigf_fallback.h
new file mode 100644
index 0000000000..b8967eda59
--- /dev/null
+++ b/sysdeps/aarch64/fpu/sv_trigf_fallback.h
@@ -0,0 +1,184 @@
+/* Vectorised fallback for Single-Precision SVE trig functions.
+
+   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/>.  */
+
+/* Row i uses q = i - 3 and stores a 4-term binary32 expansion of
+   frac((2/pi) * 2^(8*q)), biased into [0.5, 1.5). With x_reduced exponent in
+   [32, 39], ph.hi then lands on a multiple of 2^8 and contributes no bits to
+   q mod 4. Tables are split into groups to improve gather performance.  */
+static const struct trigf_fallback_data
+{
+  float d0[16];
+  float d1[16];
+  float d2[16];
+  float d3[16];
+  float pio2;
+} trigf_fallback_data = {
+  .d0 = {
+    0x1p0f,          0x1.0000a2p0f, 0x1.00a2fap0f, 0x1.45f306p-1f,
+    0x1.f306dcp-1f, 0x1.06dc9cp-1f, 0x1.6e4e44p0f, 0x1.4e4416p0f,
+    0x1.44152ap0f,  0x1.1529fcp0f,  0x1.29fc28p0f, 0x1.f84ebp-1f,
+    0x1.2757d2p0f,  0x1.57d1f6p0f, 0x1.a3ea6ap-1f, 0x1.ea69bcp-1f,
+  },
+  .d1 = {
+    0x1.45f306p-25f, 0x1.f306dcp-25f,  -0x1.f246c6p-26f, 0x1.b9391p-26f,
+    0x1.391054p-26f, 0x1.1054a8p-26f,  0x1.529fc2p-28f,  -0x1.ac07b2p-25f,
+    -0x1.ec5418p-31f, 0x1.3abe9p-27f,  -0x1.505c16p-25f, -0x1.70565ap-27f,
+    -0x1.596448p-29f, -0x1.96447ep-25f, -0x1.11f924p-27f, -0x1.f924ecp-27f,
+  },
+  .d2 = {
+    0x1.b9391p-50f,  0x1.391054p-50f, -0x1.df56bp-51f,  0x1.529fc2p-52f,
+    0x1.4fe13ap-51f, -0x1.ec5418p-55f, 0x1.d5f47ep-54f, 0x1.5f47d4p-50f,
+    0x1.f534dep-56f, -0x1.596448p-53f, 0x1.a6ee06p-51f, 0x1.dc0db6p-52f,
+    0x1.b6c52cp-57f, -0x1.24eb54p-51f, -0x1.d6a66cp-52f, 0x1.5993c4p-52f,
+  },
+  .d3 = {
+    0x1.529fc2p-76f,  0x1.4fe13ap-75f,  -0x1.ec5418p-79f, 0x1.d5f47ep-78f,
+    0x1.7d1f54p-76f, 0x1.f534dep-80f,  -0x1.65912p-79f,  0x1.a6ee06p-75f,
+    -0x1.f924ecp-83f, 0x1.b6c52cp-81f, 0x1.b6295ap-76f,  0x1.4acc9ep-79f,
+    -0x1.9b0ef2p-82f, 0x1.93c43ap-76f, -0x1.de37ep-79f,  0x1.c821p-79f,
+  },
+  .pio2 = 0x1.921fb6p+0f,
+};
+
+/* Error-free multiplication using double-float computation via TwoProd.
+   hi is the rounded product, lo is the exact FMA residual.  */
+static inline svfloat32x2_t
+two_prod (svbool_t pg, svfloat32_t a, svfloat32_t b)
+{
+  svfloat32_t hi = svmul_x (pg, a, b);
+  svfloat32_t lo = svnmls_x (pg, hi, a, b);
+  return svcreate2 (hi, lo);
+}
+
+/* Error-free sum using double-float computation via FastTwoSum, which
+   requires |a| >= |b|. hi is the rounded sum and lo recovers the low-order
+   bits lost by that rounding.  */
+static inline svfloat32x2_t
+fast_two_sum (svbool_t pg, svfloat32_t a, svfloat32_t b)
+{
+  svfloat32_t hi = svadd_x (pg, a, b);
+  svfloat32_t t = svsub_x (pg, hi, a);
+  svfloat32_t lo = svsub_x (pg, b, t);
+  return svcreate2 (hi, lo);
+}
+
+/* Gather coefficients of 2/pi for the selected rows.  */
+static inline svfloat32x4_t
+load_datablock (svbool_t pg, svuint32_t idx,
+		const struct trigf_fallback_data *d)
+{
+  idx = svand_x (pg, idx, 15);
+
+  svfloat32_t d0 = svld1_gather_index (pg, d->d0, idx);
+  svfloat32_t d1 = svld1_gather_index (pg, d->d1, idx);
+  svfloat32_t d2 = svld1_gather_index (pg, d->d2, idx);
+  svfloat32_t d3 = svld1_gather_index (pg, d->d3, idx);
+
+  return svcreate4 (d0, d1, d2, d3);
+}
+
+/* Reduce x for |x| > 0x1p8 inputs, such that:
+    x = (q + y) * (pi / 2), with y in [-1/2, 1/2]
+
+   Returns a svfloat32x2_t struct containing:
+    remainder: The remainder after reduction
+    quadrant: Quadrant of x as an integer reinterpreted as a float for packing.
+
+   Designed to be used with the SVE trig instructions.  */
+static inline svfloat32x2_t
+large_range_reduction (svbool_t pg, svfloat32_t x)
+{
+  const struct trigf_fallback_data *d = ptr_barrier (&trigf_fallback_data);
+
+  /* First, x is reduced into the range of [2^32, 2^40), by directly adjusting
+     the exponent. This ensures the leading product contribute only multiples
+     of 2^8, so the useful bits of q mod 4 are entirely contained within the
+     lower product terms.  */
+  svuint32_t ix = svreinterpret_u32 (x);
+  svint32_t x_e_m32
+      = svsub_x (pg, svreinterpret_s32 (svlsr_x (pg, ix, 23)), 127 + 32);
+
+  /* We can then use the new exponent as an index for the 2/pi table.  */
+  svuint32_t idx
+      = svreinterpret_u32 (svadd_x (pg, svasr_x (pg, x_e_m32, 3), 3));
+  svfloat32x4_t datablock = load_datablock (pg, idx, d);
+
+  /* x_e_m32 has already been split into:
+       x_e_m32 = 8 * ROW + offset
+     where ROW selected the 2/pi row above.
+
+     We want to keep the offset (x_e_m32 mod 8), and use it to produce a new
+     exponent (32 + offset) so that x_reduced is within our intended [32, 39]
+     window.  */
+  svint32_t masked
+      = svreinterpret_s32 (svand_x (pg, svreinterpret_u32 (x_e_m32), 7));
+  svuint32_t new_exponent
+      = svreinterpret_u32 (svlsl_x (pg, svadd_x (pg, masked, 127 + 32), 23));
+
+  /* Finally, we get our reduced x value, by reinserting the new exponent into
+     the original input mantissa.  */
+  svuint32_t signed_mantissa = svand_x (pg, ix, 0x807fffff);
+  svfloat32_t x_reduced
+      = svreinterpret_f32 (svorr_x (pg, new_exponent, signed_mantissa));
+
+  /* We now use the reduced x to calculate x ~= (q + y) * (pi / 2).
+     First, we multiply x_reduced by the first three chunks of the 2/pi
+     table, using double-single arithmetic to maintain a high precision
+     intermediate.  */
+  svfloat32x2_t ph = two_prod (pg, x_reduced, svget4 (datablock, 0));
+  svfloat32x2_t pm = two_prod (pg, x_reduced, svget4 (datablock, 1));
+  svfloat32x2_t pl = two_prod (pg, x_reduced, svget4 (datablock, 2));
+
+  svfloat32_t ph_lo = svget2 (ph, 1);
+  svfloat32_t pm_hi = svget2 (pm, 0);
+  svfloat32_t pm_lo = svget2 (pm, 1);
+  svfloat32_t pl_hi = svget2 (pl, 0);
+  svfloat32_t pl_lo = svget2 (pl, 1);
+
+  /* Next, we need to accumulate the results together to get an integer k for
+     our quadrant. However, ph.hi will always be a multiple of 2^8, so it
+     cannot affect k mod 4. pm.lo and pl will always be sufficiently small
+     that they cannot affect the integer portion of the result. Therefore we
+     only need to sum ph.lo and pm.hi when computing k mod 4. Rounding sum_hi
+     chooses the nearest quadrant in pi/2 units, which is the control value
+     expected by the trig instructions.  */
+  svfloat32_t sum_hi = svadd_x (pg, ph_lo, pm_hi);
+  svfloat32_t kd = svrinta_x (pg, sum_hi);
+
+  /* To compute the remainder, we need to remove k from the pi/2-scaled
+     value and accumulate the remaining terms as a double-single remainder.  */
+  svfloat32_t y_hi = svadd_x (pg, svsub_x (pg, ph_lo, kd), pm_hi);
+  svfloat32x2_t y_mid = fast_two_sum (pg, pm_lo, pl_hi);
+
+  /* The low portion of x_reduced * D3 has no meaningful contribution to the
+    result, so a simple FMA is sufficient.  */
+  svfloat32_t y_lo = svmla_x (pg, pl_lo, x_reduced, svget4 (datablock, 3));
+
+  /* We then accumulate the final hi/lo remainders.  */
+  y_hi = svadd_x (pg, y_hi, svget2 (y_mid, 0));
+  y_lo = svadd_x (pg, y_lo, svget2 (y_mid, 1));
+
+  /* Multiply the accumulated remainders by pi/2, and adding gives a single
+     final remainder.  */
+  svfloat32_t remainder = svmla_x (pg, svmul_x (pg, y_lo, sv_f32 (d->pio2)),
+				   y_hi, sv_f32 (d->pio2));
+  svint32_t quadrant = svcvt_s32_x (pg, kd);
+  /* Reinterpret quadrant into a float to pack into struct for return.  */
+  return svcreate2 (remainder, svreinterpret_f32 (quadrant));
+}