[libc] [llvm] [libc][math] Refactor exp2 implementation to header-only in src/__support/math folder. (PR #161297)
Muhammad Bassiouni via llvm-commits
llvm-commits at lists.llvm.org
Wed Oct 1 01:24:18 PDT 2025
https://github.com/bassiounix updated https://github.com/llvm/llvm-project/pull/161297
>From a4cf62e28727ac1275787a38ab0ed6d2c70b8711 Mon Sep 17 00:00:00 2001
From: bassiounix <muhammad.m.bassiouni at gmail.com>
Date: Tue, 30 Sep 2025 02:28:48 +0300
Subject: [PATCH] [libc][math] Refactor exp2 implementation to header-only in
src/__support/math folder.
---
libc/shared/math.h | 1 +
libc/shared/math/exp2.h | 23 +
libc/src/__support/math/CMakeLists.txt | 31 ++
.../math/common_constants.h} | 47 +-
libc/src/__support/math/exp2.h | 425 ++++++++++++++++++
libc/src/math/generic/CMakeLists.txt | 56 +--
libc/src/math/generic/common_constants.h | 73 ---
libc/src/math/generic/exp2.cpp | 398 +---------------
libc/src/math/generic/expm1.cpp | 5 +-
libc/src/math/generic/expm1f.cpp | 3 +-
libc/src/math/generic/log.cpp | 4 +-
libc/src/math/generic/log10.cpp | 5 +-
libc/src/math/generic/log10f.cpp | 3 +-
libc/src/math/generic/log1p.cpp | 4 +-
libc/src/math/generic/log1pf.cpp | 4 +-
libc/src/math/generic/log2.cpp | 5 +-
libc/src/math/generic/log2f.cpp | 5 +-
libc/src/math/generic/log_range_reduction.h | 3 +-
libc/src/math/generic/logf.cpp | 3 +-
libc/src/math/generic/pow.cpp | 5 +-
libc/src/math/generic/powf.cpp | 6 +-
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 1 +
.../llvm-project-overlay/libc/BUILD.bazel | 85 ++--
24 files changed, 620 insertions(+), 576 deletions(-)
create mode 100644 libc/shared/math/exp2.h
rename libc/src/{math/generic/common_constants.cpp => __support/math/common_constants.h} (95%)
create mode 100644 libc/src/__support/math/exp2.h
delete mode 100644 libc/src/math/generic/common_constants.h
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 4b2a0d8c712ad..924d0cb4f3bfa 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -47,6 +47,7 @@
#include "math/exp10f16.h"
#include "math/exp10m1f.h"
#include "math/exp10m1f16.h"
+#include "math/exp2.h"
#include "math/expf.h"
#include "math/expf16.h"
#include "math/frexpf.h"
diff --git a/libc/shared/math/exp2.h b/libc/shared/math/exp2.h
new file mode 100644
index 0000000000000..6f1e143376a8a
--- /dev/null
+++ b/libc/shared/math/exp2.h
@@ -0,0 +1,23 @@
+//===-- Shared exp2 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_EXP2_H
+#define LLVM_LIBC_SHARED_MATH_EXP2_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/exp2.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::exp2;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_EXP2_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 98f9bb42f91f4..4130fdf02078b 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -373,6 +373,15 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ common_constants
+ HDRS
+ common_constants.h
+ DEPENDS
+ libc.src.__support.macros.config
+ libc.src.__support.number_pair
+)
+
add_header_library(
cos
HDRS
@@ -704,6 +713,28 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ exp2
+ HDRS
+ exp2.h
+ DEPENDS
+ .common_constants
+ .exp_utils
+ libc.src.__support.CPP.bit
+ libc.src.__support.CPP.optional
+ libc.src.__support.FPUtil.dyadic_float
+ 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.triple_double
+ libc.src.__support.integer_literals
+ libc.src.__support.macros.optimization
+ libc.src.errno.errno
+)
+
add_header_library(
exp10
HDRS
diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/__support/math/common_constants.h
similarity index 95%
rename from libc/src/math/generic/common_constants.cpp
rename to libc/src/__support/math/common_constants.h
index 2a15df243b5a4..53abbfeef3412 100644
--- a/libc/src/math/generic/common_constants.cpp
+++ b/libc/src/__support/math/common_constants.h
@@ -6,12 +6,29 @@
//
//===----------------------------------------------------------------------===//
-#include "common_constants.h"
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_COMMON_CONSTANTS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_COMMON_CONSTANTS_H
+
#include "src/__support/macros/config.h"
#include "src/__support/number_pair.h"
namespace LIBC_NAMESPACE_DECL {
+namespace common_constants_internal {
+
+// log(2) generated by Sollya with:
+// > a = 2^-43 * nearestint(2^43*log(2));
+// LSB = 2^-43 is chosen so that e_x * LOG_2_HI is exact for -1075 < e_x < 1024.
+static constexpr double LOG_2_HI = 0x1.62e42fefa38p-1; // LSB = 2^-43
+// > b = round(log10(2) - a, D, RN);
+static constexpr double LOG_2_LO = 0x1.ef35793c7673p-45; // LSB = 2^-97
+
+// Minimax polynomial for (log(1 + x) - x)/x^2, generated by sollya with:
+// > P = fpminimax((log(1 + x) - x)/x^2, 5, [|D...|], [-2^-8, 2^-7]);
+constexpr double LOG_COEFFS[6] = {-0x1.fffffffffffffp-2, 0x1.5555555554a9bp-2,
+ -0x1.0000000094567p-2, 0x1.99999dcc9823cp-3,
+ -0x1.55550ac2e537ap-3, 0x1.21a02c4e624d7p-3};
+
// Range reduction constants for logarithms.
// r(0) = 1, r(127) = 0.5
// r(k) = 2^-8 * ceil(2^8 * (1 - 2^-8) / (1 + k*2^-7))
@@ -19,7 +36,7 @@ namespace LIBC_NAMESPACE_DECL {
// precision, and -2^-8 <= v < 2^-7.
// TODO(lntue): Add reference to how the constants are derived after the
// resulting paper is ready.
-alignas(8) const float R[128] = {
+alignas(8) static constexpr float R[128] = {
0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.e8p-1,
0x1.e4p-1, 0x1.ep-1, 0x1.dep-1, 0x1.dap-1, 0x1.d6p-1, 0x1.d4p-1, 0x1.dp-1,
0x1.ccp-1, 0x1.cap-1, 0x1.c6p-1, 0x1.c4p-1, 0x1.cp-1, 0x1.bep-1, 0x1.bap-1,
@@ -40,7 +57,7 @@ alignas(8) const float R[128] = {
0x1.0ap-1, 0x1.08p-1, 0x1.08p-1, 0x1.06p-1, 0x1.06p-1, 0x1.04p-1, 0x1.04p-1,
0x1.02p-1, 0x1.0p-1};
-const double RD[128] = {
+static constexpr double RD[128] = {
0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.e8p-1,
0x1.e4p-1, 0x1.ep-1, 0x1.dep-1, 0x1.dap-1, 0x1.d6p-1, 0x1.d4p-1, 0x1.dp-1,
0x1.ccp-1, 0x1.cap-1, 0x1.c6p-1, 0x1.c4p-1, 0x1.cp-1, 0x1.bep-1, 0x1.bap-1,
@@ -65,7 +82,7 @@ const double RD[128] = {
// available.
// Generated by Sollya with the formula: CD[i] = RD[i]*(1 + i*2^-7) - 1
// for RD[i] defined on the table above.
-const double CD[128] = {
+static constexpr double CD[128] = {
0.0, -0x1p-14, -0x1p-12, -0x1.2p-11, -0x1p-10, -0x1.9p-10,
-0x1.2p-9, -0x1.88p-9, -0x1p-8, -0x1.9p-11, -0x1.fp-10, -0x1.9cp-9,
-0x1p-12, -0x1.cp-10, -0x1.bp-9, -0x1.5p-11, -0x1.4p-9, 0x1p-14,
@@ -90,7 +107,7 @@ const double CD[128] = {
-0x1p-14, -0x1p-8,
};
-const double LOG_R[128] = {
+static constexpr double LOG_R[128] = {
0x0.0000000000000p0, 0x1.010157588de71p-7, 0x1.0205658935847p-6,
0x1.8492528c8cabfp-6, 0x1.0415d89e74444p-5, 0x1.466aed42de3eap-5,
0x1.894aa149fb343p-5, 0x1.ccb73cdddb2ccp-5, 0x1.08598b59e3a07p-4,
@@ -135,7 +152,7 @@ const double LOG_R[128] = {
0x1.5707a26bb8c66p-1, 0x1.5af405c3649ep-1, 0x1.5af405c3649ep-1,
0x1.5ee82aa24192p-1, 0x0.000000000000p0};
-const double LOG2_R[128] = {
+static constexpr double LOG2_R[128] = {
0x0.0000000000000p+0, 0x1.72c7ba20f7327p-7, 0x1.743ee861f3556p-6,
0x1.184b8e4c56af8p-5, 0x1.77394c9d958d5p-5, 0x1.d6ebd1f1febfep-5,
0x1.1bb32a600549dp-4, 0x1.4c560fe68af88p-4, 0x1.7d60496cfbb4cp-4,
@@ -188,7 +205,7 @@ const double LOG2_R[128] = {
// print("{", -c, ",", -b, "},");
// };
// We replace LOG_R[0] with log10(1.0) == 0.0
-alignas(16) const NumberPair<double> LOG_R_DD[128] = {
+alignas(16) static constexpr NumberPair<double> LOG_R_DD[128] = {
{0.0, 0.0},
{-0x1.0c76b999d2be8p-46, 0x1.010157589p-7},
{-0x1.3dc5b06e2f7d2p-45, 0x1.0205658938p-6},
@@ -324,7 +341,7 @@ alignas(16) const NumberPair<double> LOG_R_DD[128] = {
// Output range:
// [-0x1.3ffcp-15, 0x1.3e3dp-15]
// We store S2[i] = 2^16 (r(i - 2^6) - 1).
-alignas(8) const int S2[193] = {
+alignas(8) static constexpr int S2[193] = {
0x101, 0xfd, 0xf9, 0xf5, 0xf1, 0xed, 0xe9, 0xe5, 0xe1,
0xdd, 0xd9, 0xd5, 0xd1, 0xcd, 0xc9, 0xc5, 0xc1, 0xbd,
0xb9, 0xb4, 0xb0, 0xac, 0xa8, 0xa4, 0xa0, 0x9c, 0x98,
@@ -348,7 +365,7 @@ alignas(8) const int S2[193] = {
-0x1cd, -0x1d1, -0x1d5, -0x1d9, -0x1dd, -0x1e0, -0x1e4, -0x1e8, -0x1ec,
-0x1f0, -0x1f4, -0x1f8, -0x1fc};
-const double R2[193] = {
+static constexpr double R2[193] = {
0x1.0101p0, 0x1.00fdp0, 0x1.00f9p0, 0x1.00f5p0, 0x1.00f1p0,
0x1.00edp0, 0x1.00e9p0, 0x1.00e5p0, 0x1.00e1p0, 0x1.00ddp0,
0x1.00d9p0, 0x1.00d5p0, 0x1.00d1p0, 0x1.00cdp0, 0x1.00c9p0,
@@ -395,7 +412,7 @@ const double R2[193] = {
// Output range:
// [-0x1.01928p-22 , 0x1p-22]
// We store S[i] = 2^21 (r(i - 80) - 1).
-alignas(8) const int S3[161] = {
+alignas(8) static constexpr int S3[161] = {
0x50, 0x4f, 0x4e, 0x4d, 0x4c, 0x4b, 0x4a, 0x49, 0x48, 0x47, 0x46,
0x45, 0x44, 0x43, 0x42, 0x41, 0x40, 0x3f, 0x3e, 0x3d, 0x3c, 0x3b,
0x3a, 0x39, 0x38, 0x37, 0x36, 0x35, 0x34, 0x33, 0x32, 0x31, 0x30,
@@ -418,7 +435,7 @@ alignas(8) const int S3[161] = {
// Output range:
// [-0x1.0002143p-29 , 0x1p-29]
// We store S[i] = 2^28 (r(i - 65) - 1).
-alignas(8) const int S4[130] = {
+alignas(8) static constexpr int S4[130] = {
0x41, 0x40, 0x3f, 0x3e, 0x3d, 0x3c, 0x3b, 0x3a, 0x39, 0x38, 0x37,
0x36, 0x35, 0x34, 0x33, 0x32, 0x31, 0x30, 0x2f, 0x2e, 0x2d, 0x2c,
0x2b, 0x2a, 0x29, 0x28, 0x27, 0x26, 0x25, 0x24, 0x23, 0x22, 0x21,
@@ -439,7 +456,7 @@ alignas(8) const int S4[130] = {
// Table is generated with Sollya as follow:
// > display = hexadecimal;
// > for i from -104 to 89 do { D(exp(i)); };
-const double EXP_M1[195] = {
+static constexpr double EXP_M1[195] = {
0x1.f1e6b68529e33p-151, 0x1.525be4e4e601dp-149, 0x1.cbe0a45f75eb1p-148,
0x1.3884e838aea68p-146, 0x1.a8c1f14e2af5dp-145, 0x1.20a717e64a9bdp-143,
0x1.8851d84118908p-142, 0x1.0a9bdfb02d240p-140, 0x1.6a5bea046b42ep-139,
@@ -511,7 +528,7 @@ const double EXP_M1[195] = {
// Table is generated with Sollya as follow:
// > display = hexadecimal;
// > for i from 0 to 127 do { D(exp(i / 128)); };
-const double EXP_M2[128] = {
+static constexpr double EXP_M2[128] = {
0x1.0000000000000p0, 0x1.0202015600446p0, 0x1.04080ab55de39p0,
0x1.06122436410ddp0, 0x1.08205601127edp0, 0x1.0a32a84e9c1f6p0,
0x1.0c49236829e8cp0, 0x1.0e63cfa7ab09dp0, 0x1.1082b577d34edp0,
@@ -557,4 +574,8 @@ const double EXP_M2[128] = {
0x1.568bb722dd593p1, 0x1.593b7d72305bbp1,
};
+} // namespace common_constants_internal
+
} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_COMMON_CONSTANTS_H
diff --git a/libc/src/__support/math/exp2.h b/libc/src/__support/math/exp2.h
new file mode 100644
index 0000000000000..7eaa465f93841
--- /dev/null
+++ b/libc/src/__support/math/exp2.h
@@ -0,0 +1,425 @@
+//===-- Implementation header for exp2 --------------------------*- 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_EXP2_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_H
+
+#include "common_constants.h" // Lookup tables EXP2_MID1 and EXP_M2.
+#include "exp_constants.h"
+#include "exp_utils.h" // ziv_test_denorm.
+#include "src/__support/CPP/bit.h"
+#include "src/__support/CPP/optional.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/dyadic_float.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/triple_double.h"
+#include "src/__support/common.h"
+#include "src/__support/integer_literals.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace exp2_internal {
+
+using namespace common_constants_internal;
+
+using fputil::DoubleDouble;
+using fputil::TripleDouble;
+using Float128 = typename fputil::DyadicFloat<128>;
+
+using LIBC_NAMESPACE::operator""_u128;
+
+// Error bounds:
+// Errors when using double precision.
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+constexpr double ERR_D = 0x1.0p-63;
+#else
+constexpr double ERR_D = 0x1.8p-63;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Errors when using double-double precision.
+constexpr double ERR_DD = 0x1.0p-100;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Polynomial approximations with double precision. Generated by Sollya with:
+// > P = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
+// > P;
+// Error bounds:
+// | output - (2^dx - 1) / dx | < 1.5 * 2^-52.
+LIBC_INLINE static double poly_approx_d(double dx) {
+ // dx^2
+ double dx2 = dx * dx;
+ double c0 =
+ fputil::multiply_add(dx, 0x1.ebfbdff82c58ep-3, 0x1.62e42fefa39efp-1);
+ double c1 =
+ fputil::multiply_add(dx, 0x1.3b2aba7a95a89p-7, 0x1.c6b08e8fc0c0ep-5);
+ double p = fputil::multiply_add(dx2, c1, c0);
+ return p;
+}
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Polynomial approximation with double-double precision. Generated by Solya
+// with:
+// > P = fpminimax((2^x - 1)/x, 5, [|DD...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
+// Error bounds:
+// | output - 2^(dx) | < 2^-101
+LIBC_INLINE static constexpr DoubleDouble
+poly_approx_dd(const DoubleDouble &dx) {
+ // Taylor polynomial.
+ constexpr DoubleDouble COEFFS[] = {
+ {0, 0x1p0},
+ {0x1.abc9e3b39824p-56, 0x1.62e42fefa39efp-1},
+ {-0x1.5e43a53e4527bp-57, 0x1.ebfbdff82c58fp-3},
+ {-0x1.d37963a9444eep-59, 0x1.c6b08d704a0cp-5},
+ {0x1.4eda1a81133dap-62, 0x1.3b2ab6fba4e77p-7},
+ {-0x1.c53fd1ba85d14p-64, 0x1.5d87fe7a265a5p-10},
+ {0x1.d89250b013eb8p-70, 0x1.430912f86cb8ep-13},
+ };
+
+ DoubleDouble p = fputil::polyeval(dx, COEFFS[0], COEFFS[1], COEFFS[2],
+ COEFFS[3], COEFFS[4], COEFFS[5], COEFFS[6]);
+ return p;
+}
+
+// Polynomial approximation with 128-bit precision:
+// Return exp(dx) ~ 1 + a0 * dx + a1 * dx^2 + ... + a6 * dx^7
+// For |dx| < 2^-13 + 2^-30:
+// | output - exp(dx) | < 2^-126.
+LIBC_INLINE static constexpr Float128 poly_approx_f128(const Float128 &dx) {
+ constexpr Float128 COEFFS_128[]{
+ {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0
+ {Sign::POS, -128, 0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128},
+ {Sign::POS, -128, 0x3d7f7bff'058b1d50'de2d60dd'9c9a1d9f_u128},
+ {Sign::POS, -132, 0xe35846b8'2505fc59'9d3b15d9'e7fb6897_u128},
+ {Sign::POS, -134, 0x9d955b7d'd273b94e'184462f6'bcd2b9e7_u128},
+ {Sign::POS, -137, 0xaec3ff3c'53398883'39ea1bb9'64c51a89_u128},
+ {Sign::POS, -138, 0x2861225f'345c396a'842c5341'8fa8ae61_u128},
+ {Sign::POS, -144, 0xffe5fe2d'109a319d'7abeb5ab'd5ad2079_u128},
+ };
+
+ Float128 p = fputil::polyeval(dx, COEFFS_128[0], COEFFS_128[1], COEFFS_128[2],
+ COEFFS_128[3], COEFFS_128[4], COEFFS_128[5],
+ COEFFS_128[6], COEFFS_128[7]);
+ return p;
+}
+
+// Compute 2^(x) using 128-bit precision.
+// TODO(lntue): investigate triple-double precision implementation for this
+// step.
+LIBC_INLINE static constexpr Float128 exp2_f128(double x, int hi, int idx1,
+ int idx2) {
+ Float128 dx = Float128(x);
+
+ // TODO: Skip recalculating exp_mid1 and exp_mid2.
+ Float128 exp_mid1 =
+ fputil::quick_add(Float128(EXP2_MID1[idx1].hi),
+ fputil::quick_add(Float128(EXP2_MID1[idx1].mid),
+ Float128(EXP2_MID1[idx1].lo)));
+
+ Float128 exp_mid2 =
+ fputil::quick_add(Float128(EXP2_MID2[idx2].hi),
+ fputil::quick_add(Float128(EXP2_MID2[idx2].mid),
+ Float128(EXP2_MID2[idx2].lo)));
+
+ Float128 exp_mid = fputil::quick_mul(exp_mid1, exp_mid2);
+
+ Float128 p = poly_approx_f128(dx);
+
+ Float128 r = fputil::quick_mul(exp_mid, p);
+
+ r.exponent += hi;
+
+ return r;
+}
+
+// Compute 2^x with double-double precision.
+LIBC_INLINE static DoubleDouble
+exp2_double_double(double x, const DoubleDouble &exp_mid) {
+ DoubleDouble dx({0, x});
+
+ // Degree-6 polynomial approximation in double-double precision.
+ // | p - 2^x | < 2^-103.
+ DoubleDouble p = poly_approx_dd(dx);
+
+ // Error bounds: 2^-102.
+ DoubleDouble r = fputil::quick_mult(exp_mid, p);
+
+ return r;
+}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// When output is denormal.
+LIBC_INLINE static double exp2_denorm(double x) {
+ // Range reduction.
+ int k =
+ static_cast<int>(cpp::bit_cast<uint64_t>(x + 0x1.8000'0000'4p21) >> 19);
+ double kd = static_cast<double>(k);
+
+ uint32_t idx1 = (k >> 6) & 0x3f;
+ uint32_t idx2 = k & 0x3f;
+
+ int hi = k >> 12;
+
+ DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi};
+ DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi};
+ DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2);
+
+ // |dx| < 2^-13 + 2^-30.
+ double dx = fputil::multiply_add(kd, -0x1.0p-12, x); // exact
+
+ double mid_lo = dx * exp_mid.hi;
+
+ // Approximate (2^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4.
+ double p = poly_approx_d(dx);
+
+ double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, exp_mid.hi, lo, ERR_D)
+ .value();
+#else
+ if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D);
+ LIBC_LIKELY(r.has_value()))
+ return r.value();
+
+ // Use double-double
+ DoubleDouble r_dd = exp2_double_double(dx, exp_mid);
+
+ if (auto r = ziv_test_denorm(hi, r_dd.hi, r_dd.lo, ERR_DD);
+ LIBC_LIKELY(r.has_value()))
+ return r.value();
+
+ // Use 128-bit precision
+ Float128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
+
+ return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+// Check for exceptional cases when:
+// * log2(1 - 2^-54) < x < log2(1 + 2^-53)
+// * x >= 1024
+// * x <= -1022
+// * x is inf or nan
+LIBC_INLINE static constexpr double set_exceptional(double x) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ uint64_t x_u = xbits.uintval();
+ uint64_t x_abs = xbits.abs().uintval();
+
+ // |x| < log2(1 + 2^-53)
+ if (x_abs <= 0x3ca71547652b82fd) {
+ // 2^(x) ~ 1 + x/2
+ return fputil::multiply_add(x, 0.5, 1.0);
+ }
+
+ // x <= -1022 || x >= 1024 or inf/nan.
+ if (x_u > 0xc08ff00000000000) {
+ // x <= -1075 or -inf/nan
+ if (x_u >= 0xc090cc0000000000) {
+ // exp(-Inf) = 0
+ if (xbits.is_inf())
+ return 0.0;
+
+ // exp(nan) = nan
+ if (xbits.is_nan())
+ return x;
+
+ if (fputil::quick_get_round() == FE_UPWARD)
+ return FPBits::min_subnormal().get_val();
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_UNDERFLOW);
+ return 0.0;
+ }
+
+ return exp2_denorm(x);
+ }
+
+ // x >= 1024 or +inf/nan
+ // x is finite
+ if (x_u < 0x7ff0'0000'0000'0000ULL) {
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+ return FPBits::max_normal().get_val();
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ }
+ // x is +inf or nan
+ return x + FPBits::inf().get_val();
+}
+
+} // namespace exp2_internal
+
+LIBC_INLINE static constexpr double exp2(double x) {
+ using namespace exp2_internal;
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ uint64_t x_u = xbits.uintval();
+
+ // x < -1022 or x >= 1024 or log2(1 - 2^-54) < x < log2(1 + 2^-53).
+ if (LIBC_UNLIKELY(x_u > 0xc08ff00000000000 ||
+ (x_u <= 0xbc971547652b82fe && x_u >= 0x4090000000000000) ||
+ x_u <= 0x3ca71547652b82fd)) {
+ return set_exceptional(x);
+ }
+
+ // Now -1075 < x <= log2(1 - 2^-54) or log2(1 + 2^-53) < x < 1024
+
+ // Range reduction:
+ // Let x = (hi + mid1 + mid2) + lo
+ // in which:
+ // hi is an integer
+ // mid1 * 2^6 is an integer
+ // mid2 * 2^12 is an integer
+ // then:
+ // 2^(x) = 2^hi * 2^(mid1) * 2^(mid2) * 2^(lo).
+ // With this formula:
+ // - multiplying by 2^hi is exact and cheap, simply by adding the exponent
+ // field.
+ // - 2^(mid1) and 2^(mid2) are stored in 2 x 64-element tables.
+ // - 2^(lo) ~ 1 + a0*lo + a1 * lo^2 + ...
+ //
+ // We compute (hi + mid1 + mid2) together by perform the rounding on x * 2^12.
+ // Since |x| < |-1075)| < 2^11,
+ // |x * 2^12| < 2^11 * 2^12 < 2^23,
+ // So we can fit the rounded result round(x * 2^12) in int32_t.
+ // Thus, the goal is to be able to use an additional addition and fixed width
+ // shift to get an int32_t representing round(x * 2^12).
+ //
+ // Assuming int32_t using 2-complement representation, since the mantissa part
+ // of a double precision is unsigned with the leading bit hidden, if we add an
+ // extra constant C = 2^e1 + 2^e2 with e1 > e2 >= 2^25 to the product, the
+ // part that are < 2^e2 in resulted mantissa of (x*2^12*L2E + C) can be
+ // considered as a proper 2-complement representations of x*2^12.
+ //
+ // One small problem with this approach is that the sum (x*2^12 + C) in
+ // double precision is rounded to the least significant bit of the dorminant
+ // factor C. In order to minimize the rounding errors from this addition, we
+ // want to minimize e1. Another constraint that we want is that after
+ // shifting the mantissa so that the least significant bit of int32_t
+ // corresponds to the unit bit of (x*2^12*L2E), the sign is correct without
+ // any adjustment. So combining these 2 requirements, we can choose
+ // C = 2^33 + 2^32, so that the sign bit corresponds to 2^31 bit, and hence
+ // after right shifting the mantissa, the resulting int32_t has correct sign.
+ // With this choice of C, the number of mantissa bits we need to shift to the
+ // right is: 52 - 33 = 19.
+ //
+ // Moreover, since the integer right shifts are equivalent to rounding down,
+ // we can add an extra 0.5 so that it will become round-to-nearest, tie-to-
+ // +infinity. So in particular, we can compute:
+ // hmm = x * 2^12 + C,
+ // where C = 2^33 + 2^32 + 2^-1, then if
+ // k = int32_t(lower 51 bits of double(x * 2^12 + C) >> 19),
+ // the reduced argument:
+ // lo = x - 2^-12 * k is bounded by:
+ // |lo| <= 2^-13 + 2^-12*2^-19
+ // = 2^-13 + 2^-31.
+ //
+ // Finally, notice that k only uses the mantissa of x * 2^12, so the
+ // exponent 2^12 is not needed. So we can simply define
+ // C = 2^(33 - 12) + 2^(32 - 12) + 2^(-13 - 12), and
+ // k = int32_t(lower 51 bits of double(x + C) >> 19).
+
+ // Rounding errors <= 2^-31.
+ int k =
+ static_cast<int>(cpp::bit_cast<uint64_t>(x + 0x1.8000'0000'4p21) >> 19);
+ double kd = static_cast<double>(k);
+
+ uint32_t idx1 = (k >> 6) & 0x3f;
+ uint32_t idx2 = k & 0x3f;
+
+ int hi = k >> 12;
+
+ DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi};
+ DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi};
+ DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2);
+
+ // |dx| < 2^-13 + 2^-30.
+ double dx = fputil::multiply_add(kd, -0x1.0p-12, x); // exact
+
+ // We use the degree-4 polynomial to approximate 2^(lo):
+ // 2^(lo) ~ 1 + a0 * lo + a1 * lo^2 + a2 * lo^3 + a3 * lo^4 = 1 + lo * P(lo)
+ // So that the errors are bounded by:
+ // |P(lo) - (2^lo - 1)/lo| < |lo|^4 / 64 < 2^(-13 * 4) / 64 = 2^-58
+ // Let P_ be an evaluation of P where all intermediate computations are in
+ // double precision. Using either Horner's or Estrin's schemes, the evaluated
+ // errors can be bounded by:
+ // |P_(lo) - P(lo)| < 2^-51
+ // => |lo * P_(lo) - (2^lo - 1) | < 2^-64
+ // => 2^(mid1 + mid2) * |lo * P_(lo) - expm1(lo)| < 2^-63.
+ // Since we approximate
+ // 2^(mid1 + mid2) ~ exp_mid.hi + exp_mid.lo,
+ // We use the expression:
+ // (exp_mid.hi + exp_mid.lo) * (1 + dx * P_(dx)) ~
+ // ~ exp_mid.hi + (exp_mid.hi * dx * P_(dx) + exp_mid.lo)
+ // with errors bounded by 2^-63.
+
+ double mid_lo = dx * exp_mid.hi;
+
+ // Approximate (2^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4.
+ double p = poly_approx_d(dx);
+
+ double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // To multiply by 2^hi, a fast way is to simply add hi to the exponent
+ // field.
+ int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
+ double r =
+ cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(exp_mid.hi + lo));
+ return r;
+#else
+ double upper = exp_mid.hi + (lo + ERR_D);
+ double lower = exp_mid.hi + (lo - ERR_D);
+
+ if (LIBC_LIKELY(upper == lower)) {
+ // To multiply by 2^hi, a fast way is to simply add hi to the exponent
+ // field.
+ int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
+ double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper));
+ return r;
+ }
+
+ // Use double-double
+ DoubleDouble r_dd = exp2_double_double(dx, exp_mid);
+
+ double upper_dd = r_dd.hi + (r_dd.lo + ERR_DD);
+ double lower_dd = r_dd.hi + (r_dd.lo - ERR_DD);
+
+ if (LIBC_LIKELY(upper_dd == lower_dd)) {
+ // To multiply by 2^hi, a fast way is to simply add hi to the exponent
+ // field.
+ int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
+ double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper_dd));
+ return r;
+ }
+
+ // Use 128-bit precision
+ Float128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
+
+ return static_cast<double>(r_f128);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 99c1b08326d53..28ea475f09656 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1448,21 +1448,7 @@ add_entrypoint_object(
HDRS
../exp2.h
DEPENDS
- .common_constants
- libc.src.__support.CPP.bit
- libc.src.__support.CPP.optional
- libc.src.__support.FPUtil.dyadic_float
- 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.triple_double
- libc.src.__support.integer_literals
- libc.src.__support.macros.optimization
- libc.src.__support.math.exp_utils
- libc.src.errno.errno
+ libc.src.__support.math.exp2
)
add_header_library(
@@ -1613,7 +1599,6 @@ add_entrypoint_object(
HDRS
../expm1.h
DEPENDS
- .common_constants
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.dyadic_float
libc.src.__support.FPUtil.fenv_impl
@@ -1624,6 +1609,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.triple_double
libc.src.__support.integer_literals
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
libc.src.errno.errno
)
@@ -1634,7 +1620,6 @@ add_entrypoint_object(
HDRS
../expm1f.h
DEPENDS
- .common_constants
libc.src.__support.FPUtil.basic_operations
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -1643,6 +1628,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.rounding_mode
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
libc.src.errno.errno
)
@@ -1673,7 +1659,6 @@ add_entrypoint_object(
HDRS
../powf.h
DEPENDS
- .common_constants
.exp2f_impl
libc.src.__support.math.exp10f
libc.src.__support.CPP.bit
@@ -1685,6 +1670,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.sqrt
libc.src.__support.FPUtil.triple_double
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
libc.src.errno.errno
)
@@ -1695,7 +1681,6 @@ add_entrypoint_object(
HDRS
../pow.h
DEPENDS
- .common_constants
libc.hdr.errno_macros
libc.hdr.fenv_macros
libc.src.__support.CPP.bit
@@ -1707,6 +1692,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2043,26 +2029,14 @@ add_entrypoint_object(
libc.src.__support.macros.properties.types
)
-add_object_library(
- common_constants
- HDRS
- common_constants.h
- SRCS
- common_constants.cpp
- DEPENDS
- libc.src.__support.math.exp_constants
- libc.src.__support.math.acosh_float_constants
- libc.src.__support.number_pair
-)
-
add_header_library(
log_range_reduction
HDRS
log_range_reduction.h
DEPENDS
- .common_constants
- libc.src.__support.uint128
libc.src.__support.FPUtil.dyadic_float
+ libc.src.__support.math.common_constants
+ libc.src.__support.uint128
)
add_entrypoint_object(
@@ -2072,7 +2046,6 @@ add_entrypoint_object(
HDRS
../log10.h
DEPENDS
- .common_constants
.log_range_reduction
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.dyadic_float
@@ -2082,6 +2055,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.integer_literals
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2091,12 +2065,12 @@ add_entrypoint_object(
HDRS
../log10f.h
DEPENDS
- .common_constants
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.fma
libc.src.__support.FPUtil.polyeval
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2126,7 +2100,6 @@ add_entrypoint_object(
HDRS
../log1p.h
DEPENDS
- .common_constants
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.dyadic_float
libc.src.__support.FPUtil.fenv_impl
@@ -2135,6 +2108,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.integer_literals
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2144,13 +2118,13 @@ add_entrypoint_object(
HDRS
../log1pf.h
DEPENDS
- .common_constants
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.fma
libc.src.__support.FPUtil.polyeval
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2160,7 +2134,6 @@ add_entrypoint_object(
HDRS
../log2.h
DEPENDS
- .common_constants
.log_range_reduction
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.dyadic_float
@@ -2170,6 +2143,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.integer_literals
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2179,13 +2153,13 @@ add_entrypoint_object(
HDRS
../log2f.h
DEPENDS
- .common_constants
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.fma
libc.src.__support.FPUtil.polyeval
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2215,7 +2189,6 @@ add_entrypoint_object(
HDRS
../log.h
DEPENDS
- .common_constants
.log_range_reduction
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.dyadic_float
@@ -2225,6 +2198,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.integer_literals
libc.src.__support.macros.optimization
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
@@ -2234,7 +2208,6 @@ add_entrypoint_object(
HDRS
../logf.h
DEPENDS
- .common_constants
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -2242,6 +2215,7 @@ add_entrypoint_object(
libc.src.__support.FPUtil.polyeval
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.cpu_features
+ libc.src.__support.math.common_constants
)
add_entrypoint_object(
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
deleted file mode 100644
index 9ee31f09baa8a..0000000000000
--- a/libc/src/math/generic/common_constants.h
+++ /dev/null
@@ -1,73 +0,0 @@
-//===-- Common constants for math functions ---------------------*- 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_MATH_GENERIC_COMMON_CONSTANTS_H
-#define LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H
-
-#include "src/__support/FPUtil/triple_double.h"
-#include "src/__support/macros/config.h"
-#include "src/__support/math/acosh_float_constants.h"
-#include "src/__support/math/exp_constants.h"
-#include "src/__support/number_pair.h"
-
-namespace LIBC_NAMESPACE_DECL {
-
-// Lookup table for range reduction constants r for logarithms.
-extern const float R[128];
-
-// Lookup table for range reduction constants r for logarithms.
-extern const double RD[128];
-
-// Lookup table for compensated constants for exact range reduction when FMA
-// instructions are not available.
-extern const double CD[128];
-
-// Lookup table for -log(r)
-extern const double LOG_R[128];
-extern const NumberPair<double> LOG_R_DD[128];
-
-// Lookup table for -log2(r)
-extern const double LOG2_R[128];
-
-// Minimax polynomial for (log(1 + x) - x)/x^2, generated by sollya with:
-// > P = fpminimax((log(1 + x) - x)/x^2, 5, [|D...|], [-2^-8, 2^-7]);
-constexpr double LOG_COEFFS[6] = {-0x1.fffffffffffffp-2, 0x1.5555555554a9bp-2,
- -0x1.0000000094567p-2, 0x1.99999dcc9823cp-3,
- -0x1.55550ac2e537ap-3, 0x1.21a02c4e624d7p-3};
-
-// Logarithm Range Reduction - Step 2, 3, and 4.
-extern const int S2[193];
-extern const int S3[161];
-extern const int S4[130];
-
-extern const double R2[193];
-
-// log(2) generated by Sollya with:
-// > a = 2^-43 * nearestint(2^43*log(2));
-// LSB = 2^-43 is chosen so that e_x * LOG_2_HI is exact for -1075 < e_x < 1024.
-constexpr double LOG_2_HI = 0x1.62e42fefa38p-1; // LSB = 2^-43
-// > b = round(log10(2) - a, D, RN);
-constexpr double LOG_2_LO = 0x1.ef35793c7673p-45; // LSB = 2^-97
-
-// Lookup table for exp(m) with m = -104, ..., 89.
-// -104 = floor(log(single precision's min denormal))
-// 89 = ceil(log(single precision's max normal))
-// Table is generated with Sollya as follow:
-// > display = hexadecimal;
-// > for i from -104 to 89 do { D(exp(i)); };
-extern const double EXP_M1[195];
-
-// Lookup table for exp(m * 2^(-7)) with m = 0, ..., 127.
-// Table is generated with Sollya as follow:
-// > display = hexadecimal;
-// > for i from 0 to 127 do { D(exp(i / 128)); };
-extern const double EXP_M2[128];
-
-} // namespace LIBC_NAMESPACE_DECL
-
-#endif // LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H
diff --git a/libc/src/math/generic/exp2.cpp b/libc/src/math/generic/exp2.cpp
index 154154f2b90c3..20e1ff5fa43a7 100644
--- a/libc/src/math/generic/exp2.cpp
+++ b/libc/src/math/generic/exp2.cpp
@@ -7,404 +7,10 @@
//===----------------------------------------------------------------------===//
#include "src/math/exp2.h"
-#include "common_constants.h" // Lookup tables EXP2_MID1 and EXP_M2.
-#include "src/__support/CPP/bit.h"
-#include "src/__support/CPP/optional.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/FPUtil/double_double.h"
-#include "src/__support/FPUtil/dyadic_float.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/triple_double.h"
-#include "src/__support/common.h"
-#include "src/__support/integer_literals.h"
-#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "src/__support/math/exp_utils.h" // ziv_test_denorm.
+#include "src/__support/math/exp2.h"
namespace LIBC_NAMESPACE_DECL {
-using fputil::DoubleDouble;
-using fputil::TripleDouble;
-using Float128 = typename fputil::DyadicFloat<128>;
-
-using LIBC_NAMESPACE::operator""_u128;
-
-// Error bounds:
-// Errors when using double precision.
-#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-constexpr double ERR_D = 0x1.0p-63;
-#else
-constexpr double ERR_D = 0x1.8p-63;
-#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Errors when using double-double precision.
-constexpr double ERR_DD = 0x1.0p-100;
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-namespace {
-
-// Polynomial approximations with double precision. Generated by Sollya with:
-// > P = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
-// > P;
-// Error bounds:
-// | output - (2^dx - 1) / dx | < 1.5 * 2^-52.
-LIBC_INLINE double poly_approx_d(double dx) {
- // dx^2
- double dx2 = dx * dx;
- double c0 =
- fputil::multiply_add(dx, 0x1.ebfbdff82c58ep-3, 0x1.62e42fefa39efp-1);
- double c1 =
- fputil::multiply_add(dx, 0x1.3b2aba7a95a89p-7, 0x1.c6b08e8fc0c0ep-5);
- double p = fputil::multiply_add(dx2, c1, c0);
- return p;
-}
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Polynomial approximation with double-double precision. Generated by Solya
-// with:
-// > P = fpminimax((2^x - 1)/x, 5, [|DD...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
-// Error bounds:
-// | output - 2^(dx) | < 2^-101
-DoubleDouble poly_approx_dd(const DoubleDouble &dx) {
- // Taylor polynomial.
- constexpr DoubleDouble COEFFS[] = {
- {0, 0x1p0},
- {0x1.abc9e3b39824p-56, 0x1.62e42fefa39efp-1},
- {-0x1.5e43a53e4527bp-57, 0x1.ebfbdff82c58fp-3},
- {-0x1.d37963a9444eep-59, 0x1.c6b08d704a0cp-5},
- {0x1.4eda1a81133dap-62, 0x1.3b2ab6fba4e77p-7},
- {-0x1.c53fd1ba85d14p-64, 0x1.5d87fe7a265a5p-10},
- {0x1.d89250b013eb8p-70, 0x1.430912f86cb8ep-13},
- };
-
- DoubleDouble p = fputil::polyeval(dx, COEFFS[0], COEFFS[1], COEFFS[2],
- COEFFS[3], COEFFS[4], COEFFS[5], COEFFS[6]);
- return p;
-}
-
-// Polynomial approximation with 128-bit precision:
-// Return exp(dx) ~ 1 + a0 * dx + a1 * dx^2 + ... + a6 * dx^7
-// For |dx| < 2^-13 + 2^-30:
-// | output - exp(dx) | < 2^-126.
-Float128 poly_approx_f128(const Float128 &dx) {
- constexpr Float128 COEFFS_128[]{
- {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0
- {Sign::POS, -128, 0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128},
- {Sign::POS, -128, 0x3d7f7bff'058b1d50'de2d60dd'9c9a1d9f_u128},
- {Sign::POS, -132, 0xe35846b8'2505fc59'9d3b15d9'e7fb6897_u128},
- {Sign::POS, -134, 0x9d955b7d'd273b94e'184462f6'bcd2b9e7_u128},
- {Sign::POS, -137, 0xaec3ff3c'53398883'39ea1bb9'64c51a89_u128},
- {Sign::POS, -138, 0x2861225f'345c396a'842c5341'8fa8ae61_u128},
- {Sign::POS, -144, 0xffe5fe2d'109a319d'7abeb5ab'd5ad2079_u128},
- };
-
- Float128 p = fputil::polyeval(dx, COEFFS_128[0], COEFFS_128[1], COEFFS_128[2],
- COEFFS_128[3], COEFFS_128[4], COEFFS_128[5],
- COEFFS_128[6], COEFFS_128[7]);
- return p;
-}
-
-// Compute 2^(x) using 128-bit precision.
-// TODO(lntue): investigate triple-double precision implementation for this
-// step.
-Float128 exp2_f128(double x, int hi, int idx1, int idx2) {
- Float128 dx = Float128(x);
-
- // TODO: Skip recalculating exp_mid1 and exp_mid2.
- Float128 exp_mid1 =
- fputil::quick_add(Float128(EXP2_MID1[idx1].hi),
- fputil::quick_add(Float128(EXP2_MID1[idx1].mid),
- Float128(EXP2_MID1[idx1].lo)));
-
- Float128 exp_mid2 =
- fputil::quick_add(Float128(EXP2_MID2[idx2].hi),
- fputil::quick_add(Float128(EXP2_MID2[idx2].mid),
- Float128(EXP2_MID2[idx2].lo)));
-
- Float128 exp_mid = fputil::quick_mul(exp_mid1, exp_mid2);
-
- Float128 p = poly_approx_f128(dx);
-
- Float128 r = fputil::quick_mul(exp_mid, p);
-
- r.exponent += hi;
-
- return r;
-}
-
-// Compute 2^x with double-double precision.
-DoubleDouble exp2_double_double(double x, const DoubleDouble &exp_mid) {
- DoubleDouble dx({0, x});
-
- // Degree-6 polynomial approximation in double-double precision.
- // | p - 2^x | < 2^-103.
- DoubleDouble p = poly_approx_dd(dx);
-
- // Error bounds: 2^-102.
- DoubleDouble r = fputil::quick_mult(exp_mid, p);
-
- return r;
-}
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-// When output is denormal.
-double exp2_denorm(double x) {
- // Range reduction.
- int k =
- static_cast<int>(cpp::bit_cast<uint64_t>(x + 0x1.8000'0000'4p21) >> 19);
- double kd = static_cast<double>(k);
-
- uint32_t idx1 = (k >> 6) & 0x3f;
- uint32_t idx2 = k & 0x3f;
-
- int hi = k >> 12;
-
- DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi};
- DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi};
- DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2);
-
- // |dx| < 2^-13 + 2^-30.
- double dx = fputil::multiply_add(kd, -0x1.0p-12, x); // exact
-
- double mid_lo = dx * exp_mid.hi;
-
- // Approximate (2^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4.
- double p = poly_approx_d(dx);
-
- double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- return ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, exp_mid.hi, lo, ERR_D)
- .value();
-#else
- if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D);
- LIBC_LIKELY(r.has_value()))
- return r.value();
-
- // Use double-double
- DoubleDouble r_dd = exp2_double_double(dx, exp_mid);
-
- if (auto r = ziv_test_denorm(hi, r_dd.hi, r_dd.lo, ERR_DD);
- LIBC_LIKELY(r.has_value()))
- return r.value();
-
- // Use 128-bit precision
- Float128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
-
- return static_cast<double>(r_f128);
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-}
-
-// Check for exceptional cases when:
-// * log2(1 - 2^-54) < x < log2(1 + 2^-53)
-// * x >= 1024
-// * x <= -1022
-// * x is inf or nan
-double set_exceptional(double x) {
- using FPBits = typename fputil::FPBits<double>;
- FPBits xbits(x);
-
- uint64_t x_u = xbits.uintval();
- uint64_t x_abs = xbits.abs().uintval();
-
- // |x| < log2(1 + 2^-53)
- if (x_abs <= 0x3ca71547652b82fd) {
- // 2^(x) ~ 1 + x/2
- return fputil::multiply_add(x, 0.5, 1.0);
- }
-
- // x <= -1022 || x >= 1024 or inf/nan.
- if (x_u > 0xc08ff00000000000) {
- // x <= -1075 or -inf/nan
- if (x_u >= 0xc090cc0000000000) {
- // exp(-Inf) = 0
- if (xbits.is_inf())
- return 0.0;
-
- // exp(nan) = nan
- if (xbits.is_nan())
- return x;
-
- if (fputil::quick_get_round() == FE_UPWARD)
- return FPBits::min_subnormal().get_val();
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_UNDERFLOW);
- return 0.0;
- }
-
- return exp2_denorm(x);
- }
-
- // x >= 1024 or +inf/nan
- // x is finite
- if (x_u < 0x7ff0'0000'0000'0000ULL) {
- int rounding = fputil::quick_get_round();
- if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
- return FPBits::max_normal().get_val();
-
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_OVERFLOW);
- }
- // x is +inf or nan
- return x + FPBits::inf().get_val();
-}
-
-} // namespace
-
-LLVM_LIBC_FUNCTION(double, exp2, (double x)) {
- using FPBits = typename fputil::FPBits<double>;
- FPBits xbits(x);
-
- uint64_t x_u = xbits.uintval();
-
- // x < -1022 or x >= 1024 or log2(1 - 2^-54) < x < log2(1 + 2^-53).
- if (LIBC_UNLIKELY(x_u > 0xc08ff00000000000 ||
- (x_u <= 0xbc971547652b82fe && x_u >= 0x4090000000000000) ||
- x_u <= 0x3ca71547652b82fd)) {
- return set_exceptional(x);
- }
-
- // Now -1075 < x <= log2(1 - 2^-54) or log2(1 + 2^-53) < x < 1024
-
- // Range reduction:
- // Let x = (hi + mid1 + mid2) + lo
- // in which:
- // hi is an integer
- // mid1 * 2^6 is an integer
- // mid2 * 2^12 is an integer
- // then:
- // 2^(x) = 2^hi * 2^(mid1) * 2^(mid2) * 2^(lo).
- // With this formula:
- // - multiplying by 2^hi is exact and cheap, simply by adding the exponent
- // field.
- // - 2^(mid1) and 2^(mid2) are stored in 2 x 64-element tables.
- // - 2^(lo) ~ 1 + a0*lo + a1 * lo^2 + ...
- //
- // We compute (hi + mid1 + mid2) together by perform the rounding on x * 2^12.
- // Since |x| < |-1075)| < 2^11,
- // |x * 2^12| < 2^11 * 2^12 < 2^23,
- // So we can fit the rounded result round(x * 2^12) in int32_t.
- // Thus, the goal is to be able to use an additional addition and fixed width
- // shift to get an int32_t representing round(x * 2^12).
- //
- // Assuming int32_t using 2-complement representation, since the mantissa part
- // of a double precision is unsigned with the leading bit hidden, if we add an
- // extra constant C = 2^e1 + 2^e2 with e1 > e2 >= 2^25 to the product, the
- // part that are < 2^e2 in resulted mantissa of (x*2^12*L2E + C) can be
- // considered as a proper 2-complement representations of x*2^12.
- //
- // One small problem with this approach is that the sum (x*2^12 + C) in
- // double precision is rounded to the least significant bit of the dorminant
- // factor C. In order to minimize the rounding errors from this addition, we
- // want to minimize e1. Another constraint that we want is that after
- // shifting the mantissa so that the least significant bit of int32_t
- // corresponds to the unit bit of (x*2^12*L2E), the sign is correct without
- // any adjustment. So combining these 2 requirements, we can choose
- // C = 2^33 + 2^32, so that the sign bit corresponds to 2^31 bit, and hence
- // after right shifting the mantissa, the resulting int32_t has correct sign.
- // With this choice of C, the number of mantissa bits we need to shift to the
- // right is: 52 - 33 = 19.
- //
- // Moreover, since the integer right shifts are equivalent to rounding down,
- // we can add an extra 0.5 so that it will become round-to-nearest, tie-to-
- // +infinity. So in particular, we can compute:
- // hmm = x * 2^12 + C,
- // where C = 2^33 + 2^32 + 2^-1, then if
- // k = int32_t(lower 51 bits of double(x * 2^12 + C) >> 19),
- // the reduced argument:
- // lo = x - 2^-12 * k is bounded by:
- // |lo| <= 2^-13 + 2^-12*2^-19
- // = 2^-13 + 2^-31.
- //
- // Finally, notice that k only uses the mantissa of x * 2^12, so the
- // exponent 2^12 is not needed. So we can simply define
- // C = 2^(33 - 12) + 2^(32 - 12) + 2^(-13 - 12), and
- // k = int32_t(lower 51 bits of double(x + C) >> 19).
-
- // Rounding errors <= 2^-31.
- int k =
- static_cast<int>(cpp::bit_cast<uint64_t>(x + 0x1.8000'0000'4p21) >> 19);
- double kd = static_cast<double>(k);
-
- uint32_t idx1 = (k >> 6) & 0x3f;
- uint32_t idx2 = k & 0x3f;
-
- int hi = k >> 12;
-
- DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi};
- DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi};
- DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2);
-
- // |dx| < 2^-13 + 2^-30.
- double dx = fputil::multiply_add(kd, -0x1.0p-12, x); // exact
-
- // We use the degree-4 polynomial to approximate 2^(lo):
- // 2^(lo) ~ 1 + a0 * lo + a1 * lo^2 + a2 * lo^3 + a3 * lo^4 = 1 + lo * P(lo)
- // So that the errors are bounded by:
- // |P(lo) - (2^lo - 1)/lo| < |lo|^4 / 64 < 2^(-13 * 4) / 64 = 2^-58
- // Let P_ be an evaluation of P where all intermediate computations are in
- // double precision. Using either Horner's or Estrin's schemes, the evaluated
- // errors can be bounded by:
- // |P_(lo) - P(lo)| < 2^-51
- // => |lo * P_(lo) - (2^lo - 1) | < 2^-64
- // => 2^(mid1 + mid2) * |lo * P_(lo) - expm1(lo)| < 2^-63.
- // Since we approximate
- // 2^(mid1 + mid2) ~ exp_mid.hi + exp_mid.lo,
- // We use the expression:
- // (exp_mid.hi + exp_mid.lo) * (1 + dx * P_(dx)) ~
- // ~ exp_mid.hi + (exp_mid.hi * dx * P_(dx) + exp_mid.lo)
- // with errors bounded by 2^-63.
-
- double mid_lo = dx * exp_mid.hi;
-
- // Approximate (2^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4.
- double p = poly_approx_d(dx);
-
- double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // To multiply by 2^hi, a fast way is to simply add hi to the exponent
- // field.
- int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
- double r =
- cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(exp_mid.hi + lo));
- return r;
-#else
- double upper = exp_mid.hi + (lo + ERR_D);
- double lower = exp_mid.hi + (lo - ERR_D);
-
- if (LIBC_LIKELY(upper == lower)) {
- // To multiply by 2^hi, a fast way is to simply add hi to the exponent
- // field.
- int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
- double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper));
- return r;
- }
-
- // Use double-double
- DoubleDouble r_dd = exp2_double_double(dx, exp_mid);
-
- double upper_dd = r_dd.hi + (r_dd.lo + ERR_DD);
- double lower_dd = r_dd.hi + (r_dd.lo - ERR_DD);
-
- if (LIBC_LIKELY(upper_dd == lower_dd)) {
- // To multiply by 2^hi, a fast way is to simply add hi to the exponent
- // field.
- int64_t exp_hi = static_cast<int64_t>(hi) << FPBits::FRACTION_LEN;
- double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper_dd));
- return r;
- }
-
- // Use 128-bit precision
- Float128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
-
- return static_cast<double>(r_f128);
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-}
+LLVM_LIBC_FUNCTION(double, exp2, (double x)) { return math::exp2(x); }
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/expm1.cpp b/libc/src/math/generic/expm1.cpp
index c360554a8af8f..a3d0c1aa5261c 100644
--- a/libc/src/math/generic/expm1.cpp
+++ b/libc/src/math/generic/expm1.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/expm1.h"
-#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
@@ -22,6 +21,8 @@
#include "src/__support/integer_literals.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/math/common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
+#include "src/__support/math/exp_constants.h"
#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
#define LIBC_MATH_EXPM1_SKIP_ACCURATE_PASS
@@ -59,6 +60,8 @@ constexpr double MLOG_2_EXP2_M12_LO = 0x1.b0e2633fe0685p-79;
namespace {
+using namespace common_constants_internal;
+
// Polynomial approximations with double precision:
// Return expm1(dx) / x ~ 1 + dx / 2 + dx^2 / 6 + dx^3 / 24.
// For |dx| < 2^-13 + 2^-30:
diff --git a/libc/src/math/generic/expm1f.cpp b/libc/src/math/generic/expm1f.cpp
index b2967e2516197..72c8aa358d618 100644
--- a/libc/src/math/generic/expm1f.cpp
+++ b/libc/src/math/generic/expm1f.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/expm1f.h"
-#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "src/__support/FPUtil/BasicOperations.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FMA.h"
@@ -20,10 +19,12 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+#include "src/__support/math/common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
+ using namespace common_constants_internal;
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
diff --git a/libc/src/math/generic/log.cpp b/libc/src/math/generic/log.cpp
index 0cd4424ee0baf..66ce0591e78ef 100644
--- a/libc/src/math/generic/log.cpp
+++ b/libc/src/math/generic/log.cpp
@@ -18,8 +18,8 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "common_constants.h"
#include "log_range_reduction.h"
+#include "src/__support/math/common_constants.h"
namespace LIBC_NAMESPACE_DECL {
@@ -30,6 +30,8 @@ using LIBC_NAMESPACE::operator""_u128;
namespace {
+using namespace common_constants_internal;
+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// A simple upper bound for the error of e_x * log(2) - log(r).
constexpr double HI_ERR = 0x1.0p-85;
diff --git a/libc/src/math/generic/log10.cpp b/libc/src/math/generic/log10.cpp
index 1c4e559ba083c..95f24fac01453 100644
--- a/libc/src/math/generic/log10.cpp
+++ b/libc/src/math/generic/log10.cpp
@@ -18,8 +18,8 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "common_constants.h"
#include "log_range_reduction.h"
+#include "src/__support/math/common_constants.h"
namespace LIBC_NAMESPACE_DECL {
@@ -30,6 +30,8 @@ using LIBC_NAMESPACE::operator""_u128;
namespace {
+using namespace common_constants_internal;
+
constexpr fputil::DoubleDouble LOG10_E = {0x1.95355baaafad3p-57,
0x1.bcb7b1526e50ep-2};
@@ -739,6 +741,7 @@ double log10_accurate(int e_x, int index, double m_x) {
} // namespace
LLVM_LIBC_FUNCTION(double, log10, (double x)) {
+ using namespace common_constants_internal;
using FPBits_t = typename fputil::FPBits<double>;
FPBits_t xbits(x);
diff --git a/libc/src/math/generic/log10f.cpp b/libc/src/math/generic/log10f.cpp
index 81e7cdbcfe5a1..6b9cc5d3e0dbf 100644
--- a/libc/src/math/generic/log10f.cpp
+++ b/libc/src/math/generic/log10f.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/log10f.h"
-#include "common_constants.h" // Lookup table for (1/f)
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FMA.h"
#include "src/__support/FPUtil/FPBits.h"
@@ -18,6 +17,7 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/common_constants.h" // Lookup table for (1/f)
// This is an algorithm for log10(x) in single precision which is
// correctly rounded for all rounding modes, based on the implementation of
@@ -104,6 +104,7 @@ static constexpr double LOG10_R[128] = {
0x1.30cb3a7bb3625p-2, 0x1.34413509f79ffp-2};
LLVM_LIBC_FUNCTION(float, log10f, (float x)) {
+ using namespace common_constants_internal;
constexpr double LOG10_2 = 0x1.34413509f79ffp-2;
using FPBits = typename fputil::FPBits<float>;
diff --git a/libc/src/math/generic/log1p.cpp b/libc/src/math/generic/log1p.cpp
index 09f465a6ba774..1595981ac309b 100644
--- a/libc/src/math/generic/log1p.cpp
+++ b/libc/src/math/generic/log1p.cpp
@@ -18,7 +18,7 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "common_constants.h"
+#include "src/__support/math/common_constants.h"
namespace LIBC_NAMESPACE_DECL {
@@ -29,6 +29,8 @@ using LIBC_NAMESPACE::operator""_u128;
namespace {
+using namespace common_constants_internal;
+
// R1[i] = 2^-8 * nearestint( 2^8 / (1 + i * 2^-7) )
constexpr double R1[129] = {
0x1p0, 0x1.fcp-1, 0x1.f8p-1, 0x1.f4p-1, 0x1.fp-1, 0x1.ecp-1, 0x1.eap-1,
diff --git a/libc/src/math/generic/log1pf.cpp b/libc/src/math/generic/log1pf.cpp
index 16b1b34a6c944..f0289c2320a24 100644
--- a/libc/src/math/generic/log1pf.cpp
+++ b/libc/src/math/generic/log1pf.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/log1pf.h"
-#include "common_constants.h" // Lookup table for (1/f) and log(f)
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FMA.h"
#include "src/__support/FPUtil/FPBits.h"
@@ -18,6 +17,8 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/acosh_float_constants.h"
+#include "src/__support/math/common_constants.h" // Lookup table for (1/f) and log(f)
// This is an algorithm for log10(x) in single precision which is
// correctly rounded for all rounding modes.
@@ -38,6 +39,7 @@ namespace internal {
// We don't need to treat denormal and 0
LIBC_INLINE float log(double x) {
using namespace acoshf_internal;
+ using namespace common_constants_internal;
constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
using FPBits = typename fputil::FPBits<double>;
diff --git a/libc/src/math/generic/log2.cpp b/libc/src/math/generic/log2.cpp
index 27ca2fc350f17..f0c0ae302592f 100644
--- a/libc/src/math/generic/log2.cpp
+++ b/libc/src/math/generic/log2.cpp
@@ -18,8 +18,8 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
-#include "common_constants.h"
#include "log_range_reduction.h"
+#include "src/__support/math/common_constants.h"
namespace LIBC_NAMESPACE_DECL {
@@ -30,6 +30,8 @@ using LIBC_NAMESPACE::operator""_u128;
namespace {
+using namespace common_constants_internal;
+
constexpr fputil::DoubleDouble LOG2_E = {0x1.777d0ffda0d24p-56,
0x1.71547652b82fep0};
@@ -859,6 +861,7 @@ double log2_accurate(int e_x, int index, double m_x) {
} // namespace
LLVM_LIBC_FUNCTION(double, log2, (double x)) {
+ using namespace common_constants_internal;
using FPBits_t = typename fputil::FPBits<double>;
FPBits_t xbits(x);
diff --git a/libc/src/math/generic/log2f.cpp b/libc/src/math/generic/log2f.cpp
index cff718eec2169..7353f0302ab3e 100644
--- a/libc/src/math/generic/log2f.cpp
+++ b/libc/src/math/generic/log2f.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/log2f.h"
-#include "common_constants.h" // Lookup table for (1/f)
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
@@ -15,7 +14,8 @@
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/math/common_constants.h" // Lookup table for (1/f)
// This is a correctly-rounded algorithm for log2(x) in single precision with
// round-to-nearest, tie-to-even mode from the RLIBM project at:
@@ -55,6 +55,7 @@
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(float, log2f, (float x)) {
+ using namespace common_constants_internal;
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
diff --git a/libc/src/math/generic/log_range_reduction.h b/libc/src/math/generic/log_range_reduction.h
index 8c94230ee7577..7484506ddf68f 100644
--- a/libc/src/math/generic/log_range_reduction.h
+++ b/libc/src/math/generic/log_range_reduction.h
@@ -9,9 +9,9 @@
#ifndef LLVM_LIBC_SRC_MATH_GENERIC_LOG_RANGE_REDUCTION_H
#define LLVM_LIBC_SRC_MATH_GENERIC_LOG_RANGE_REDUCTION_H
-#include "common_constants.h"
#include "src/__support/FPUtil/dyadic_float.h"
#include "src/__support/macros/config.h"
+#include "src/__support/math/common_constants.h"
#include "src/__support/uint128.h"
namespace LIBC_NAMESPACE_DECL {
@@ -36,6 +36,7 @@ struct LogRR {
LIBC_INLINE fputil::DyadicFloat<128>
log_range_reduction(double m_x, const LogRR &log_table,
fputil::DyadicFloat<128> &sum) {
+ using namespace common_constants_internal;
using Float128 = typename fputil::DyadicFloat<128>;
using MType = typename Float128::MantissaType;
diff --git a/libc/src/math/generic/logf.cpp b/libc/src/math/generic/logf.cpp
index e8d2ba2cfe175..4d2947dc14db0 100644
--- a/libc/src/math/generic/logf.cpp
+++ b/libc/src/math/generic/logf.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/logf.h"
-#include "common_constants.h" // Lookup table for (1/f) and log(f)
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
@@ -17,6 +16,7 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/common_constants.h" // Lookup table for (1/f) and log(f)
// This is an algorithm for log(x) in single precision which is correctly
// rounded for all rounding modes, based on the implementation of log(x) from
@@ -53,6 +53,7 @@
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(float, logf, (float x)) {
+ using namespace common_constants_internal;
constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
using FPBits = typename fputil::FPBits<float>;
diff --git a/libc/src/math/generic/pow.cpp b/libc/src/math/generic/pow.cpp
index 43e99a7acf690..c9f685b82fcb5 100644
--- a/libc/src/math/generic/pow.cpp
+++ b/libc/src/math/generic/pow.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/math/pow.h"
-#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
@@ -21,6 +20,8 @@
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/math/common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
+#include "src/__support/math/exp_constants.h" // Lookup tables EXP_M1 and EXP_M2.
namespace LIBC_NAMESPACE_DECL {
@@ -28,6 +29,8 @@ using fputil::DoubleDouble;
namespace {
+using namespace common_constants_internal;
+
// Constants for log2(x) range reduction, generated by Sollya with:
// > for i from 0 to 127 do {
// r = 2^-8 * ceil( 2^8 * (1 - 2^(-8)) / (1 + i*2^-7) );
diff --git a/libc/src/math/generic/powf.cpp b/libc/src/math/generic/powf.cpp
index a45ef511c9bad..12246e9878758 100644
--- a/libc/src/math/generic/powf.cpp
+++ b/libc/src/math/generic/powf.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#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/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
@@ -15,10 +14,13 @@
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/sqrt.h" // Speedup for powf(x, 1/2) = sqrtf(x)
+#include "src/__support/FPUtil/triple_double.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/math/common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "src/__support/math/exp10f.h" // Speedup for powf(10, y) = exp10f(y)
+#include "src/__support/math/exp_constants.h"
#include "exp2f_impl.h" // Speedup for powf(2, y) = exp2f(y)
@@ -29,6 +31,8 @@ using fputil::TripleDouble;
namespace {
+using namespace common_constants_internal;
+
#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
alignas(16) constexpr DoubleDouble LOG2_R_DD[128] = {
{0.0, 0.0},
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index ea4634cbe7f9f..040f6356b2577 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -40,6 +40,7 @@ add_fp_unittest(
libc.src.__support.math.exp10m1f16
libc.src.__support.math.erff
libc.src.__support.math.exp
+ libc.src.__support.math.exp2
libc.src.__support.math.exp10
libc.src.__support.math.exp10f
libc.src.__support.math.exp10f16
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 17221932927b0..ef2e7b8d71456 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -80,6 +80,7 @@ TEST(LlvmLibcSharedMathTest, AllDouble) {
EXPECT_FP_EQ(0x1p+0, LIBC_NAMESPACE::shared::cos(0.0));
EXPECT_FP_EQ(0x0p+0, LIBC_NAMESPACE::shared::dsqrtl(0.0));
EXPECT_FP_EQ(0x1p+0, LIBC_NAMESPACE::shared::exp(0.0));
+ EXPECT_FP_EQ(0x1p+0, LIBC_NAMESPACE::shared::exp2(0.0));
EXPECT_FP_EQ(0x1p+0, LIBC_NAMESPACE::shared::exp10(0.0));
}
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index e57d9dea036dd..7107c85e9768f 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -2156,25 +2156,14 @@ libc_function(
########################### math support library ###############################
-libc_support_library(
- name = "common_constants",
- srcs = ["src/math/generic/common_constants.cpp"],
- hdrs = ["src/math/generic/common_constants.h"],
- deps = [
- ":__support_math_acosh_float_constants",
- ":__support_math_exp_constants",
- ":__support_number_pair",
- ],
-)
-
libc_support_library(
name = "log_range_reduction",
hdrs = ["src/math/generic/log_range_reduction.h"],
deps = [
":__support_common",
":__support_fputil_dyadic_float",
+ ":__support_math_common_constants",
":__support_uint128",
- ":common_constants",
],
)
@@ -2189,8 +2178,8 @@ libc_support_library(
":__support_fputil_polyeval",
":__support_fputil_rounding_mode",
":__support_macros_optimization",
+ ":__support_math_common_constants",
":__support_math_exp10f_utils",
- ":common_constants",
],
)
@@ -2557,6 +2546,16 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_common_constants",
+ hdrs = ["src/__support/math/common_constants.h"],
+ deps = [
+ ":__support_math_acosh_float_constants",
+ ":__support_math_exp_constants",
+ ":__support_number_pair",
+ ],
+)
+
libc_support_library(
name = "__support_math_cos",
hdrs = ["src/__support/math/cos.h"],
@@ -2632,8 +2631,8 @@ libc_support_library(
":__support_fputil_polyeval",
":__support_fputil_rounding_mode",
":__support_macros_optimization",
+ ":__support_math_common_constants",
":__support_sincosf_utils",
- ":common_constants",
],
)
@@ -2878,6 +2877,24 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_exp2",
+ hdrs = ["src/__support/math/exp2.h"],
+ deps = [
+ ":__support_fputil_double_double",
+ ":__support_fputil_dyadic_float",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer",
+ ":__support_fputil_polyeval",
+ ":__support_fputil_rounding_mode",
+ ":__support_fputil_triple_double",
+ ":__support_integer_literals",
+ ":__support_macros_optimization",
+ ":__support_math_common_constants",
+ ":__support_math_exp_utils",
+ ],
+)
+
libc_support_library(
name = "__support_math_exp10",
hdrs = ["src/__support/math/exp10.h"],
@@ -3652,17 +3669,7 @@ libc_math_function(
libc_math_function(
name = "exp2",
additional_deps = [
- ":__support_fputil_double_double",
- ":__support_fputil_dyadic_float",
- ":__support_fputil_multiply_add",
- ":__support_fputil_nearest_integer",
- ":__support_fputil_polyeval",
- ":__support_fputil_rounding_mode",
- ":__support_fputil_triple_double",
- ":__support_integer_literals",
- ":__support_macros_optimization",
- ":__support_math_exp_utils",
- ":common_constants",
+ ":__support_math_exp2",
],
)
@@ -3706,7 +3713,7 @@ libc_math_function(
":__support_fputil_triple_double",
":__support_integer_literals",
":__support_macros_optimization",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -3720,7 +3727,7 @@ libc_math_function(
":__support_fputil_rounding_mode",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4233,7 +4240,7 @@ libc_math_function(
":__support_integer_literals",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
":log_range_reduction",
],
)
@@ -4246,7 +4253,7 @@ libc_math_function(
":__support_fputil_polyeval",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4268,7 +4275,7 @@ libc_math_function(
":__support_integer_literals",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
":log_range_reduction",
],
)
@@ -4281,7 +4288,7 @@ libc_math_function(
":__support_fputil_polyeval",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4303,7 +4310,7 @@ libc_math_function(
":__support_integer_literals",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4315,7 +4322,7 @@ libc_math_function(
":__support_fputil_polyeval",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4330,7 +4337,7 @@ libc_math_function(
":__support_integer_literals",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":common_constants",
+ ":__support_math_common_constants",
":log_range_reduction",
],
)
@@ -4342,7 +4349,7 @@ libc_math_function(
":__support_fputil_multiply_add",
":__support_fputil_polyeval",
":__support_macros_optimization",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4488,7 +4495,7 @@ libc_math_function(
":__support_fputil_nearest_integer",
":__support_fputil_polyeval",
":__support_fputil_sqrt",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4503,7 +4510,7 @@ libc_math_function(
":__support_fputil_triple_double",
":__support_macros_optimization",
":__support_math_exp10f",
- ":common_constants",
+ ":__support_math_common_constants",
":exp2f_impl",
],
)
@@ -4664,7 +4671,7 @@ libc_math_function(
":__support_fputil_rounding_mode",
":__support_macros_optimization",
":__support_math_sinhfcoshf_utils",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
@@ -4771,7 +4778,7 @@ libc_math_function(
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
":__support_math_exp10f_utils",
- ":common_constants",
+ ":__support_math_common_constants",
],
)
More information about the llvm-commits
mailing list