[libc] [llvm] [libc][math] Refactor exp10f implementation to header-only in src/__support/math folder. (PR #148405)

Muhammad Bassiouni via llvm-commits llvm-commits at lists.llvm.org
Wed Jul 16 12:57:27 PDT 2025


https://github.com/bassiounix updated https://github.com/llvm/llvm-project/pull/148405

>From 6282f05226fa87566f042a71479fe194c6d95065 Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Fri, 11 Jul 2025 03:12:38 +0300
Subject: [PATCH 1/6] [libc][math] Refactor exp implementation to header-only
 in src/__support/math folder.

---
 libc/shared/math.h                       | 1 +
 libc/src/math/generic/common_constants.h | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/libc/shared/math.h b/libc/shared/math.h
index b37aa46820523..ac422f36bb2bb 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -21,5 +21,6 @@
 #include "math/ldexpf.h"
 #include "math/ldexpf128.h"
 #include "math/ldexpf16.h"
+#include "math/exp.h"
 
 #endif // LLVM_LIBC_SHARED_MATH_H
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index 291816a7889ad..a1d6de3c87240 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -13,6 +13,8 @@
 #include "src/__support/macros/config.h"
 #include "src/__support/math/exp_constants.h"
 #include "src/__support/number_pair.h"
+#include "src/__support/math/exp_constants.h"
+
 
 namespace LIBC_NAMESPACE_DECL {
 

>From a4dc3df8357d91aec8d03368c5814cdc1053144d Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Fri, 11 Jul 2025 03:19:04 +0300
Subject: [PATCH 2/6] fix style

---
 libc/shared/math.h                       | 1 -
 libc/src/math/generic/common_constants.h | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/libc/shared/math.h b/libc/shared/math.h
index ac422f36bb2bb..b37aa46820523 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -21,6 +21,5 @@
 #include "math/ldexpf.h"
 #include "math/ldexpf128.h"
 #include "math/ldexpf16.h"
-#include "math/exp.h"
 
 #endif // LLVM_LIBC_SHARED_MATH_H
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index a1d6de3c87240..80eee06b2d4d9 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -14,7 +14,7 @@
 #include "src/__support/math/exp_constants.h"
 #include "src/__support/number_pair.h"
 #include "src/__support/math/exp_constants.h"
-
+#include "src/__support/number_pair.h"
 
 namespace LIBC_NAMESPACE_DECL {
 

>From 72de600c2b207bb88f310605b08b8b521a7d3c83 Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Sun, 13 Jul 2025 00:40:16 +0300
Subject: [PATCH 3/6] [libc][math] Refactor exp10 implementation to header-only
 in src/__support/math folder.

---
 libc/src/math/generic/common_constants.h | 2 --
 1 file changed, 2 deletions(-)

diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index 80eee06b2d4d9..291816a7889ad 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -13,8 +13,6 @@
 #include "src/__support/macros/config.h"
 #include "src/__support/math/exp_constants.h"
 #include "src/__support/number_pair.h"
-#include "src/__support/math/exp_constants.h"
-#include "src/__support/number_pair.h"
 
 namespace LIBC_NAMESPACE_DECL {
 

>From 46f7684bfd488fb6d0c1dbed68e05ad6b553e7da Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Fri, 11 Jul 2025 03:12:38 +0300
Subject: [PATCH 4/6] [libc][math] Refactor exp implementation to header-only
 in src/__support/math folder.

---
 libc/shared/math.h                       | 1 +
 libc/src/math/generic/common_constants.h | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/libc/shared/math.h b/libc/shared/math.h
index b37aa46820523..ac422f36bb2bb 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -21,5 +21,6 @@
 #include "math/ldexpf.h"
 #include "math/ldexpf128.h"
 #include "math/ldexpf16.h"
+#include "math/exp.h"
 
 #endif // LLVM_LIBC_SHARED_MATH_H
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index 291816a7889ad..a1d6de3c87240 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -13,6 +13,8 @@
 #include "src/__support/macros/config.h"
 #include "src/__support/math/exp_constants.h"
 #include "src/__support/number_pair.h"
+#include "src/__support/math/exp_constants.h"
+
 
 namespace LIBC_NAMESPACE_DECL {
 

>From dc6914c3d122b5cdda42870f4ac7a32540fa40de Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Fri, 11 Jul 2025 03:19:04 +0300
Subject: [PATCH 5/6] fix style

---
 libc/shared/math.h                       | 1 -
 libc/src/math/generic/common_constants.h | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/libc/shared/math.h b/libc/shared/math.h
index ac422f36bb2bb..b37aa46820523 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -21,6 +21,5 @@
 #include "math/ldexpf.h"
 #include "math/ldexpf128.h"
 #include "math/ldexpf16.h"
-#include "math/exp.h"
 
 #endif // LLVM_LIBC_SHARED_MATH_H
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index a1d6de3c87240..80eee06b2d4d9 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -14,7 +14,7 @@
 #include "src/__support/math/exp_constants.h"
 #include "src/__support/number_pair.h"
 #include "src/__support/math/exp_constants.h"
-
+#include "src/__support/number_pair.h"
 
 namespace LIBC_NAMESPACE_DECL {
 

>From 392da0bfb0088cb609f4f078a0ff80add0f9e4c1 Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Sun, 13 Jul 2025 04:42:34 +0300
Subject: [PATCH 6/6] [libc][math] Refactor exp10f implementation to
 header-only in src/__support/math folder.

---
 libc/shared/math.h                            |   1 +
 libc/shared/math/exp10f.h                     |  24 +++
 libc/src/__support/math/CMakeLists.txt        |  28 +++
 .../exp10f_impl.h => __support/math/exp10f.h} |  17 +-
 libc/src/__support/math/exp10f_utils.h        | 171 ++++++++++++++++++
 libc/src/math/generic/CMakeLists.txt          |  33 +---
 libc/src/math/generic/atanhf.cpp              |   1 +
 libc/src/math/generic/coshf.cpp               |   2 +-
 libc/src/math/generic/exp10f.cpp              |   7 +-
 libc/src/math/generic/explogxf.h              | 156 +---------------
 libc/src/math/generic/powf.cpp                |   7 +-
 libc/src/math/generic/sinhf.cpp               |   1 +
 .../llvm-project-overlay/libc/BUILD.bazel     |  53 ++++--
 13 files changed, 281 insertions(+), 220 deletions(-)
 create mode 100644 libc/shared/math/exp10f.h
 rename libc/src/{math/generic/exp10f_impl.h => __support/math/exp10f.h} (91%)
 create mode 100644 libc/src/__support/math/exp10f_utils.h

diff --git a/libc/shared/math.h b/libc/shared/math.h
index b37aa46820523..2ae7c1d58ae10 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -13,6 +13,7 @@
 
 #include "math/exp.h"
 #include "math/exp10.h"
+#include "math/exp10f.h"
 #include "math/expf.h"
 #include "math/expf16.h"
 #include "math/frexpf.h"
diff --git a/libc/shared/math/exp10f.h b/libc/shared/math/exp10f.h
new file mode 100644
index 0000000000000..556e78ab3b7a7
--- /dev/null
+++ b/libc/shared/math/exp10f.h
@@ -0,0 +1,24 @@
+//===-- Shared exp10f function -----------------------------------*- C++
+//-*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_EXP10F_H
+#define LLVM_LIBC_SHARED_MATH_EXP10F_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/exp10f.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::exp10f;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_EXP10F_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 0bfc996c44fc8..ad36679409f89 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -170,3 +170,31 @@ add_header_library(
     libc.src.__support.integer_literals
     libc.src.__support.macros.optimization
 )
+
+add_header_library(
+  exp10f_utils
+  HDRS
+    exp10f_utils.h
+  DEPENDS
+    libc.src.__support.FPUtil.basic_operations
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.nearest_integer
+    libc.src.__support.FPUtil.polyeval
+    libc.src.__support.common
+    libc.src.__support.math.exp_utils
+)
+
+add_header_library(
+  exp10f
+  HDRS
+    exp10f.h
+  DEPENDS
+    .exp10f_utils
+    libc.src.__support.macros.config
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.optimization
+)
diff --git a/libc/src/math/generic/exp10f_impl.h b/libc/src/__support/math/exp10f.h
similarity index 91%
rename from libc/src/math/generic/exp10f_impl.h
rename to libc/src/__support/math/exp10f.h
index 975fd01a0a25c..807b4f0d6c109 100644
--- a/libc/src/math/generic/exp10f_impl.h
+++ b/libc/src/__support/math/exp10f.h
@@ -1,4 +1,4 @@
-//===-- Single-precision 10^x function ------------------------------------===//
+//===-- Implementation header for exp10f ------------------------*- C++ -*-===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
@@ -6,22 +6,21 @@
 //
 //===----------------------------------------------------------------------===//
 
-#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H
-#define LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H
 
-#include "explogxf.h"
+#include "exp10f_utils.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/rounding_mode.h"
-#include "src/__support/common.h"
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
 
 namespace LIBC_NAMESPACE_DECL {
-namespace generic {
+namespace math {
 
-LIBC_INLINE float exp10f(float x) {
+static constexpr float exp10f(float x) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
@@ -132,7 +131,7 @@ LIBC_INLINE float exp10f(float x) {
   return static_cast<float>(multiply_add(p, lo2 * rr.mh, c0 * rr.mh));
 }
 
-} // namespace generic
+} // namespace math
 } // namespace LIBC_NAMESPACE_DECL
 
-#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H
diff --git a/libc/src/__support/math/exp10f_utils.h b/libc/src/__support/math/exp10f_utils.h
new file mode 100644
index 0000000000000..c491872cc725c
--- /dev/null
+++ b/libc/src/__support/math/exp10f_utils.h
@@ -0,0 +1,171 @@
+//===-- Common utils for exp10f ---------------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+struct ExpBase {
+  // Base = e
+  static constexpr int MID_BITS = 5;
+  static constexpr int MID_MASK = (1 << MID_BITS) - 1;
+  // log2(e) * 2^5
+  static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS);
+  // High and low parts of -log(2) * 2^(-5)
+  static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS);
+  static constexpr double M_LOGB_2_LO =
+      -0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS);
+  // Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya
+  // with:
+  // > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN));
+  static constexpr int64_t EXP_2_MID[1 << MID_BITS] = {
+      0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f,
+      0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa,
+      0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715,
+      0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d,
+      0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429,
+      0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74,
+      0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db,
+      0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d,
+      0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c,
+      0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f,
+      0x3ffea4afa2a490da, 0x3fff50765b6e4540,
+  };
+
+  // Approximating e^dx with degree-5 minimax polynomial generated by Sollya:
+  // > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]);
+  // Then:
+  //   e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5.
+  static constexpr double COEFFS[4] = {
+      0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5,
+      0x1.11112a0e34bdbp-7};
+
+  LIBC_INLINE static double powb_lo(double dx) {
+    using fputil::multiply_add;
+    double dx2 = dx * dx;
+    double c0 = 1.0 + dx;
+    // c1 = COEFFS[0] + COEFFS[1] * dx
+    double c1 = multiply_add(dx, ExpBase::COEFFS[1], ExpBase::COEFFS[0]);
+    // c2 = COEFFS[2] + COEFFS[3] * dx
+    double c2 = multiply_add(dx, ExpBase::COEFFS[3], ExpBase::COEFFS[2]);
+    // r = c4 + c5 * dx^4
+    //   = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7
+    return fputil::polyeval(dx2, c0, c1, c2);
+  }
+};
+
+struct Exp10Base : public ExpBase {
+  // log2(10) * 2^5
+  static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS);
+  // High and low parts of -log10(2) * 2^(-5).
+  // Notice that since |x * log2(10)| < 150:
+  //   |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13
+  // So when the FMA instructions are not available, in order for the product
+  //   k * M_LOGB_2_HI
+  // to be exact, we only store the high part of log10(2) up to 38 bits
+  // (= 53 - 15) of precision.
+  // It is generated by Sollya with:
+  // > round(log10(2), 44, RN);
+  static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS);
+  // > round(log10(2) - 0x1.34413509f8p-2, D, RN);
+  static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS);
+
+  // Approximating 10^dx with degree-5 minimax polynomial generated by Sollya:
+  // > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]);
+  // Then:
+  //   10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
+  static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1,
+                                       0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0,
+                                       0x1.1429e74a98f43p-1};
+
+  static double powb_lo(double dx) {
+    using fputil::multiply_add;
+    double dx2 = dx * dx;
+    // c0 = 1 + COEFFS[0] * dx
+    double c0 = multiply_add(dx, Exp10Base::COEFFS[0], 1.0);
+    // c1 = COEFFS[1] + COEFFS[2] * dx
+    double c1 = multiply_add(dx, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]);
+    // c2 = COEFFS[3] + COEFFS[4] * dx
+    double c2 = multiply_add(dx, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]);
+    // r = c0 + dx^2 * (c1 + c2 * dx^2)
+    //   = c0 + c1 * dx^2 + c2 * dx^4
+    //   = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
+    return fputil::polyeval(dx2, c0, c1, c2);
+  }
+};
+
+constexpr int LOG_P1_BITS = 6;
+constexpr int LOG_P1_SIZE = 1 << LOG_P1_BITS;
+
+// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
+extern const double LOG_P1_LOG2[LOG_P1_SIZE];
+
+// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
+extern const double LOG_P1_1_OVER[LOG_P1_SIZE];
+
+// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
+// K_LOG2_ODD starts from x^3
+extern const double K_LOG2_ODD[4];
+extern const double K_LOG2_EVEN[4];
+
+// Output of range reduction for exp_b: (2^(mid + hi), lo)
+// where:
+//   b^x = 2^(mid + hi) * b^lo
+struct exp_b_reduc_t {
+  double mh; // 2^(mid + hi)
+  double lo;
+};
+
+// The function correctly calculates b^x value with at least float precision
+// in a limited range.
+// Range reduction:
+//   b^x = 2^(hi + mid) * b^lo
+// where:
+//   x = (hi + mid) * log_b(2) + lo
+//   hi is an integer,
+//   0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer
+//   -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1)
+// Base class needs to provide the following constants:
+//   - MID_BITS    : number of bits after decimal points used for mid
+//   - MID_MASK    : 2^MID_BITS - 1, mask to extract mid bits
+//   - LOG2_B      : log2(b) * 2^MID_BITS for scaling
+//   - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS)
+//   - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS)
+//   - EXP_2_MID   : look up table for bit fields of 2^mid
+// Return:
+//   { 2^(hi + mid), lo }
+template <class Base>
+LIBC_INLINE static constexpr exp_b_reduc_t exp_b_range_reduc(float x) {
+  double xd = static_cast<double>(x);
+  // kd = round((hi + mid) * log2(b) * 2^MID_BITS)
+  double kd = fputil::nearest_integer(Base::LOG2_B * xd);
+  // k = round((hi + mid) * log2(b) * 2^MID_BITS)
+  int k = static_cast<int>(kd);
+  // hi = floor(kd * 2^(-MID_BITS))
+  // exp_hi = shift hi to the exponent field of double precision.
+  uint64_t exp_hi = static_cast<uint64_t>(k >> Base::MID_BITS)
+                    << fputil::FPBits<double>::FRACTION_LEN;
+  // mh = 2^hi * 2^mid
+  // mh_bits = bit field of mh
+  uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi;
+  double mh = fputil::FPBits<double>(mh_bits).get_val();
+  // dx = lo = x - (hi + mid) * log(2)
+  double dx = fputil::multiply_add(
+      kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd));
+  return {mh, dx};
+}
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 352c2ad4ab22a..ebdd3d7235141 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -358,7 +358,6 @@ add_entrypoint_object(
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fma
-    libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.polyeval
     libc.src.__support.macros.optimization
 )
