[libc-commits] [libc] 962d16f - [libc][mathvec] Vectorise exp2f and exp10f (#211365)
via libc-commits
libc-commits at lists.llvm.org
Tue Aug 4 07:03:47 PDT 2026
Author: DylanFleming-arm
Date: 2026-08-04T15:03:42+01:00
New Revision: 962d16fd97793d22c42227e7fab397f62475c18b
URL: https://github.com/llvm/llvm-project/commit/962d16fd97793d22c42227e7fab397f62475c18b
DIFF: https://github.com/llvm/llvm-project/commit/962d16fd97793d22c42227e7fab397f62475c18b.diff
LOG: [libc][mathvec] Vectorise exp2f and exp10f (#211365)
Replaces exp2f and exp10f with fully vectorised implementations.
Includes refactoring of expf polynomial evalulation into expf_utils.h
Added:
Modified:
libc/src/__support/CPP/CMakeLists.txt
libc/src/__support/CPP/simd.h
libc/src/__support/mathvec/CMakeLists.txt
libc/src/__support/mathvec/exp10f.h
libc/src/__support/mathvec/exp2f.h
libc/src/__support/mathvec/expf.h
libc/src/__support/mathvec/expf_utils.h
Removed:
################################################################################
diff --git a/libc/src/__support/CPP/CMakeLists.txt b/libc/src/__support/CPP/CMakeLists.txt
index e11d823333b25..8256caebfe754 100644
--- a/libc/src/__support/CPP/CMakeLists.txt
+++ b/libc/src/__support/CPP/CMakeLists.txt
@@ -246,6 +246,7 @@ add_header_library(
simd.h
DEPENDS
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
.utility
.tuple
)
diff --git a/libc/src/__support/CPP/simd.h b/libc/src/__support/CPP/simd.h
index 93fecc71d36d0..0a7f1eefbd438 100644
--- a/libc/src/__support/CPP/simd.h
+++ b/libc/src/__support/CPP/simd.h
@@ -22,6 +22,7 @@
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
#include <stddef.h>
@@ -208,9 +209,22 @@ LIBC_INLINE constexpr static simd<T, N> abs(simd<T, N> x) {
return __builtin_elementwise_abs(x);
}
template <typename T, size_t N>
-LIBC_INLINE constexpr static simd<T, N> fma(simd<T, N> x, simd<T, N> y,
- simd<T, N> z) {
- return __builtin_elementwise_fma(x, y, z);
+LIBC_INLINE constexpr static simd<T, N> multiply_add(simd<T, N> x, simd<T, N> y,
+ simd<T, N> z) {
+#if __has_builtin(__builtin_elementwise_fma)
+#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
+ if constexpr (cpp::is_same_v<T, float>)
+ return __builtin_elementwise_fma(x, y, z);
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ if constexpr (cpp::is_same_v<T, double>)
+ return __builtin_elementwise_fma(x, y, z);
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+#endif // __has_builtin(__builtin_elementwise_fma)
+
+ simd<T, N> mul = x * y;
+ return mul + z;
}
template <typename T, size_t N>
LIBC_INLINE constexpr static simd<T, N> ceil(simd<T, N> x) {
diff --git a/libc/src/__support/mathvec/CMakeLists.txt b/libc/src/__support/mathvec/CMakeLists.txt
index 5fcb2feb2f4a8..6e5f3bbddc1be 100644
--- a/libc/src/__support/mathvec/CMakeLists.txt
+++ b/libc/src/__support/mathvec/CMakeLists.txt
@@ -122,8 +122,10 @@ add_header_library(
HDRS
exp2f.h
DEPENDS
+ .expf_utils
libc.src.__support.CPP.simd
- libc.src.__support.math.exp2f
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.common
)
add_header_library(
@@ -140,8 +142,10 @@ add_header_library(
HDRS
exp10f.h
DEPENDS
+ .expf_utils
libc.src.__support.CPP.simd
- libc.src.__support.math.exp10f
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.common
)
add_header_library(
@@ -289,4 +293,4 @@ add_header_library(
DEPENDS
libc.src.__support.CPP.simd
libc.src.__support.math.tanpif
-)
\ No newline at end of file
+)
diff --git a/libc/src/__support/mathvec/exp10f.h b/libc/src/__support/mathvec/exp10f.h
index 25a62fcbda9e4..d22baf16b40ac 100644
--- a/libc/src/__support/mathvec/exp10f.h
+++ b/libc/src/__support/mathvec/exp10f.h
@@ -14,19 +14,55 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP10F_H
#define LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP10F_H
+#include "expf_utils.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/exp10f.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/common.h"
namespace LIBC_NAMESPACE_DECL {
namespace mathvec {
+template <size_t N>
+LIBC_INLINE static cpp::simd<double, N> inline_exp10(cpp::simd<double, N> x) {
+ constexpr cpp::simd<double, N> shift = 0x1.800000000ffc0p+46;
+ constexpr cpp::simd<double, N> log2_10 = 0x1.a934f0979a371p+1;
+
+ // Rounds n to be (x * log2(10)) to the nearest multiple of 1/64
+ // While preparing z to be used as an index for the lookup in eval_exp
+ cpp::simd<double, N> z = cpp::multiply_add(x, log2_10, shift);
+ cpp::simd<double, N> n = z - shift;
+
+ constexpr cpp::simd<double, N> log2_10_lo = 0x1.7f2495fb7fa6dp-53;
+ constexpr cpp::simd<double, N> negative_ln2 = -0x1.62e42fefa39efp-1;
+
+ // Using 10^x = e^(x * ln(10)) = e^(x * log2(10) * ln(2))
+ // Reduces x into r as:
+ // r = (x * log2(10) - n) * ln(2), with |r| <= ln2/128,
+ // Computed as (n - x * log2(10)) * -ln(2) to avoid negating n.
+ cpp::simd<double, N> r;
+ r = cpp::multiply_add(-x, log2_10, n);
+ r = cpp::multiply_add(-x, log2_10_lo, r);
+ r = r * negative_ln2;
+
+ return eval_exp(r, z);
+}
+
template <size_t N>
LIBC_INLINE cpp::simd<float, N> exp10f(cpp::simd<float, N> x) {
- return cpp::map(x, [](float a) { return math::exp10f(a); });
+ using FPBits = typename fputil::FPBits<float>;
+
+ cpp::simd<bool, N> is_inf = x >= 0x1.344136p+5f;
+ cpp::simd<bool, N> is_zero = x <= -0x1.693c6bp+5f;
+ cpp::simd<bool, N> is_special = is_inf | is_zero;
+
+ cpp::simd<float, N> special_res = is_inf ? FPBits::inf().get_val() : 0.0f;
+
+ cpp::simd<double, N> x_d = cpp::simd_cast<double, float, N>(x);
+ cpp::simd<double, N> y = inline_exp10(x_d);
+ cpp::simd<float, N> ret = cpp::simd_cast<float, double, N>(y);
+
+ return is_special ? special_res : ret;
}
} // namespace mathvec
diff --git a/libc/src/__support/mathvec/exp2f.h b/libc/src/__support/mathvec/exp2f.h
index fb9c79890a726..848b9647dd612 100644
--- a/libc/src/__support/mathvec/exp2f.h
+++ b/libc/src/__support/mathvec/exp2f.h
@@ -14,19 +14,51 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP2F_H
#define LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP2F_H
+#include "expf_utils.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/exp2f.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/common.h"
namespace LIBC_NAMESPACE_DECL {
namespace mathvec {
+template <size_t N>
+LIBC_INLINE static cpp::simd<double, N> inline_exp2(cpp::simd<double, N> x) {
+ constexpr cpp::simd<double, N> shift = 0x1.800000000ffc0p+46;
+
+ // Rounds n to be x to the nearest multiple of 1/64
+ // While preparing z to be used as an index for the lookup in eval_exp
+ cpp::simd<double, N> z = x + shift;
+ cpp::simd<double, N> n = z - shift;
+
+ // ln(2), with a bias of +768 FP64 ULP, which is enough to remove the
+ // hard to round cases when casting back to FP32.
+ constexpr cpp::simd<double, N> biased_ln2 = 0x1.62e42fefa3cefp-1;
+
+ // Using 2^x = e^(x * ln2),
+ // Reduces x into r as:
+ // r = (x - n) * ln(2), with |r| <= ln2/128.
+ // then computes 2^x = 2^n * exp(r).
+ cpp::simd<double, N> r = (x - n) * biased_ln2;
+ return eval_exp(r, z);
+}
+
template <size_t N>
LIBC_INLINE cpp::simd<float, N> exp2f(cpp::simd<float, N> x) {
- return cpp::map(x, [](float a) { return math::exp2f(a); });
+ using FPBits = typename fputil::FPBits<float>;
+
+ cpp::simd<bool, N> is_inf = x >= 0x1p7;
+ cpp::simd<bool, N> is_zero = x <= -0x1.2cp7f;
+ cpp::simd<bool, N> is_special = is_inf | is_zero;
+
+ cpp::simd<float, N> special_res = is_inf ? FPBits::inf().get_val() : 0.0f;
+
+ cpp::simd<double, N> x_d = cpp::simd_cast<double, float, N>(x);
+ cpp::simd<double, N> y = inline_exp2(x_d);
+ cpp::simd<float, N> ret = cpp::simd_cast<float, double, N>(y);
+
+ return is_special ? special_res : ret;
}
} // namespace mathvec
diff --git a/libc/src/__support/mathvec/expf.h b/libc/src/__support/mathvec/expf.h
index c27de2527e5fa..95f5a5dc76874 100644
--- a/libc/src/__support/mathvec/expf.h
+++ b/libc/src/__support/mathvec/expf.h
@@ -28,7 +28,7 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp(cpp::simd<double, N> x) {
// inv_ln2 = round(1/log(2), D, RN);
constexpr cpp::simd<double, N> inv_ln2 = 0x1.71547652b82fep+0;
- cpp::simd<double, N> z = shift + x * inv_ln2;
+ cpp::simd<double, N> z = cpp::multiply_add(x, inv_ln2, shift);
cpp::simd<double, N> n = z - shift;
// ln2_hi = round(log(2), D, RN);
@@ -37,40 +37,18 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp(cpp::simd<double, N> x) {
constexpr cpp::simd<double, N> ln2_lo = 0x1.abc9e3b39803fp-56;
cpp::simd<double, N> r = x;
- r = r - n * ln2_hi;
- r = r - n * ln2_lo;
-
- // Coefficients of exp approximation, generated by Sollya with:
- // poly = 1 + x;
- // for i from 2 to 5 do {
- // r = remez(exp(x)-poly(x), 5-i, [-log(2)/128;log(2)/128], x^i, 1e-10);
- // c = coeff(roundcoefficients(r, [|D ...|]), 0);
- // poly = poly + x^i*c;
- // c;
- // };
- constexpr cpp::simd<double, N> c0 = 0x1.fffffffffdbcep-2;
- constexpr cpp::simd<double, N> c1 = 0x1.55555555543c2p-3;
- constexpr cpp::simd<double, N> c2 = 0x1.555573c64f2e3p-5;
- constexpr cpp::simd<double, N> c3 = 0x1.111126b4eff73p-7;
-
- /* y = exp(r) - 1 ~= r + C0 r^2 + C1 r^3 + C2 r^4 + C3 r^5. */
- cpp::simd<double, N> r2 = r * r;
- cpp::simd<double, N> p01 = c0 + r * c1;
- cpp::simd<double, N> p23 = c2 + r * c3;
- cpp::simd<double, N> p04 = p01 + r2 * p23;
- cpp::simd<double, N> y = r + p04 * r2;
-
- cpp::simd<uint64_t, N> u = cpp::bit_cast<cpp::simd<uint64_t, N>>(z);
- cpp::simd<double, N> s = exp_lookup(u);
- return s + s * y;
+ r = cpp::multiply_add(-n, ln2_hi, r);
+ r = cpp::multiply_add(-n, ln2_lo, r);
+
+ return eval_exp(r, z);
}
template <size_t N>
LIBC_INLINE cpp::simd<float, N> expf(cpp::simd<float, N> x) {
using FPBits = typename fputil::FPBits<float>;
- cpp::simd<bool, N> is_inf = x >= 0x1.62e38p+9;
- cpp::simd<bool, N> is_zero = x <= -0x1.628c2ap+9;
+ cpp::simd<bool, N> is_inf = x >= 0x1.62e43p+6f;
+ cpp::simd<bool, N> is_zero = x <= -0x1.9fe36ap+6f;
cpp::simd<bool, N> is_special = is_inf | is_zero;
cpp::simd<float, N> special_res = is_inf ? FPBits::inf().get_val() : 0.0f;
diff --git a/libc/src/__support/mathvec/expf_utils.h b/libc/src/__support/mathvec/expf_utils.h
index c23e286f6563e..8ce3d3c3b5e6e 100644
--- a/libc/src/__support/mathvec/expf_utils.h
+++ b/libc/src/__support/mathvec/expf_utils.h
@@ -7,7 +7,7 @@
//===----------------------------------------------------------------------===//
///
/// \file
-/// This file contains utility functions for single-precision SIMD exp.
+/// This file contains utility functions for single-precision SIMD exp/2/10.
///
//===----------------------------------------------------------------------===//
@@ -32,6 +32,39 @@ LIBC_INLINE static cpp::simd<double, N> exp_lookup(cpp::simd<uint64_t, N> u) {
return cpp::bit_cast<cpp::simd<double, N>>(result);
}
+template <size_t N>
+LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
+ cpp::simd<double, N> z) {
+ // Coefficients of exp approximation, generated by Sollya with:
+ // poly = 1 + x;
+ // for i from 2 to 5 do {
+ // r = remez(exp(x)-poly(x), 5-i, [-log(2)/128;log(2)/128], x^i, 1e-10);
+ // c = coeff(roundcoefficients(r, [|D ...|]), 0);
+ // poly = poly + x^i*c;
+ // c;
+ // };
+ constexpr cpp::simd<double, N> c0 = 0x1.fffffffffdbcep-2;
+ constexpr cpp::simd<double, N> c1 = 0x1.55555555543c2p-3;
+ constexpr cpp::simd<double, N> c2 = 0x1.555573c64f2e3p-5;
+ constexpr cpp::simd<double, N> c3 = 0x1.111126b4eff73p-7;
+
+ // y = exp(r) - 1 ~= r + C0 r^2 + C1 r^3 + C2 r^4 + C3 r^5.
+ cpp::simd<double, N> r2 = r * r;
+ cpp::simd<double, N> p01 = cpp::multiply_add(c1, r, c0);
+ cpp::simd<double, N> p23 = cpp::multiply_add(c3, r, c2);
+ cpp::simd<double, N> p04 = cpp::multiply_add(p23, r2, p01);
+ cpp::simd<double, N> y = cpp::multiply_add(p04, r2, r);
+
+ // Table lookup for 2^n, where n is a multiple of 1/64
+ cpp::simd<uint64_t, N> u = cpp::bit_cast<cpp::simd<uint64_t, N>>(z);
+ cpp::simd<double, N> s = exp_lookup(u);
+
+ // e^x = 2^n * exp(r)
+ // Since y = exp(r) - 1, e^x = 2^n * (1 + y)
+ // Or as an FMA: e^x = 2^n + (2^n * y)
+ return cpp::multiply_add(y, s, s);
+}
+
} // namespace mathvec
} // namespace LIBC_NAMESPACE_DECL
More information about the libc-commits
mailing list