[libc-commits] [libc] [libc][mathvec] Vectorise exp2f and exp10f (PR #211365)
via libc-commits
libc-commits at lists.llvm.org
Tue Aug 4 05:59:54 PDT 2026
https://github.com/DylanFleming-arm updated https://github.com/llvm/llvm-project/pull/211365
>From ee2fbd23fa00fc686475dc7b5c8d8214f28cad2c Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Fri, 17 Jul 2026 05:41:14 +0000
Subject: [PATCH 1/6] [libc][mathvec] Vectorise exp2f and exp10f
Replaces exp2f and exp10f with fully vectorised implementations.
Includes refactoring of expf polynomial evalulation into expf_utils.h
---
libc/src/__support/mathvec/CMakeLists.txt | 8 +++-
libc/src/__support/mathvec/exp10f.h | 46 ++++++++++++++++++++---
libc/src/__support/mathvec/exp2f.h | 42 ++++++++++++++++++---
libc/src/__support/mathvec/expf.h | 24 +-----------
libc/src/__support/mathvec/expf_utils.h | 35 ++++++++++++++++-
5 files changed, 119 insertions(+), 36 deletions(-)
diff --git a/libc/src/__support/mathvec/CMakeLists.txt b/libc/src/__support/mathvec/CMakeLists.txt
index 5fcb2feb2f4a8..9d04ce02f998a 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(
diff --git a/libc/src/__support/mathvec/exp10f.h b/libc/src/__support/mathvec/exp10f.h
index 25a62fcbda9e4..2414190e5fc07 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 = (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 = n - (x * log2_10);
+ r = r - (x * log2_10_lo);
+ 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 > 0x1p+7f;
+ cpp::simd<bool, N> is_zero = x < -0x1p+7f;
+ 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..22cdc4a6ac145 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 > 0x1p+9f;
+ cpp::simd<bool, N> is_zero = x < -0x1p+9f;
+ 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..3df1cf5f1108b 100644
--- a/libc/src/__support/mathvec/expf.h
+++ b/libc/src/__support/mathvec/expf.h
@@ -40,29 +40,7 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp(cpp::simd<double, N> 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;
+ return eval_exp(r, z);
}
template <size_t N>
diff --git a/libc/src/__support/mathvec/expf_utils.h b/libc/src/__support/mathvec/expf_utils.h
index c23e286f6563e..d00ac1282b245 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 = 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;
+
+ // 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 s + s * y;
+}
+
} // namespace mathvec
} // namespace LIBC_NAMESPACE_DECL
>From 0a1014e7d242610d091377f91e0225b70b4761a7 Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Thu, 23 Jul 2026 00:39:08 +0000
Subject: [PATCH 2/6] Refactored expf_utils.h to use multiply_add builtins
---
libc/src/__support/CPP/CMakeLists.txt | 2 ++
libc/src/__support/CPP/simd.h | 28 +++++++++++++++++++----
libc/src/__support/mathvec/CMakeLists.txt | 2 +-
libc/src/__support/mathvec/expf_utils.h | 10 ++++----
4 files changed, 31 insertions(+), 11 deletions(-)
diff --git a/libc/src/__support/CPP/CMakeLists.txt b/libc/src/__support/CPP/CMakeLists.txt
index 49cb07f329111..3732595fe212f 100644
--- a/libc/src/__support/CPP/CMakeLists.txt
+++ b/libc/src/__support/CPP/CMakeLists.txt
@@ -238,7 +238,9 @@ add_header_library(
HDRS
simd.h
DEPENDS
+ libc.src.__support.FPUtil.multiply_add
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 88d1ea4091a4e..ea28d39d1c84c 100644
--- a/libc/src/__support/CPP/simd.h
+++ b/libc/src/__support/CPP/simd.h
@@ -19,9 +19,11 @@
#include "src/__support/CPP/tuple.h"
#include "src/__support/CPP/type_traits.h"
#include "src/__support/CPP/utility/integer_sequence.h"
+#include "src/__support/FPUtil/multiply_add.h"
#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,11 +210,6 @@ 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);
-}
-template <typename T, size_t N>
LIBC_INLINE constexpr static simd<T, N> ceil(simd<T, N> x) {
return __builtin_elementwise_ceil(x);
}
@@ -415,6 +412,27 @@ LIBC_INLINE static cpp::simd<T, N> map(cpp::simd<T, N> v, F f) {
// TODO: where expressions, scalar overloads, ABI types.
} // namespace cpp
+
+namespace fputil {
+#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
+template <size_t N>
+LIBC_INLINE constexpr cpp::simd<float, N> multiply_add(cpp::simd<float, N> x,
+ cpp::simd<float, N> y,
+ cpp::simd<float, N> z) {
+ return __builtin_elementwise_fma(x, y, z);
+}
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+template <size_t N>
+LIBC_INLINE constexpr cpp::simd<double, N>
+multiply_add(cpp::simd<double, N> x, cpp::simd<double, N> y,
+ cpp::simd<double, N> z) {
+ return __builtin_elementwise_fma(x, y, z);
+}
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+} // namespace fputil
+
} // namespace LIBC_NAMESPACE_DECL
#endif // LIBC_HAS_VECTOR_TYPE
diff --git a/libc/src/__support/mathvec/CMakeLists.txt b/libc/src/__support/mathvec/CMakeLists.txt
index 9d04ce02f998a..6e5f3bbddc1be 100644
--- a/libc/src/__support/mathvec/CMakeLists.txt
+++ b/libc/src/__support/mathvec/CMakeLists.txt
@@ -293,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/expf_utils.h b/libc/src/__support/mathvec/expf_utils.h
index d00ac1282b245..cd674d9ded381 100644
--- a/libc/src/__support/mathvec/expf_utils.h
+++ b/libc/src/__support/mathvec/expf_utils.h
@@ -50,10 +50,10 @@ LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
// 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<double, N> p01 = fputil::multiply_add(c1, r, c0);
+ cpp::simd<double, N> p23 = fputil::multiply_add(c3, r, c2);
+ cpp::simd<double, N> p04 = fputil::multiply_add(p23, r2, p01);
+ cpp::simd<double, N> y = fputil::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);
@@ -62,7 +62,7 @@ LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
// 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 s + s * y;
+ return fputil::multiply_add(y, s, s);
}
} // namespace mathvec
>From a445d0c5b439d6fb03ba03441e60c895b41172d8 Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Thu, 23 Jul 2026 01:50:20 +0000
Subject: [PATCH 3/6] Fixed ABI error from falling back to default multiply_add
Reverted the change removing simd.h's fma, instead adding guards
to it so non-fma targets default to x * y + z, rather than the
builtin.
Replaced all calls in expf and exp10f with cpp::fma calls.
---
libc/src/__support/CPP/CMakeLists.txt | 1 -
libc/src/__support/CPP/simd.h | 39 +++++++++++--------------
libc/src/__support/mathvec/exp10f.h | 6 ++--
libc/src/__support/mathvec/expf.h | 6 ++--
libc/src/__support/mathvec/expf_utils.h | 10 +++----
5 files changed, 28 insertions(+), 34 deletions(-)
diff --git a/libc/src/__support/CPP/CMakeLists.txt b/libc/src/__support/CPP/CMakeLists.txt
index 3732595fe212f..baf8e47095fa4 100644
--- a/libc/src/__support/CPP/CMakeLists.txt
+++ b/libc/src/__support/CPP/CMakeLists.txt
@@ -238,7 +238,6 @@ add_header_library(
HDRS
simd.h
DEPENDS
- libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.cpu_features
.utility
diff --git a/libc/src/__support/CPP/simd.h b/libc/src/__support/CPP/simd.h
index ea28d39d1c84c..18d8a6f97bc14 100644
--- a/libc/src/__support/CPP/simd.h
+++ b/libc/src/__support/CPP/simd.h
@@ -19,7 +19,6 @@
#include "src/__support/CPP/tuple.h"
#include "src/__support/CPP/type_traits.h"
#include "src/__support/CPP/utility/integer_sequence.h"
-#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
@@ -210,6 +209,23 @@ 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) {
+#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)
+
+ return x * y + z;
+}
+template <typename T, size_t N>
LIBC_INLINE constexpr static simd<T, N> ceil(simd<T, N> x) {
return __builtin_elementwise_ceil(x);
}
@@ -412,27 +428,6 @@ LIBC_INLINE static cpp::simd<T, N> map(cpp::simd<T, N> v, F f) {
// TODO: where expressions, scalar overloads, ABI types.
} // namespace cpp
-
-namespace fputil {
-#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
-template <size_t N>
-LIBC_INLINE constexpr cpp::simd<float, N> multiply_add(cpp::simd<float, N> x,
- cpp::simd<float, N> y,
- cpp::simd<float, N> z) {
- return __builtin_elementwise_fma(x, y, z);
-}
-#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
-
-#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-template <size_t N>
-LIBC_INLINE constexpr cpp::simd<double, N>
-multiply_add(cpp::simd<double, N> x, cpp::simd<double, N> y,
- cpp::simd<double, N> z) {
- return __builtin_elementwise_fma(x, y, z);
-}
-#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-} // namespace fputil
-
} // namespace LIBC_NAMESPACE_DECL
#endif // LIBC_HAS_VECTOR_TYPE
diff --git a/libc/src/__support/mathvec/exp10f.h b/libc/src/__support/mathvec/exp10f.h
index 2414190e5fc07..41c625039752b 100644
--- a/libc/src/__support/mathvec/exp10f.h
+++ b/libc/src/__support/mathvec/exp10f.h
@@ -30,7 +30,7 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp10(cpp::simd<double, N> x) {
// 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 = (x * log2_10) + shift;
+ cpp::simd<double, N> z = cpp::fma(x, log2_10, shift);
cpp::simd<double, N> n = z - shift;
constexpr cpp::simd<double, N> log2_10_lo = 0x1.7f2495fb7fa6dp-53;
@@ -41,8 +41,8 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp10(cpp::simd<double, N> x) {
// 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 = n - (x * log2_10);
- r = r - (x * log2_10_lo);
+ r = cpp::fma(-x, log2_10, n);
+ r = cpp::fma(-x, log2_10_lo, r);
r = r * negative_ln2;
return eval_exp(r, z);
diff --git a/libc/src/__support/mathvec/expf.h b/libc/src/__support/mathvec/expf.h
index 3df1cf5f1108b..68c4b2139bc6e 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::fma(x, inv_ln2, shift);
cpp::simd<double, N> n = z - shift;
// ln2_hi = round(log(2), D, RN);
@@ -37,8 +37,8 @@ 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;
+ r = cpp::fma(-n, ln2_hi, r);
+ r = cpp::fma(-n, ln2_lo, r);
return eval_exp(r, z);
}
diff --git a/libc/src/__support/mathvec/expf_utils.h b/libc/src/__support/mathvec/expf_utils.h
index cd674d9ded381..336464d55455b 100644
--- a/libc/src/__support/mathvec/expf_utils.h
+++ b/libc/src/__support/mathvec/expf_utils.h
@@ -50,10 +50,10 @@ LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
// 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 = fputil::multiply_add(c1, r, c0);
- cpp::simd<double, N> p23 = fputil::multiply_add(c3, r, c2);
- cpp::simd<double, N> p04 = fputil::multiply_add(p23, r2, p01);
- cpp::simd<double, N> y = fputil::multiply_add(p04, r2, r);
+ cpp::simd<double, N> p01 = cpp::fma(c1, r, c0);
+ cpp::simd<double, N> p23 = cpp::fma(c3, r, c2);
+ cpp::simd<double, N> p04 = cpp::fma(p23, r2, p01);
+ cpp::simd<double, N> y = cpp::fma(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);
@@ -62,7 +62,7 @@ LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
// 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 fputil::multiply_add(y, s, s);
+ return cpp::fma(y, s, s);
}
} // namespace mathvec
>From 3f511cd632b4f702a66f5fca13be60c25c20fcdb Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Fri, 24 Jul 2026 20:03:49 +0000
Subject: [PATCH 4/6] Renamed simd.h fma to multiply_add to align with scalar
---
libc/src/__support/CPP/simd.h | 4 ++--
libc/src/__support/mathvec/exp10f.h | 6 +++---
libc/src/__support/mathvec/expf.h | 6 +++---
libc/src/__support/mathvec/expf_utils.h | 10 +++++-----
4 files changed, 13 insertions(+), 13 deletions(-)
diff --git a/libc/src/__support/CPP/simd.h b/libc/src/__support/CPP/simd.h
index 18d8a6f97bc14..bfa0cb5ef81bf 100644
--- a/libc/src/__support/CPP/simd.h
+++ b/libc/src/__support/CPP/simd.h
@@ -209,8 +209,8 @@ 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) {
+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>)
diff --git a/libc/src/__support/mathvec/exp10f.h b/libc/src/__support/mathvec/exp10f.h
index 41c625039752b..a3b8b5ece61dc 100644
--- a/libc/src/__support/mathvec/exp10f.h
+++ b/libc/src/__support/mathvec/exp10f.h
@@ -30,7 +30,7 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp10(cpp::simd<double, N> x) {
// 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::fma(x, log2_10, shift);
+ 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;
@@ -41,8 +41,8 @@ LIBC_INLINE static cpp::simd<double, N> inline_exp10(cpp::simd<double, N> x) {
// 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::fma(-x, log2_10, n);
- r = cpp::fma(-x, log2_10_lo, 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);
diff --git a/libc/src/__support/mathvec/expf.h b/libc/src/__support/mathvec/expf.h
index 68c4b2139bc6e..752d127b4d544 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 = cpp::fma(x, inv_ln2, shift);
+ 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,8 +37,8 @@ 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 = cpp::fma(-n, ln2_hi, r);
- r = cpp::fma(-n, ln2_lo, r);
+ r = cpp::multiply_add(-n, ln2_hi, r);
+ r = cpp::multiply_add(-n, ln2_lo, r);
return eval_exp(r, z);
}
diff --git a/libc/src/__support/mathvec/expf_utils.h b/libc/src/__support/mathvec/expf_utils.h
index 336464d55455b..8ce3d3c3b5e6e 100644
--- a/libc/src/__support/mathvec/expf_utils.h
+++ b/libc/src/__support/mathvec/expf_utils.h
@@ -50,10 +50,10 @@ LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
// 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::fma(c1, r, c0);
- cpp::simd<double, N> p23 = cpp::fma(c3, r, c2);
- cpp::simd<double, N> p04 = cpp::fma(p23, r2, p01);
- cpp::simd<double, N> y = cpp::fma(p04, r2, 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);
@@ -62,7 +62,7 @@ LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
// 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::fma(y, s, s);
+ return cpp::multiply_add(y, s, s);
}
} // namespace mathvec
>From 487f56fe5a9c66ea90e18d32b1de23deec9e0a5d Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Fri, 31 Jul 2026 18:23:41 +0000
Subject: [PATCH 5/6] Added temporary to multiply_add to prevent fusing on
non-FMA targets
---
libc/src/__support/CPP/simd.h | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/libc/src/__support/CPP/simd.h b/libc/src/__support/CPP/simd.h
index bfa0cb5ef81bf..9a91a32ffa6b0 100644
--- a/libc/src/__support/CPP/simd.h
+++ b/libc/src/__support/CPP/simd.h
@@ -223,7 +223,8 @@ LIBC_INLINE constexpr static simd<T, N> multiply_add(simd<T, N> x, simd<T, N> y,
#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
#endif // __has_builtin(__builtin_elementwise_fma)
- return x * y + z;
+ 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) {
>From fc186a41c282b45600b1bfa3b11f316eea672b35 Mon Sep 17 00:00:00 2001
From: Dylan Fleming <Dylan.Fleming at arm.com>
Date: Tue, 4 Aug 2026 12:57:58 +0000
Subject: [PATCH 6/6] Updated zero/inf bounds to be exact overflow/underflow
points
---
libc/src/__support/mathvec/exp10f.h | 4 ++--
libc/src/__support/mathvec/exp2f.h | 4 ++--
libc/src/__support/mathvec/expf.h | 4 ++--
3 files changed, 6 insertions(+), 6 deletions(-)
diff --git a/libc/src/__support/mathvec/exp10f.h b/libc/src/__support/mathvec/exp10f.h
index a3b8b5ece61dc..d22baf16b40ac 100644
--- a/libc/src/__support/mathvec/exp10f.h
+++ b/libc/src/__support/mathvec/exp10f.h
@@ -52,8 +52,8 @@ template <size_t N>
LIBC_INLINE cpp::simd<float, N> exp10f(cpp::simd<float, N> x) {
using FPBits = typename fputil::FPBits<float>;
- cpp::simd<bool, N> is_inf = x > 0x1p+7f;
- cpp::simd<bool, N> is_zero = x < -0x1p+7f;
+ 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;
diff --git a/libc/src/__support/mathvec/exp2f.h b/libc/src/__support/mathvec/exp2f.h
index 22cdc4a6ac145..848b9647dd612 100644
--- a/libc/src/__support/mathvec/exp2f.h
+++ b/libc/src/__support/mathvec/exp2f.h
@@ -48,8 +48,8 @@ template <size_t N>
LIBC_INLINE cpp::simd<float, N> exp2f(cpp::simd<float, N> x) {
using FPBits = typename fputil::FPBits<float>;
- cpp::simd<bool, N> is_inf = x > 0x1p+9f;
- cpp::simd<bool, N> is_zero = x < -0x1p+9f;
+ 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;
diff --git a/libc/src/__support/mathvec/expf.h b/libc/src/__support/mathvec/expf.h
index 752d127b4d544..95f5a5dc76874 100644
--- a/libc/src/__support/mathvec/expf.h
+++ b/libc/src/__support/mathvec/expf.h
@@ -47,8 +47,8 @@ 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;
More information about the libc-commits
mailing list