[libc-commits] [libc] [libc][mathvec] Vectorise cosf (PR #224290)
via libc-commits
libc-commits at lists.llvm.org
Thu Sep 17 06:13:56 PDT 2026
https://github.com/DylanFleming-arm created https://github.com/llvm/llvm-project/pull/224290
Replaces loop over scalar cosf with a fully vectorised implementation.
>From e643276205983a81cfd42be20e9ab6d7dd930964 Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Thu, 17 Sep 2026 13:05:17 +0000
Subject: [PATCH] [libc][mathvec] Vectorise cosf
Replaces loop over scalar cosf with a fully vectorised implementation.
---
libc/src/__support/mathvec/CMakeLists.txt | 3 +-
libc/src/__support/mathvec/cosf.h | 124 +++++++++++++++++-
libc/src/__support/mathvec/trig_reductionf.h | 29 ++--
.../__support/mathvec/trig_reductionf_nofma.h | 33 +++--
4 files changed, 154 insertions(+), 35 deletions(-)
diff --git a/libc/src/__support/mathvec/CMakeLists.txt b/libc/src/__support/mathvec/CMakeLists.txt
index 431ad79f10501..dba7b6baba3d6 100644
--- a/libc/src/__support/mathvec/CMakeLists.txt
+++ b/libc/src/__support/mathvec/CMakeLists.txt
@@ -87,7 +87,8 @@ add_header_library(
cosf.h
DEPENDS
libc.src.__support.CPP.simd
- libc.src.__support.math.cosf
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.macros.properties.cpu_features
)
add_header_library(
diff --git a/libc/src/__support/mathvec/cosf.h b/libc/src/__support/mathvec/cosf.h
index b6ea0c848afd7..0d3cbe2922898 100644
--- a/libc/src/__support/mathvec/cosf.h
+++ b/libc/src/__support/mathvec/cosf.h
@@ -15,18 +15,132 @@
#define LLVM_LIBC_SRC___SUPPORT_MATHVEC_COSF_H
#include "src/__support/CPP/simd.h"
-#define LIBC_MATH_HAS_NO_ERRNO
-#define LIBC_MATH_HAS_NO_EXCEPT
-#define LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
-#include "src/__support/math/cosf.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/macros/properties/cpu_features.h"
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+#include "src/__support/mathvec/trig_reductionf.h"
+#else
+#include "src/__support/mathvec/trig_reductionf_nofma.h"
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
namespace LIBC_NAMESPACE_DECL {
namespace mathvec {
+template <size_t N>
+LIBC_INLINE static cpp::simd<double, N> cospif_poly(cpp::simd<double, N> r) {
+ // Approximate cos(pi * r) as (1/4 - r^2) * P(r^2) for |r| <= 0.5.
+ // These coefficients aren't produced directly via sollya, but rather
+ // are fine-tuned by iterative adjustment to remove hard to round cases.
+ // TODO: Create a tool to deterministically reproduce these coefficients.
+ // see https://github.com/llvm/llvm-project/issues/220984
+ constexpr cpp::simd<double, N> c0 = 0x1p2;
+ constexpr cpp::simd<double, N> c1 = -0x1.de9e64df22f07p1;
+ constexpr cpp::simd<double, N> c2 = 0x1.472be1223eeadp0;
+ constexpr cpp::simd<double, N> c3 = -0x1.d4fcd82b511ebp-3;
+ constexpr cpp::simd<double, N> c4 = 0x1.9f05c866286ddp-6;
+ constexpr cpp::simd<double, N> c5 = -0x1.f308c07837812p-10;
+ constexpr cpp::simd<double, N> c6 = 0x1.b22004be59c73p-14;
+ constexpr cpp::simd<double, N> c7 = -0x1.14bd54572305fp-18;
+ constexpr cpp::simd<double, N> f = 0.25;
+
+ cpp::simd<double, N> r2 = r * r;
+ cpp::simd<double, N> r4 = r2 * r2;
+ cpp::simd<double, N> p01 = cpp::multiply_add(r2, c1, c0);
+ cpp::simd<double, N> p23 = cpp::multiply_add(r2, c3, c2);
+ cpp::simd<double, N> p45 = cpp::multiply_add(r2, c5, c4);
+ cpp::simd<double, N> p67 = cpp::multiply_add(r2, c7, c6);
+ cpp::simd<double, N> p47 = cpp::multiply_add(r4, p67, p45);
+ cpp::simd<double, N> p27 = cpp::multiply_add(r4, p47, p23);
+ cpp::simd<double, N> p07 = cpp::multiply_add(r4, p27, p01);
+
+ cpp::simd<double, N> factor = cpp::multiply_add(-r, r, f);
+
+ return factor * p07;
+}
+
+// Correct specific cases which aren't able to correctly round from the normal
+// codepath.
+template <size_t N>
+LIBC_INLINE static cpp::simd<float, N>
+repair_hard_to_round(cpp::simd<float, N> x, cpp::simd<float, N> y) {
+ y = (x == 0x1.fcd9eep+38f) ? 0x1.3371d2p-17f : y;
+ y = (x == 0x1.3170f0p+63f) ? 0x1.fe2976p-1f : y;
+ y = (x == 0x1.15313ep+69f) ? 0x1.f52484p-20f : y;
+
+#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ y = (x == 0x1.8db252p+25f) ? -0x1.527a0ap-5f : y;
+ y = (x == 0x1.aff73cp+32f) ? -0x1.d59e7ap-23f : y;
+ y = (x == 0x1.47d0fep+34f) ? -0x1.149dbp-29f : y;
+ y = (x == 0x1.455500p+51f) ? 0x1.115d7ep-1f : y;
+ y = (x == 0x1.407c74p+100f) ? 0x1.a4a122p-19f : y;
+#endif
+
+ return y;
+}
+
+// Due to specific lane correction being expensive, we can perform a fast
+// estimated check to determine if we should branch to specific hard to round
+// case correction. Will return true for all hard to round cases, but can
+// produce false positives. False positive rate in [0x1.90bdfap+20f, inf]:
+// With FMA: 52/901,226,755 ~= 1/2^24
+// Without FMA: 7,040,869/901,226,755 ~= 1/2^7
+template <size_t N>
+LIBC_INLINE static cpp::simd<bool, N>
+is_maybe_hard_to_round(cpp::simd<float, N> x) {
+ cpp::simd<uint32_t, N> x_bits = cpp::bit_cast<cpp::simd<uint32_t>>(x);
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ return ((x_bits * 0xd6c4de5f) & 0xefffda56) == 0x46048000;
+#else
+ return ((x_bits * 0x30d04f5f) & 0x8b002240) == 0x8002000;
+#endif
+}
+
template <size_t N>
LIBC_INLINE cpp::simd<float, N> cosf(cpp::simd<float, N> x) {
- return cpp::map(x, [](float a) { return math::cosf(a); });
+ using FPBits = typename fputil::FPBits<float>;
+ cpp::simd<double, N> x_d = cpp::simd_cast<double>(x);
+ cpp::simd<float, N> ax = cpp::abs(x);
+
+ // If all lanes are below pi/2, we can skip the reduction entirely.
+ cpp::simd<bool, N> is_small = ax <= 0x1.921fb6p+0f;
+ if (cpp::all_of(is_small)) {
+ constexpr cpp::simd<double, N> inv_pi = 0x1.45f306dc9c883p-2;
+ return cpp::simd_cast<float>(cospif_poly(x_d * inv_pi));
+ }
+
+ // Computes the main reduction pass
+ Reduction<N> reduce = fast_reduction(x_d);
+
+ // Large inputs require a more involved reduction, as well as inf handling.
+ cpp::simd<bool, N> has_large_reduction = ax > 0x1.90bdfap+20f;
+ if (LIBC_UNLIKELY(cpp::any_of(has_large_reduction))) {
+ cpp::simd<bool, N> is_finite = ax < FPBits::inf().get_val();
+ Reduction<N> large_reduce = large_reduction(x_d);
+ reduce.r = has_large_reduction ? large_reduce.r : reduce.r;
+ reduce.r = is_finite ? reduce.r : FPBits::quiet_nan().get_val();
+ reduce.odd = has_large_reduction ? large_reduce.odd : reduce.odd;
+ }
+
+ // Both reduction paths feed into a single polynomial evaluation + sign
+ // correction.
+ cpp::simd<float, N> poly = cpp::simd_cast<float>(cospif_poly(reduce.r));
+ cpp::simd<uint32_t, N> sign = cpp::simd_cast<uint32_t>(reduce.odd) << 31;
+
+ cpp::simd<float, N> y = cpp::bit_cast<cpp::simd<float>>(
+ cpp::bit_cast<cpp::simd<uint32_t>>(poly) ^ sign);
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ y = (ax == 0x1.d2eb54p+10) ? -0x1.633eb4p-13 : y;
+#endif
+ if (LIBC_UNLIKELY(cpp::any_of(has_large_reduction))) {
+ if (LIBC_UNLIKELY(cpp::any_of(is_maybe_hard_to_round(ax))))
+ return repair_hard_to_round(ax, y);
+ }
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return y;
}
} // namespace mathvec
diff --git a/libc/src/__support/mathvec/trig_reductionf.h b/libc/src/__support/mathvec/trig_reductionf.h
index dac2506683400..fe976bc972161 100644
--- a/libc/src/__support/mathvec/trig_reductionf.h
+++ b/libc/src/__support/mathvec/trig_reductionf.h
@@ -40,19 +40,19 @@ LIBC_INLINE static Reduction<N> fast_reduction(cpp::simd<double, N> x) {
cpp::simd<double, N> t = shift - z;
// r = x/pi - k
- cpp::simd<double, N> r;
- r = cpp::multiply_add(x, inv_pi, t);
+ cpp::simd<double, N> r = cpp::multiply_add(x, inv_pi, t);
r = cpp::multiply_add(x, inv_pi_tail, r);
return {r, cpp::bit_cast<cpp::simd<int64_t, N>>(z)};
}
-// Two-double expansions of 2^(8*q) / pi reduced modulo an even integer,
-// q = 3..12. Padded to a 16 length array for easier access.
+// Two-double expansions of 2^(8*q) / pi reduced modulo an even integer.
+// Entries 0..14 correspond to q = 0..14, and entry 15 corresponds to q = -1
+// after masking. Entry 13 and 14 act as padding for inactive lanes.
LIBC_INLINE_VAR constexpr double INV_PI_HI[16] = {
- 0,
- 0,
- 0,
+ 0x1.45f306dc9c883p-2,
+ -0x1.067c91b1bbeadp-1,
+ 0x1.836e4e44152ap-1,
-0x1.236377d5ac07bp-2,
-0x1.b1bbead603d8bp-1,
-0x1.bbead603d8a83p-1,
@@ -65,13 +65,13 @@ LIBC_INLINE_VAR constexpr double INV_PI_HI[16] = {
0x1.f534ddc0db629p-1,
0,
0,
- 0,
+ 0x1.45f306dc9c883p-10,
};
LIBC_INLINE_VAR constexpr double INV_PI_LO[16] = {
- 0,
- 0,
- 0,
+ -0x1.6b01ec55ff4d6p-56,
+ -0x1.80f62a0b82b2dp-55,
+ -0x1.ec54170565912p-56,
-0x1.505c1596447e5p-58,
0x1.f47d4d377036ep-55,
0x1.f534ddc0db629p-57,
@@ -84,7 +84,7 @@ LIBC_INLINE_VAR constexpr double INV_PI_LO[16] = {
0x1.664f10e4107f9p-55,
0,
0,
- 0,
+ -0x1.6b01e7abff3d6p-64,
};
// Reduces non-negative large finite inputs x >= 0x1p49.
@@ -100,9 +100,8 @@ LIBC_INLINE static Reduction<N> large_reduction(cpp::simd<double, N> x) {
cpp::simd<int64_t, N> q =
cpp::simd_cast<int64_t>(ix >> 55) - cpp::simd<int64_t, N>(131);
- // While sufficiently large x will always produce q within [3, 12],
- // not all input lanes are guaranteed to require the large reduction,
- // so we mask with 15 to keep all values within bounds.
+ // For float inputs requiring large reduction, q is in [-1, 12]. Masking
+ // maps q = -1 to entry 15 and keeps inactive lanes within bounds.
cpp::simd<int64_t, N> idx = q & cpp::simd<int64_t, N>(15);
cpp::simd<double, N> c_hi =
cpp::gather<cpp::simd<double, N>>(true, idx, INV_PI_HI);
diff --git a/libc/src/__support/mathvec/trig_reductionf_nofma.h b/libc/src/__support/mathvec/trig_reductionf_nofma.h
index ac2f029922c5b..9d4ba0fa57d4a 100644
--- a/libc/src/__support/mathvec/trig_reductionf_nofma.h
+++ b/libc/src/__support/mathvec/trig_reductionf_nofma.h
@@ -55,12 +55,15 @@ LIBC_INLINE static Reduction<N> fast_reduction(cpp::simd<double, N> x) {
return {r, cpp::bit_cast<cpp::simd<int64_t, N>>(z)};
}
-// Three-double expansions of 2^(8*q) / pi, reduced modulo an even integer,
-// q = 3..12. Padded to 16 length arrays for safe masked indexing.
+// Three-double expansions of 2^(8*q) / pi reduced modulo an even integer.
+// Entries 0..14 correspond to q = 0..14, and entry 15 corresponds to q = -1
+// after masking. Entry 13 and 14 act as padding for inactive lanes.
+// The first two parts have 29 significant bits, such that the fp64 product
+// is exact when multiplied by an fp32 input.
LIBC_INLINE_VAR constexpr double INV_PI_HI[16] = {
- 0,
- 0,
- 0,
+ 0x1.45f306ep-2,
+ -0x1.067c91bp-1,
+ 0x1.836e4e4p-1,
-0x1.236377d000000p-2,
-0x1.b1bbead000000p-1,
-0x1.bbead60000000p-1,
@@ -73,12 +76,13 @@ LIBC_INLINE_VAR constexpr double INV_PI_HI[16] = {
0x1.f534ddc000000p-1,
0,
0,
+ 0x1.45f306ep-10,
};
LIBC_INLINE_VAR constexpr double INV_PI_MID[16] = {
- 0,
- 0,
- 0,
+ -0x1.b1bbeadp-33,
+ -0x1.bbead6p-33,
+ 0x1.054a7f1p-31,
-0x1.6b01ec5000000p-32,
-0x1.80f62a1000000p-31,
-0x1.ec54170000000p-32,
@@ -91,12 +95,13 @@ LIBC_INLINE_VAR constexpr double INV_PI_MID[16] = {
0x1.b6c52b3000000p-34,
0,
0,
+ -0x1.b1bbeadp-41,
};
LIBC_INLINE_VAR constexpr double INV_PI_LO[16] = {
- 0,
- 0,
- 0,
+ -0x1.80f62a0b82b2dp-63,
+ -0x1.ec54170565912p-64,
+ -0x1.8a82e0acb223fp-61,
-0x1.05c1596447e50p-62,
0x1.1f534ddc0db80p-61,
-0x1.596447e493ae0p-62,
@@ -109,6 +114,7 @@ LIBC_INLINE_VAR constexpr double INV_PI_LO[16] = {
0x1.3c439041fe400p-65,
0,
0,
+ -0x1.80f62a0b82b2dp-71,
};
// Reduces non-negative large finite inputs x >= 0x1p49.
@@ -124,9 +130,8 @@ LIBC_INLINE static Reduction<N> large_reduction(cpp::simd<double, N> x) {
cpp::simd<int64_t, N> q =
cpp::simd_cast<int64_t>(ix >> 55) - cpp::simd<int64_t, N>(131);
- // While sufficiently large x will always produce q within [3, 12],
- // not all input lanes are guaranteed to require the large reduction,
- // so we mask with 15 to keep all values within bounds.
+ // For float inputs requiring large reduction, q is in [-1, 12]. Masking
+ // maps q = -1 to entry 15 and keeps inactive lanes within bounds.
cpp::simd<int64_t, N> idx = q & cpp::simd<int64_t, N>(15);
cpp::simd<double, N> c_hi =
cpp::gather<cpp::simd<double, N>>(true, idx, INV_PI_HI);
More information about the libc-commits
mailing list