@@ -448,7 +447,6 @@ add_entrypoint_object(
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.fma
-    libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.polyeval
     libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.macros.optimization
@@ -1461,21 +1459,6 @@ add_entrypoint_object(
     libc.src.errno.errno
 )
 
-add_header_library(
-  exp10f_impl
-  HDRS
-    exp10f_impl.h
-  DEPENDS
-    .explogxf
-    libc.src.__support.FPUtil.fenv_impl
-    libc.src.__support.FPUtil.fp_bits
-    libc.src.__support.FPUtil.multiply_add
-    libc.src.__support.FPUtil.rounding_mode
-    libc.src.__support.macros.optimization
-    libc.src.__support.common
-    libc.src.errno.errno
-)
-
 add_entrypoint_object(
   exp10f
   SRCS
@@ -1483,7 +1466,8 @@ add_entrypoint_object(
   HDRS
     ../exp10f.h
   DEPENDS
-    .exp10f_impl
+    libc.src.__support.math.exp10f
+    libc.src.errno.errno
 )
 
 add_entrypoint_object(
@@ -1620,17 +1604,15 @@ add_entrypoint_object(
     ../powf.h
   DEPENDS
     .common_constants
-    .exp10f_impl
     .exp2f_impl
     .explogxf
+    libc.src.__support.math.exp10f
     libc.src.__support.CPP.bit
-    libc.src.__support.CPP.optional
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
-    libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.FPUtil.sqrt
     libc.src.__support.FPUtil.triple_double
     libc.src.__support.macros.optimization
@@ -3792,13 +3774,9 @@ add_object_library(
     explogxf.cpp
   DEPENDS
     .common_constants
-    libc.src.__support.FPUtil.basic_operations
-    libc.src.__support.FPUtil.fenv_impl
-    libc.src.__support.FPUtil.multiply_add
-    libc.src.__support.FPUtil.nearest_integer
-    libc.src.__support.FPUtil.polyeval
-    libc.src.__support.common
     libc.src.__support.math.exp_utils
+    libc.src.__support.math.exp10f_utils
+    libc.src.__support.macros.properties.cpu_features
     libc.src.errno.errno
 )
 
@@ -3981,6 +3959,7 @@ add_entrypoint_object(
   DEPENDS
     .explogxf
     libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.macros.optimization
 )
 
diff --git a/libc/src/math/generic/atanhf.cpp b/libc/src/math/generic/atanhf.cpp
index 2149314d2f676..f6fde766ef785 100644
--- a/libc/src/math/generic/atanhf.cpp
+++ b/libc/src/math/generic/atanhf.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/atanhf.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
diff --git a/libc/src/math/generic/coshf.cpp b/libc/src/math/generic/coshf.cpp
index c869f7d9dec5f..9f87564d524a6 100644
--- a/libc/src/math/generic/coshf.cpp
+++ b/libc/src/math/generic/coshf.cpp
@@ -7,8 +7,8 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/coshf.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/rounding_mode.h"
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
diff --git a/libc/src/math/generic/exp10f.cpp b/libc/src/math/generic/exp10f.cpp
index 5284c380f52ec..b2d4f097bc7ce 100644
--- a/libc/src/math/generic/exp10f.cpp
+++ b/libc/src/math/generic/exp10f.cpp
@@ -7,12 +7,11 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/exp10f.h"
-#include "src/__support/common.h"
-#include "src/__support/macros/config.h"
-#include "src/math/generic/exp10f_impl.h"
+
+#include "src/__support/math/exp10f.h"
 
 namespace LIBC_NAMESPACE_DECL {
 
-LLVM_LIBC_FUNCTION(float, exp10f, (float x)) { return generic::exp10f(x); }
+LLVM_LIBC_FUNCTION(float, exp10f, (float x)) { return math::exp10f(x); }
 
 } // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h
index 5ae1457ca780e..d9f39a0d6080f 100644
--- a/libc/src/math/generic/explogxf.h
+++ b/libc/src/math/generic/explogxf.h
@@ -10,166 +10,14 @@
 #define LLVM_LIBC_SRC_MATH_GENERIC_EXPLOGXF_H
 
 #include "common_constants.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/FPUtil/nearest_integer.h"
+
 #include "src/__support/common.h"
-#include "src/__support/macros/config.h"
 #include "src/__support/macros/properties/cpu_features.h"
-
+#include "src/__support/math/exp10f_utils.h"
 #include "src/__support/math/exp_utils.h"
 
 namespace LIBC_NAMESPACE_DECL {
 
-struct ExpBase {
-  // Base = e
-  static constexpr int MID_BITS = 5;
-  static constexpr int MID_MASK = (1 << MID_BITS) - 1;
-  // log2(e) * 2^5
-  static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS);
-  // High and low parts of -log(2) * 2^(-5)
-  static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS);
-  static constexpr double M_LOGB_2_LO =
-      -0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS);
-  // Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya
-  // with:
-  // > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN));
-  static constexpr int64_t EXP_2_MID[1 << MID_BITS] = {
-      0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f,
-      0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa,
-      0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715,
-      0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d,
-      0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429,
-      0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74,
-      0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db,
-      0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d,
-      0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c,
-      0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f,
-      0x3ffea4afa2a490da, 0x3fff50765b6e4540,
-  };
-
-  // Approximating e^dx with degree-5 minimax polynomial generated by Sollya:
-  // > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]);
-  // Then:
-  //   e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5.
-  static constexpr double COEFFS[4] = {
-      0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5,
-      0x1.11112a0e34bdbp-7};
-
-  LIBC_INLINE static double powb_lo(double dx) {
-    using fputil::multiply_add;
-    double dx2 = dx * dx;
-    double c0 = 1.0 + dx;
-    // c1 = COEFFS[0] + COEFFS[1] * dx
-    double c1 = multiply_add(dx, ExpBase::COEFFS[1], ExpBase::COEFFS[0]);
-    // c2 = COEFFS[2] + COEFFS[3] * dx
-    double c2 = multiply_add(dx, ExpBase::COEFFS[3], ExpBase::COEFFS[2]);
-    // r = c4 + c5 * dx^4
-    //   = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7
-    return fputil::polyeval(dx2, c0, c1, c2);
-  }
-};
-
-struct Exp10Base : public ExpBase {
-  // log2(10) * 2^5
-  static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS);
-  // High and low parts of -log10(2) * 2^(-5).
-  // Notice that since |x * log2(10)| < 150:
-  //   |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13
-  // So when the FMA instructions are not available, in order for the product
-  //   k * M_LOGB_2_HI
-  // to be exact, we only store the high part of log10(2) up to 38 bits
-  // (= 53 - 15) of precision.
-  // It is generated by Sollya with:
-  // > round(log10(2), 44, RN);
-  static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS);
-  // > round(log10(2) - 0x1.34413509f8p-2, D, RN);
-  static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS);
-
-  // Approximating 10^dx with degree-5 minimax polynomial generated by Sollya:
-  // > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]);
-  // Then:
-  //   10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
-  static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1,
-                                       0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0,
-                                       0x1.1429e74a98f43p-1};
-
-  static double powb_lo(double dx) {
-    using fputil::multiply_add;
-    double dx2 = dx * dx;
-    // c0 = 1 + COEFFS[0] * dx
-    double c0 = multiply_add(dx, Exp10Base::COEFFS[0], 1.0);
-    // c1 = COEFFS[1] + COEFFS[2] * dx
-    double c1 = multiply_add(dx, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]);
-    // c2 = COEFFS[3] + COEFFS[4] * dx
-    double c2 = multiply_add(dx, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]);
-    // r = c0 + dx^2 * (c1 + c2 * dx^2)
-    //   = c0 + c1 * dx^2 + c2 * dx^4
-    //   = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
-    return fputil::polyeval(dx2, c0, c1, c2);
-  }
-};
-
-constexpr int LOG_P1_BITS = 6;
-constexpr int LOG_P1_SIZE = 1 << LOG_P1_BITS;
-
-// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
-extern const double LOG_P1_LOG2[LOG_P1_SIZE];
-
-// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
-extern const double LOG_P1_1_OVER[LOG_P1_SIZE];
-
-// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
-// K_LOG2_ODD starts from x^3
-extern const double K_LOG2_ODD[4];
-extern const double K_LOG2_EVEN[4];
-
-// Output of range reduction for exp_b: (2^(mid + hi), lo)
-// where:
-//   b^x = 2^(mid + hi) * b^lo
-struct exp_b_reduc_t {
-  double mh; // 2^(mid + hi)
-  double lo;
-};
-
-// The function correctly calculates b^x value with at least float precision
-// in a limited range.
-// Range reduction:
-//   b^x = 2^(hi + mid) * b^lo
-// where:
-//   x = (hi + mid) * log_b(2) + lo
-//   hi is an integer,
-//   0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer
-//   -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1)
-// Base class needs to provide the following constants:
-//   - MID_BITS    : number of bits after decimal points used for mid
-//   - MID_MASK    : 2^MID_BITS - 1, mask to extract mid bits
-//   - LOG2_B      : log2(b) * 2^MID_BITS for scaling
-//   - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS)
-//   - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS)
-//   - EXP_2_MID   : look up table for bit fields of 2^mid
-// Return:
-//   { 2^(hi + mid), lo }
-template <class Base> LIBC_INLINE exp_b_reduc_t exp_b_range_reduc(float x) {
-  double xd = static_cast<double>(x);
-  // kd = round((hi + mid) * log2(b) * 2^MID_BITS)
-  double kd = fputil::nearest_integer(Base::LOG2_B * xd);
-  // k = round((hi + mid) * log2(b) * 2^MID_BITS)
-  int k = static_cast<int>(kd);
-  // hi = floor(kd * 2^(-MID_BITS))
-  // exp_hi = shift hi to the exponent field of double precision.
-  uint64_t exp_hi = static_cast<uint64_t>(k >> Base::MID_BITS)
-                    << fputil::FPBits<double>::FRACTION_LEN;
-  // mh = 2^hi * 2^mid
-  // mh_bits = bit field of mh
-  uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi;
-  double mh = fputil::FPBits<double>(mh_bits).get_val();
-  // dx = lo = x - (hi + mid) * log(2)
-  double dx = fputil::multiply_add(
-      kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd));
-  return {mh, dx};
-}
-
 // The function correctly calculates sinh(x) and cosh(x) by calculating exp(x)
 // and exp(-x) simultaneously.
 // To compute e^x, we perform the following range
