[PATCH 2/4] aarch64: Vectorise SVE fp32 sin/cos special case
Dylan Fleming
Dylan.Fleming@arm.com
Thu Jun 18 11:07:06 GMT 2026
Implemented single precision vector fallback for SVE sinf/cosf.
Includes a rewrite of SVE sinf to use the trig instructions.
---
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(-)
create mode 100644 sysdeps/aarch64/fpu/sv_trigf_fallback.h
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));
+}
--
2.43.0
More information about the Libc-alpha
mailing list