[libc-commits] [libc] [libc][math] Make pow function correctly rounded for all rounding modes. (PR #222827)
via libc-commits
libc-commits at lists.llvm.org
Thu Sep 10 19:18:10 PDT 2026
llvmorg-github-actions[bot] wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-libc
Author: lntue
<details>
<summary>Changes</summary>
Fix overflow issues reported by Paul Zimmermann.
Algorithm overview:
Evaluate mainly as:
```
x^y = 2^(y * log2(x))
```
- Fast path:
- Compute `log2(x) = e_x + log2(m_x)` in double-double precision.
- Scale `y * log2(x)` and evaluate `2^(y * log2(x)) = 2^hi * 2^mid * 2^lo`
using a 64-entry lookup table for `2^mid` and a degree-5 polynomial
evaluated with Estrin's scheme for `2^lo`.
- Perform Ziv's rounding test.
- 128-bit accurate path (pow_accurate_128.h):
- If the fast path Ziv's test fails, recompute `log2` and `exp2` using
fixed-point arithmetic with Frac128 and Frac64.
- Evaluate degree-15 polynomial for `log2` and degree-11 polynomial for
exp2 with minimax coefficients generated by Sollya.
- Use Horner scheme in 2 steps: higher degree terms are evaluated
with 64-bit precision before combining into 128-bit precision.
- Detect exact and midpoint rounding boundaries using algorithm described in
Lauter, C. and Lefevre, V., "Rounding Boundary Cases for Values of the
Exponential and Power Functions," IEEE Trans. Comput. 58(8):1063-1074.
- 256-bit accurate path (pow_accurate_256.h):
- If the 128-bit result fails the Ziv's test, recompute using 256-bit fixed-point
arithmetic with Frac256.
- Evaluate degree-31 polynomial for `log2` and degree-21 polynomial for
`exp2` using Horner scheme in 3 steps (64-bit -> 128-bit -> 256-bit).
- Share the 64-entry table `EXP2_MID_FRAC256` with the 128-bit path.
- Add `Frac256` class in `libc/src/__support/frac256.h` for unsigned 256-bit
fixed-point arithmetic.
Assisted-by: Gemini was used for debugging, performance analysis, and refactoring.
---
Patch is 145.86 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/222827.diff
14 Files Affected:
- (modified) libc/docs/headers/math/index.rst (+1-1)
- (modified) libc/src/__support/CMakeLists.txt (+11)
- (modified) libc/src/__support/frac128.h (+26-1)
- (added) libc/src/__support/frac256.h (+100)
- (modified) libc/src/__support/frac64.h (+8)
- (modified) libc/src/__support/macros/attributes.h (+8)
- (modified) libc/src/__support/math/CMakeLists.txt (+12)
- (modified) libc/src/__support/math/pow.h (+346-428)
- (added) libc/src/__support/math/pow_accurate_128.h (+375)
- (added) libc/src/__support/math/pow_accurate_256.h (+818)
- (added) libc/src/__support/math/pow_fast.h (+304)
- (added) libc/src/__support/math/pow_utils.h (+671)
- (modified) libc/test/src/math/pow_test.cpp (+25-5)
- (modified) libc/test/src/math/smoke/pow_test.cpp (+24-23)
``````````diff
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index 2d5e188e804d4..67fc2e259eff9 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -331,7 +331,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
| logp1 | | | | | | | 7.12.6.14 | F.10.3.14 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
-| pow | |check| | 1 ULP | | | | | 7.12.7.5 | F.10.4.5 |
+| pow | |check| | |check| | | | | | 7.12.7.5 | F.10.4.5 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
| powi\* | | | | | | | | |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt
index b93c88260cdc9..f6d998fee2ef6 100644
--- a/libc/src/__support/CMakeLists.txt
+++ b/libc/src/__support/CMakeLists.txt
@@ -392,6 +392,17 @@ add_header_library(
frac128.h
DEPENDS
.big_int
+ .frac64
+ libc.src.__support.macros.config
+)
+
+add_header_library(
+ frac256
+ HDRS
+ frac256.h
+ DEPENDS
+ .big_int
+ .frac128
libc.src.__support.macros.config
)
diff --git a/libc/src/__support/frac128.h b/libc/src/__support/frac128.h
index 298d1d9a61d70..68b337ef50445 100644
--- a/libc/src/__support/frac128.h
+++ b/libc/src/__support/frac128.h
@@ -1,15 +1,21 @@
-//===-- 128-bit unsigned fractional type -------------------*- 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
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the declaration of 128-bit unsigned fractional type.
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC128_H
#define LLVM_LIBC_SRC___SUPPORT_FRAC128_H
#include "big_int.h"
+#include "frac64.h"
#include "src/__support/macros/config.h"
namespace LIBC_NAMESPACE_DECL {
@@ -17,6 +23,17 @@ namespace LIBC_NAMESPACE_DECL {
struct Frac128 : public UInt<128> {
using UInt<128>::UInt;
+ // Convert Frac128 number to Frac64 with round-to-nearest
+ // (using bit 63 as the rounding bit).
+ LIBC_INLINE constexpr explicit operator Frac64() const {
+ uint64_t round = val[0] >> 63;
+ return Frac64(val[1] + round);
+ }
+
+ LIBC_INLINE constexpr Frac64 to_frac64() const {
+ return static_cast<Frac64>(*this);
+ }
+
LIBC_INLINE constexpr Frac128 operator~() const {
Frac128 r{};
r.val[0] = ~val[0];
@@ -53,6 +70,14 @@ struct Frac128 : public UInt<128> {
*this = *this * other;
return *this;
}
+
+ LIBC_INLINE constexpr Frac128 operator<<(size_t s) const {
+ return Frac128((UInt<128>(*this) << s).val);
+ }
+
+ LIBC_INLINE constexpr Frac128 operator>>(size_t s) const {
+ return Frac128((UInt<128>(*this) >> s).val);
+ }
};
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/__support/frac256.h b/libc/src/__support/frac256.h
new file mode 100644
index 0000000000000..ed11ffe1f9d58
--- /dev/null
+++ b/libc/src/__support/frac256.h
@@ -0,0 +1,100 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the declaration of 256-bit unsigned fractional type.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC256_H
+#define LLVM_LIBC_SRC___SUPPORT_FRAC256_H
+
+#include "big_int.h"
+#include "frac128.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+struct Frac256 : public UInt<256> {
+ using UInt<256>::UInt;
+
+ // Convert Frac256 number to Frac128 with round-to-nearest
+ // (using bit 127 as the rounding bit).
+ LIBC_INLINE constexpr explicit operator Frac128() const {
+ uint64_t round = val[1] >> 63;
+ uint64_t lo = val[2] + round;
+ uint64_t hi = val[3] + (lo < round ? 1 : 0);
+ return Frac128({lo, hi});
+ }
+
+ LIBC_INLINE constexpr Frac128 to_frac128() const {
+ return static_cast<Frac128>(*this);
+ }
+
+ // Convert Frac256 number to Frac64 with round-to-nearest
+ // (using bit 191 as the rounding bit, i.e., bit 63 of limb 2).
+ LIBC_INLINE constexpr explicit operator Frac64() const {
+ uint64_t round = val[2] >> 63;
+ return Frac64(val[3] + round);
+ }
+
+ LIBC_INLINE constexpr Frac64 to_frac64() const {
+ return static_cast<Frac64>(*this);
+ }
+
+ LIBC_INLINE constexpr Frac256 operator~() const {
+ Frac256 r{};
+ r.val[0] = ~val[0];
+ r.val[1] = ~val[1];
+ r.val[2] = ~val[2];
+ r.val[3] = ~val[3];
+ return r;
+ }
+
+ LIBC_INLINE constexpr Frac256 operator+(const Frac256 &other) const {
+ UInt<256> r = UInt<256>(*this) + (UInt<256>(other));
+ return Frac256(r.val);
+ }
+
+ LIBC_INLINE constexpr Frac256 operator-(const Frac256 &other) const {
+ UInt<256> r = UInt<256>(*this) - (UInt<256>(other));
+ return Frac256(r.val);
+ }
+
+ LIBC_INLINE constexpr Frac256 operator*(const Frac256 &other) const {
+ UInt<256> r = UInt<256>::quick_mul_hi(UInt<256>(other));
+ return Frac256(r.val);
+ }
+
+ LIBC_INLINE constexpr Frac256 &operator+=(const Frac256 &other) {
+ *this = *this + other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac256 &operator-=(const Frac256 &other) {
+ *this = *this - other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac256 &operator*=(const Frac256 &other) {
+ *this = *this * other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac256 operator<<(size_t s) const {
+ return Frac256((UInt<256>(*this) << s).val);
+ }
+
+ LIBC_INLINE constexpr Frac256 operator>>(size_t s) const {
+ return Frac256((UInt<256>(*this) >> s).val);
+ }
+};
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_FRAC256_H
diff --git a/libc/src/__support/frac64.h b/libc/src/__support/frac64.h
index bbb0d3bb30903..350b0a8b701b5 100644
--- a/libc/src/__support/frac64.h
+++ b/libc/src/__support/frac64.h
@@ -51,6 +51,14 @@ struct Frac64 : public UInt<64> {
*this = *this * other;
return *this;
}
+
+ LIBC_INLINE constexpr Frac64 operator<<(size_t s) const {
+ return Frac64(val[0] << s);
+ }
+
+ LIBC_INLINE constexpr Frac64 operator>>(size_t s) const {
+ return Frac64(val[0] >> s);
+ }
};
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/__support/macros/attributes.h b/libc/src/__support/macros/attributes.h
index cffcabcd29bd3..8fb56659d28b3 100644
--- a/libc/src/__support/macros/attributes.h
+++ b/libc/src/__support/macros/attributes.h
@@ -29,6 +29,14 @@
#define LIBC_INLINE_ASM __asm__ __volatile__
#define LIBC_UNUSED __attribute__((unused))
+#if __has_attribute(always_inline)
+#define LIBC_ALWAYS_INLINE LIBC_INLINE __attribute__((always_inline))
+#elif defined(LIBC_COMPILER_IS_MSVC) || defined(_MSC_VER)
+#define LIBC_ALWAYS_INLINE __forceinline
+#else
+#define LIBC_ALWAYS_INLINE LIBC_INLINE
+#endif
+
#ifndef LIBC_HAS_BUILTIN_IS_CONSTANT_EVALUATED
#if (defined(LIBC_COMPILER_IS_GCC) && (LIBC_COMPILER_GCC_VER >= 900)) || \
(defined(LIBC_COMPILER_IS_CLANG) && LIBC_COMPILER_CLANG_VER >= 900)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 0871a5bf229b6..2115e390d8ef7 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -5269,20 +5269,32 @@ add_header_library(
pow
HDRS
pow.h
+ pow_fast.h
+ pow_utils.h
+ pow_accurate_128.h
+ pow_accurate_256.h
DEPENDS
.common_constants
.exp_constants
+ .exp2
libc.hdr.errno_macros
libc.hdr.fenv_macros
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.double_double
+ 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.sqrt
+ libc.src.__support.frac64
+ libc.src.__support.frac128
+ libc.src.__support.frac256
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+ libc.src.__support.uint128
)
add_header_library(
diff --git a/libc/src/__support/math/pow.h b/libc/src/__support/math/pow.h
index c03e6be271d05..8d21ff2f28362 100644
--- a/libc/src/__support/math/pow.h
+++ b/libc/src/__support/math/pow.h
@@ -1,368 +1,168 @@
-//===-- Implementation header for pow ---------------------------*- 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
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Implementation header for double-precision pow(x, y).
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_H
-#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
-#include "exp_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"
-#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/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
-#include "src/__support/FPUtil/sqrt.h" // Speedup for pow(x, 1/2) = sqrt(x)
#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"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/common_constants.h"
+#include "src/__support/math/exp_constants.h"
+#include "src/__support/math/pow_utils.h"
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#include "src/__support/math/pow_fast.h"
+#else
+
+#include "src/__support/math/pow_accurate_128.h"
+#include "src/__support/math/pow_accurate_256.h"
namespace LIBC_NAMESPACE_DECL {
namespace math {
-namespace pow_internal {
-
-using fputil::DoubleDouble;
-
-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) );
-// b = nearestint(log2(r) * 2^41) * 2^-41;
-// c = round(log2(r) - b, D, RN);
-// print("{", -c, ",", -b, "},");
-// };
-// This is the same as -log2(RD[i]), with the least significant bits of the
-// high part set to be 2^-41, so that the sum of high parts + e_x is exact in
-// double precision.
-// We also replace the first and the last ones to be 0.
-LIBC_INLINE_VAR constexpr DoubleDouble LOG2_R_DD[128] = {
- {0.0, 0.0},
- {-0x1.19b14945cf6bap-44, 0x1.72c7ba21p-7},
- {-0x1.95539356f93dcp-43, 0x1.743ee862p-6},
- {0x1.abe0a48f83604p-43, 0x1.184b8e4c5p-5},
- {0x1.635577970e04p-43, 0x1.77394c9d9p-5},
- {-0x1.401fbaaa67e3cp-45, 0x1.d6ebd1f2p-5},
- {-0x1.5b1799ceaeb51p-43, 0x1.1bb32a6008p-4},
- {0x1.7c407050799bfp-43, 0x1.4c560fe688p-4},
- {0x1.da6339da288fcp-43, 0x1.7d60496cf8p-4},
- {0x1.be4f6f22dbbadp-43, 0x1.960caf9ab8p-4},
- {-0x1.c760bc9b188c4p-45, 0x1.c7b528b71p-4},
- {0x1.164e932b2d51cp-44, 0x1.f9c95dc1dp-4},
- {0x1.924ae921f7ecap-45, 0x1.097e38ce6p-3},
- {-0x1.6d25a5b8a19b2p-44, 0x1.22dadc2ab4p-3},
- {0x1.e50a1644ac794p-43, 0x1.3c6fb650ccp-3},
- {0x1.f34baa74a7942p-43, 0x1.494f863b8cp-3},
- {-0x1.8f7aac147fdc1p-46, 0x1.633a8bf438p-3},
- {0x1.f84be19cb9578p-43, 0x1.7046031c78p-3},
- {-0x1.66cccab240e9p-46, 0x1.8a8980abfcp-3},
- {-0x1.3f7a55cd2af4cp-47, 0x1.97c1cb13c8p-3},
- {0x1.3458cde69308cp-43, 0x1.b2602497d4p-3},
- {-0x1.667f21fa8423fp-44, 0x1.bfc67a8p-3},
- {0x1.d2fe4574e09b9p-47, 0x1.dac22d3e44p-3},
- {0x1.367bde40c5e6dp-43, 0x1.e857d3d36p-3},
- {0x1.d45da26510033p-46, 0x1.01d9bbcfa6p-2},
- {-0x1.7204f55bbf90dp-44, 0x1.08bce0d96p-2},
- {-0x1.d4f1b95e0ff45p-43, 0x1.169c05364p-2},
- {0x1.c20d74c0211bfp-44, 0x1.1d982c9d52p-2},
- {0x1.ad89a083e072ap-43, 0x1.249cd2b13cp-2},
- {0x1.cd0cb4492f1bcp-43, 0x1.32bfee370ep-2},
- {-0x1.2101a9685c779p-47, 0x1.39de8e155ap-2},
- {0x1.9451cd394fe8dp-43, 0x1.4106017c3ep-2},
- {0x1.661e393a16b95p-44, 0x1.4f6fbb2cecp-2},
- {-0x1.c6d8d86531d56p-44, 0x1.56b22e6b58p-2},
- {0x1.c1c885adb21d3p-43, 0x1.5dfdcf1eeap-2},
- {0x1.3bb5921006679p-45, 0x1.6552b49986p-2},
- {0x1.1d406db502403p-43, 0x1.6cb0f6865cp-2},
- {0x1.55a63e278bad5p-43, 0x1.7b89f02cf2p-2},
- {-0x1.66ae2a7ada553p-49, 0x1.8304d90c12p-2},
- {-0x1.66cccab240e9p-45, 0x1.8a8980abfcp-2},
- {-0x1.62404772a151dp-45, 0x1.921800924ep-2},
- {0x1.ac9bca36fd02ep-44, 0x1.99b072a96cp-2},
- {0x1.4bc302ffa76fbp-43, 0x1.a8ff97181p-2},
- {0x1.01fea1ec47c71p-43, 0x1.b0b67f4f46p-2},
- {-0x1.f20203b3186a6p-43, 0x1.b877c57b1cp-2},
- {-0x1.2642415d47384p-45, 0x1.c043859e3p-2},
- {-0x1.bc76a2753b99bp-50, 0x1.c819dc2d46p-2},
- {-0x1.da93ae3a5f451p-43, 0x1.cffae611aep-2},
- {-0x1.50e785694a8c6p-43, 0x1.d7e6c0abc4p-2},
- {0x1.c56138c894641p-43, 0x1.dfdd89d586p-2},
- {0x1.5669df6a2b592p-43, 0x1.e7df5fe538p-2},
- {-0x1.ea92d9e0e8ac2p-48, 0x1.efec61b012p-2},
- {0x1.a0331af2e6feap-43, 0x1.f804ae8d0cp-2},
- {0x1.9518ce032f41dp-48, 0x1.0014332bep-1},
- {-0x1.b3b3864c60011p-44, 0x1.042bd4b9a8p-1},
- {-0x1.103e8f00d41c8p-45, 0x1.08494c66b9p-1},
- {0x1.65be75cc3da17p-43, 0x1.0c6caaf0c5p-1},
- {0x1.3676289cd3dd4p-43, 0x1.1096015deep-1},
- {-0x1.41dfc7d7c3321p-43, 0x1.14c560fe69p-1},
- {0x1.e0cda8bd74461p-44, 0x1.18fadb6e2dp-1},
- {0x1.2a606046ad444p-44, 0x1.1d368296b5p-1},
- {0x1.f9ea977a639cp-43, 0x1.217868b0c3p-1},
- {-0x1.50520a377c7ecp-45, 0x1.25c0a0463cp-1},
- {0x1.6e3cb71b554e7p-47, 0x1.2a0f3c3407p-1},
- {-0x1.4275f1035e5e8p-48, 0x1.2e644fac05p-1},
- {-0x1.4275f1035e5e8p-48, 0x1.2e644fac05p-1},
- {-0x1.979a5db68721dp-45, 0x1.32bfee370fp-1},
- {0x1.1ee969a95f529p-43, 0x1.37222bb707p-1},
- {0x1.bb4b69336b66ep-43, 0x1.3b8b1c68fap-1},
- {0x1.d5e6a8a4fb059p-45, 0x1.3ffad4e74fp-1},
- {0x1.3106e404cabb7p-44, 0x1.44716a2c08p-1},
- {0x1.3106e404cabb7p-44, 0x1.44716a2c08p-1},
- {-0x1.9bcaf1aa4168ap-43, 0x1.48eef19318p-1},
- {0x1.1646b761c48dep-44, 0x1.4d7380dcc4p-1},
- {0x1.2f0c0bfe9dbecp-43, 0x1.51ff2e3021p-1},
- {0x1.29904613e33cp-43, 0x1.5692101d9bp-1},
- {0x1.1d406db502403p-44, 0x1.5b2c3da197p-1},
- {0x1.1d406db502403p-44, 0x1.5b2c3da197p-1},
- {-0x1.125d6cbcd1095p-44, 0x1.5fcdce2728p-1},
- {-0x1.bd9b32266d92cp-43, 0x1.6476d98adap-1},
- {0x1.54243b21709cep-44, 0x1.6927781d93p-1},
- {0x1.54243b21709cep-44, 0x1.6927781d93p-1},
- {-0x1.ce60916e52e91p-44, 0x1.6ddfc2a79p-1},
- {0x1.f1f5ae718f241p-43, 0x1.729fd26b7p-1},
- {-0x1.6eb9612e0b4f3p-43, 0x1.7767c12968p-1},
- {-0x1.6eb9612e0b4f3p-43, 0x1.7767c12968p-1},
- {0x1.fed21f9cb2cc5p-43, 0x1.7c37a9227ep-1},
- {0x1.7f5dc57266758p-43, 0x1.810fa51bf6p-1},
- {0x1.7f5dc57266758p-43, 0x1.810fa51bf6p-1},
- {0x1.5b338360c2ae2p-43, 0x1.85efd062c6p-1},
- {-0x1.96fc8f4b56502p-43, 0x1.8ad846cf37p-1},
- {-0x1.96fc8f4b56502p-43, 0x1.8ad846cf37p-1},
- {-0x1.bdc81c4db3134p-44, 0x1.8fc924c89bp-1},
- {0x1.36c101ee1344p-43, 0x1.94c287492cp-1},
- {0x1.36c101ee1344p-43, 0x1.94c287492cp-1},
- {0x1.e41fa0a62e6aep-44, 0x1.99c48be206p-1},
- {-0x1.d97ee9124773bp-46, 0x1.9ecf50bf44p-1},
- {-0x1.d97ee9124773bp-46, 0x1.9ecf50bf44p-1},
- {-0x1.3f94e00e7d6bcp-46, 0x1.a3e2f4ac44p-1},
- {-0x1.6879fa00b120ap-43, 0x1.a8ff971811p-1},
- {-0x1.6879fa00b120ap-43, 0x1.a8ff971811p-1},
- {0x1.1659d8e2d7d38p-44, 0x1.ae255819fp-1},
- {0x1.1e5e0ae0d3f8ap-43, 0x1.b35458761dp-1},
- {0x1.1e5e0ae0d3f8ap-43, 0x1.b35458761dp-1},
- {0x1.484a15babcf88p-43, 0x1.b88cb9a2abp-1},
- {0x1.484a15babcf88p-43, 0x1.b88cb9a2abp-1},
- {0x1.871a7610e40bdp-45, 0x1.bdce9dcc96p-1},
- {-0x1.2d90e5edaeceep-43, 0x1.c31a27dd01p-1},
- {-0x1.2d90e5edaeceep-43, 0x1.c31a27dd01p-1},
- {-0x1.5dd31d962d373p-43, 0x1.c86f7b7ea5p-1},
- {-0x1.5dd31d962d373p-43, 0x1.c86f7b7ea5p-1},
- {-0x1.9ad57391924a7p-43, 0x1.cdcebd2374p-1},
- {-0x1.3167ccc538261p-44, 0x1.d338120a6ep-1},
- {-0x1.3167ccc538261p-44, 0x1.d338120a6ep-1},
- {0x1.c7a4ff65ddbc9p-45, 0x1.d8aba045bp-1},
- {0x1.c7a4ff65ddbc9p-45, 0x1.d8aba045bp-1},
- {-0x1.f9ab3cf74babap-44, 0x1.de298ec0bbp-1},
- {-0x1.f9ab3cf74babap-44, 0x1.de298ec0bbp-1},
- {0x1.52842c1c1e586p-43, 0x1.e3b20546f5p-1},
- {0x1.52842c1c1e586p-43, 0x1.e3b20546f5p-1},
- {0x1.3c6764fc87b4ap-48, 0x1.e9452c8a71p-1},
- {0x1.3c6764fc87b4ap-48, 0x1.e9452c8a71p-1},
- {-0x1.a0976c0a2827dp-44, 0x1.eee32e2aedp-1},
- {-0x1.a0976c0a2827dp-44, 0x1.eee32e2aedp-1},
- {-0x1.a45314dc4fc42p-43, 0x1.f48c34bd1fp-1},
- {-0x1.a45314dc4fc42p-43, 0x1.f48c34bd1fp-1},
- {0x1.ef5d00e390ap-44, 0x1.fa406bd244p-1},
- {0.0, 1.0},
-};
-
-LIBC_INLINE bool is_odd_integer(double x) {
- using FPBits = fputil::FPBits<double>;
- FPBits xbits(x);
- uint64_t x_u = xbits.uintval();
- unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
- unsigned lsb =
- static_cast<unsigned>(cpp::countr_zero(x_u | FPBits::EXP_MASK));
- constexpr unsigned UNIT_EXPONENT =
- static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
- return (x_e + lsb == UNIT_EXPONENT);
-}
-
-LIBC_INLINE bool is_integer(double x) {
- using FPBits = fputil::FPBits<double>;
- FPBits xbits(x);
- uint64_t x_u = xbits.uintval();
- unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
- unsigned lsb =
- static_cast<unsigned>(cpp::countr_zero(x_u | FPBits::EXP_MASK));
- constexpr unsigned UNIT_EXPONENT =
- static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
- return (x_e + lsb >= UNIT_EXPONENT);
-}
-
-} // namespace pow_internal
+// Overview of the main part of pow(x, y) = x^y computations.
+//
+// Let x = 2^(e_x) * m_x > 0. Then:
+// x^y = 2^( y * log2(x) )
+// = 2^( y * ( e_x + log2(m_x) ) )
+// = 2^( e_h + e_l )
+// = 2^(e_h) * 2^(e_l)
+// where:
+// e_h = round(y * log2(x)),
+// e_l = {y * log2(x)} = y * log2(x) - e_h.
+//
+// In particular, e_h is an integer, and |e_l| <= 0.5.
+//
+// For the final result to ...
[truncated]
``````````
</details>
https://github.com/llvm/llvm-project/pull/222827
More information about the libc-commits
mailing list