diff --git a/libc/src/math/generic/powf.cpp b/libc/src/math/generic/powf.cpp
index dfdfd5d6d5760..a45ef511c9bad 100644
--- a/libc/src/math/generic/powf.cpp
+++ b/libc/src/math/generic/powf.cpp
@@ -9,20 +9,17 @@
 #include "src/math/powf.h"
 #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
 #include "src/__support/CPP/bit.h"
-#include "src/__support/CPP/optional.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
 #include "src/__support/FPUtil/double_double.h"
-#include "src/__support/FPUtil/except_value_utils.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/nearest_integer.h"
-#include "src/__support/FPUtil/rounding_mode.h"
 #include "src/__support/FPUtil/sqrt.h" // Speedup for powf(x, 1/2) = sqrtf(x)
 #include "src/__support/common.h"
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/math/exp10f.h" // Speedup for powf(10, y) = exp10f(y)
 
-#include "exp10f_impl.h" // Speedup for powf(10, y) = exp10f(y)
 #include "exp2f_impl.h"  // Speedup for powf(2, y) = exp2f(y)
 
 namespace LIBC_NAMESPACE_DECL {
@@ -781,7 +778,7 @@ LLVM_LIBC_FUNCTION(float, powf, (float x, float y)) {
       return generic::exp2f(y);
     case 0x4120'0000: // x = 10.0f
       // pow(10, y) = exp10(y)
-      return generic::exp10f(y);
+      return math::exp10f(y);
 #endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
     }
 
diff --git a/libc/src/math/generic/sinhf.cpp b/libc/src/math/generic/sinhf.cpp
index d6158fd302536..63111f84de141 100644
--- a/libc/src/math/generic/sinhf.cpp
+++ b/libc/src/math/generic/sinhf.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/sinhf.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/rounding_mode.h"
 #include "src/__support/macros/config.h"
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 8b60ca13562f6..4463d28864752 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -1999,13 +1999,13 @@ libc_support_library(
     srcs = ["src/math/generic/explogxf.cpp"],
     hdrs = ["src/math/generic/explogxf.h"],
     deps = [
-        ":__support_common",
         ":__support_fputil_fenv_impl",
         ":__support_fputil_fma",
         ":__support_fputil_multiply_add",
         ":__support_fputil_nearest_integer",
-        ":__support_fputil_polyeval",
         ":__support_math_exp_utils",
+        ":__support_math_exp10f_utils",
+        ":__support_macros_properties_cpu_features",
         ":common_constants",
     ],
 )
@@ -2050,19 +2050,6 @@ libc_support_library(
     ],
 )
 
