[libc-commits] [libc] [llvm] [libc][math] Improve accuracy for exp*f float-only implementations. (PR #227881)

[via libc-commits]

https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/227881

>From 81b7f32d299d91946021561e583e697fa6613af0 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Wed, 30 Sep 2026 21:39:55 +0000
Subject: [PATCH] [libc][math] Improve accuracy for exp*f float-only
 implementations.

---
 libc/src/__support/math/CMakeLists.txt        |  1 +
 libc/src/__support/math/exp2f_float_utils.h   | 25 ++++++++++++-------
 libc/src/__support/math/expf_float_eval.h     |  6 ++---
 .../llvm-project-overlay/libc/BUILD.bazel     |  1 +
 4 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index ec8b8d9eee4e5..c6a1ad1909097 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3104,6 +3104,7 @@ add_header_library(
   DEPENDS
     libc.src.__support.CPP.bit
     libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
     libc.src.__support.macros.config
     libc.src.__support.macros.optimization
 )
diff --git a/libc/src/__support/math/exp2f_float_utils.h b/libc/src/__support/math/exp2f_float_utils.h
index 1c96050bea37f..28bad5a92dd76 100644
--- a/libc/src/__support/math/exp2f_float_utils.h
+++ b/libc/src/__support/math/exp2f_float_utils.h
@@ -15,6 +15,7 @@
 #define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_FLOAT_UTILS_H
 
 #include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/PolyEval.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/optimization.h"
@@ -70,9 +71,8 @@ LIBC_INLINE float scale_exp2f(float y, int m) {
   return (y * scale) * 0x1.0p-25f;
 }
 
-// Degree-6 polynomial approximating 2^u - 1 = u * P(u) on [-0.5, 0.5] using
-// Estrin's scheme and exponent scaling.
-// Evaluates 2^x = 2^m * 2^u.
+// Degree-6 polynomial approximating 2^u - 1 = u * P(u) on [-0.5, 0.5] and
+// exponent scaling to evaluate 2^x = 2^m * 2^u.
 // Generated by Sollya with:
 // > display = hexadecimal;
 // > P = fpminimax((2^x - 1)/x, 5, [|SG...|], [-0.5, 0.5]);
@@ -80,6 +80,15 @@ LIBC_INLINE float scale_exp2f(float y, int m) {
 //   0x1.68cb...p-28
 // > dirtyinfnorm((1 + x * P - 2^x) / 2^x, [-0.5, 0.5]);
 //   0x1.5109...p-28
+//
+// Instead of using pure Estrin's scheme on P(u), we evaluate the degree-4
+// part p = c1 + c2 * u + c3 * u^2 + c4 * u^3 + c5 * u^4 with Estrin's
+// scheme, and then evaluate y = 1.0f + u * (c0 + u * p) with Horner's
+// scheme.
+// In pure Estrin's scheme, p0 = c0 + c1 * u, q0 = p0 + u^2 * p1, and
+// p = q0 + u^4 * p2 all round at the exponent of c0 (~2^-1), accumulating
+// rounding errors on non-FMA targets. Evaluating c1..c5 first keeps only a
+// single rounding at c0's scale.
 LIBC_INLINE float exp2f_eval(float u, int m) {
   constexpr float COEFFS[6] = {0x1.62e430p-1f,  0x1.ebfbdep-3f,
                                0x1.c6af9ep-5f,  0x1.3b2c58p-7f,
@@ -88,14 +97,12 @@ LIBC_INLINE float exp2f_eval(float u, int m) {
   float u2 = u * u;
   float u4 = u2 * u2;
 
-  float p0 = fputil::multiply_add(u, COEFFS[1], COEFFS[0]);
-  float p1 = fputil::multiply_add(u, COEFFS[3], COEFFS[2]);
-  float p2 = fputil::multiply_add(u, COEFFS[5], COEFFS[4]);
-
+  float p0 = fputil::multiply_add(u, COEFFS[2], COEFFS[1]);
+  float p1 = fputil::multiply_add(u, COEFFS[4], COEFFS[3]);
   float q0 = fputil::multiply_add(u2, p1, p0);
-  float p = fputil::multiply_add(u4, p2, q0);
+  float p = fputil::multiply_add(u4, COEFFS[5], q0);
 
-  float y = fputil::multiply_add(u, p, 1.0f);
+  float y = fputil::polyeval(u, 1.0f, COEFFS[0], p);
 
   return scale_exp2f(y, m);
 }
diff --git a/libc/src/__support/math/expf_float_eval.h b/libc/src/__support/math/expf_float_eval.h
index e479cb58d2fea..ab8ad649c539c 100644
--- a/libc/src/__support/math/expf_float_eval.h
+++ b/libc/src/__support/math/expf_float_eval.h
@@ -100,11 +100,11 @@ LIBC_INLINE float expf(float x) {
   // Constants generated by Sollya with:
   // > display = hexadecimal;
   // > LOG2_E = round(log2(exp(1)), SG, RN);
-  // > LOG_2_HI = round(log(2), 12, RN);
+  // > LOG_2_HI = round(log(2), 16, RN);
   // > LOG_2_LO = round(log(2) - LOG_2_HI, SG, RN);
   constexpr float LOG2_E = 0x1.715476p+0f;
-  constexpr float LOG_2_HI = 0x1.62ep-1f;
-  constexpr float LOG_2_LO = 0x1.0bfbe8p-15f;
+  constexpr float LOG_2_HI = 0x1.62e4p-1f;
+  constexpr float LOG_2_LO = 0x1.7f7d1cp-20f;
 
   float kf = fputil::nearest_integer(x * LOG2_E);
   int k = static_cast<int>(kf);
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index b6ec7c595fced..f2bd619962b39 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -8701,6 +8701,7 @@ libc_support_library(
         ":__support_common",
         ":__support_cpp_bit",
         ":__support_fputil_multiply_add",
+        ":__support_fputil_polyeval",
         ":__support_macros_config",
         ":__support_macros_optimization",
     ],