-libc_support_library(
-    name = "exp10f_impl",
-    hdrs = ["src/math/generic/exp10f_impl.h"],
-    deps = [
-        ":__support_fputil_fma",
-        ":__support_fputil_multiply_add",
-        ":__support_fputil_rounding_mode",
-        ":__support_macros_optimization",
-        ":common_constants",
-        ":explogxf",
-    ],
-)
-
 libc_support_library(
     name = "exp2f_impl",
     hdrs = ["src/math/generic/exp2f_impl.h"],
@@ -2263,6 +2250,33 @@ libc_support_library(
     ],
 )
 
+libc_support_library(
+    name = "__support_math_exp10f_utils",
+    hdrs = ["src/__support/math/exp10f_utils.h"],
+    deps = [
+        ":__support_fputil_basic_operations",
+        ":__support_fputil_fenv_impl",
+        ":__support_fputil_multiply_add",
+        ":__support_fputil_nearest_integer",
+        ":__support_fputil_polyeval",
+        ":__support_common",
+        ":__support_math_exp_utils",
+    ],
+)
+
+libc_support_library(
+    name = "__support_math_exp10f",
+    hdrs = ["src/__support/math/exp10f.h"],
+    deps = [
+        ":__support_math_exp10f_utils",
+        ":__support_fputil_fenv_impl",
+        ":__support_fputil_fp_bits",
+        ":__support_fputil_multiply_add",
+        ":__support_fputil_rounding_mode",
+        ":__support_macros_optimization",
+    ],
+)
+
 ############################### complex targets ################################
 
 libc_function(
@@ -2726,10 +2740,10 @@ libc_math_function(
     name = "cosf",
     additional_deps = [
         ":__support_fputil_fma",
-        ":__support_fputil_multiply_add",
         ":__support_macros_optimization",
         ":__support_macros_properties_cpu_features",
         ":sincosf_utils",
+        ":errno",
     ],
 )
 
@@ -2875,7 +2889,8 @@ libc_math_function(
 libc_math_function(
     name = "exp10f",
     additional_deps = [
-        ":exp10f_impl",
+        ":__support_math_exp10f",
+        ":errno",
     ],
 )
 
@@ -3724,14 +3739,13 @@ libc_math_function(
         ":__support_fputil_multiply_add",
         ":__support_fputil_nearest_integer",
         ":__support_fputil_polyeval",
-        ":__support_fputil_rounding_mode",
         ":__support_fputil_sqrt",
         ":__support_fputil_triple_double",
         ":__support_macros_optimization",
+        ":__support_math_exp10f",
         ":common_constants",
         ":explogxf",
         ":exp2f_impl",
-        ":exp10f_impl",
     ],
 )
 
@@ -3840,7 +3854,6 @@ libc_math_function(
     name = "sinf",
     additional_deps = [
         ":__support_fputil_fma",
-        ":__support_fputil_multiply_add",
         ":__support_fputil_polyeval",
         ":__support_fputil_rounding_mode",
         ":__support_macros_optimization",



More information about the llvm-commits mailing list