[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:32:42 PDT 2026
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/222827
>From e1e2c4deb65f3d3ab2d1fb30d2a375c889edc135 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 11 Sep 2026 01:58:20 +0000
Subject: [PATCH 1/2] [libc][math] Make pow function correctly rounded for all
rounding modes.
---
libc/docs/headers/math/index.rst | 2 +-
libc/src/__support/CMakeLists.txt | 11 +
libc/src/__support/frac128.h | 27 +-
libc/src/__support/frac256.h | 100 +++
libc/src/__support/frac64.h | 8 +
libc/src/__support/macros/attributes.h | 8 +
libc/src/__support/math/CMakeLists.txt | 12 +
libc/src/__support/math/pow.h | 774 +++++++++----------
libc/src/__support/math/pow_accurate_128.h | 375 ++++++++++
libc/src/__support/math/pow_accurate_256.h | 818 +++++++++++++++++++++
libc/src/__support/math/pow_fast.h | 304 ++++++++
libc/src/__support/math/pow_utils.h | 671 +++++++++++++++++
libc/test/src/math/pow_test.cpp | 30 +-
libc/test/src/math/smoke/pow_test.cpp | 47 +-
14 files changed, 2729 insertions(+), 458 deletions(-)
create mode 100644 libc/src/__support/frac256.h
create mode 100644 libc/src/__support/math/pow_accurate_128.h
create mode 100644 libc/src/__support/math/pow_accurate_256.h
create mode 100644 libc/src/__support/math/pow_fast.h
create mode 100644 libc/src/__support/math/pow_utils.h
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 be fit in double precision, the exponent field can be
+// bounded by:
+// -1075 <= e_h <= 1024,
+// and anything outside that can be passed to quick overflow/underflow logic.
+//
+// Since |e_l| <= 0.5:
+// 0.5 < 2^(e_l) < 2,
+// the relative error of x^y = 2^(e_h) * 2^(e_l) is about
+// ~ absolute error of e_l
+// ~ absolute error of y * log2(x)
+// ~ relative error(log2(x)) * |y|
+// ~ relative error(log2(x)) * |e_h|
+// < relative error(log2(x)) * 2^11.
+//
+// Roughly speaking, to compute x^y with relative error < 2^(-n), we will need
+// to compute log2(x) with relative error < 2^(-n - 11), approximately.
+//
+// To compute log2(x), we will perform range reduction for log2(m_x):
+// dx = r * m_x - 1.
+// Then m_x = (1 + dx) / r, and
+// log2(m_x) = log2( (1 + dx) / r )
+// = log2(1 + dx) - log2(r),
+// where -log2(r) is obtained from look up tables.
+// The computation of dx = r * m_x - 1 is exact, and we choose the size of the
+// look up table for r's such that:
+// -2^-8 <= dx < 2^-7.
+// Combining them together, we have that:
+// log2(x) = e_x - log2(r) + log2(1 + dx)
+// = e_x - log2(r) + log2(e) * (dx - dx^2/2 + dx^3/3 - ...)
+// So in the worst case, where e_x = log2(r) = 0, the relative error of log2(x)
+// computation will be:
+// ~ absolute_error(log2(1 + dx)) / |dx|.
+// So if we compute log2(1 + dx) accurately up to dx^n term, and use a
+// polynomial approximation for the dx^(n + 1) and higher terms:
+// log2(1 + dx) ~ log2(e) * (dx - dx^2/2 + ... + (-1)^(n + 1) dx^n / n) +
+// dx^(n + 1) * P(dx)
+// such that the absolute error of P(dx) is smaller than double precision ulp,
+// then the overall relative error of our log2(1 + dx) approximation is:
+// ~ ulp(dx^(n + 1)) / |dx| ~ 2^(-52) * |dx^n|
+//
+// Let's consider 3 cases:
+// - For a fast path of a correctly rounded implementation, we will want the
+// relative error of x^y bounded above by:
+// ~ 2^(-precision - 10) = 2^(-53 - 10) = 2^(-63)
+// so that the Ziv accuracy test passes with probability > 1 - 2^-10 ~ 99.99%.
+// And from the above argument, we need to compute log2(1 + dx) with relative
+// error ~ 2^(-63 - 11) = 2^(-74).
+// So we will need to compute log2(1 + dx) accurately up to dx^n such that:
+// 2^(-52) * |dx^n| ~ 2^(-74).
+// Or equivalently:
+// |dx^n| < 2^(-22).
+// With our bounds from range reduction |dx| < 2^(-7), n = 3 will be quite
+// close to our desired precision needed.
+//
+// In summary, for a fast path of correctly rounded implementation, we will
+// compute log2(e) * (dx - dx^2/2 + dx^3/3) accurately, and approximate higher
+// terms with dx^4 * P(dx), where:
+// |(log2(1 + dx) - log2(e) * (dx - dx^2/2 + dx^3/3))/dx^4 - P(dx)| < 2^-53.
+//
+// - For a strictly < 1 ULP error fast version across the full exponent range
+// (as implemented in pow_fast), we will want the relative error of x^y
+// bounded above by:
+// ~ 2^(-precision) = 2^(-53).
+// So the relative error needed to approximate log2(1 + dx) is:
+// ~ 2^(-53 - 11) = 2^(-64),
+// and we will need to compute log2(1 + dx) accurately up to dx^n such that:
+// 2^(-52) * |dx^n| ~ 2^(-64),
+// or equivalently:
+// |dx^n| < 2^-12.
+// With |dx| < 2^(-7), n = 2 satisfies |dx^2| < 2^(-14) < 2^(-12).
+// Hence, we compute log2(e) * (dx - dx^2/2) accurately in DoubleDouble, and
+// approximate higher terms with dx^3 * P(dx), where:
+// |(log2(1 + dx) - log2(e) * (dx - dx^2/2))/dx^3 - P(dx)| < 2^-53.
+// This guarantees error < 1 ULP everywhere (empirically bounded by 0.75 ULP).
+//
+// - If we were to choose n = 1 (computing only log2(e) * dx accurately, and
+// approximating higher terms with dx^2 * P(dx)), the relative error of
+// log2(1 + dx) is only bounded by:
+// ~ 2^(-52) * |dx| < 2^(-52) * 2^(-7) = 2^(-59).
+// The resulting relative error of x^y is about:
+// ~ 2^(-59) * |e_h|.
+// While this achieves < 1 ULP for typical inputs with |e_h| <= 64, near
+// extreme exponent boundaries (|e_h| ~ 2^10 and |dx| ~ 2^-7) the error can
+// drift up to ~7.8 ULPs (bounded theoretically by ~16 ULPs).
LIBC_INLINE double pow(double x, double y) {
using namespace pow_internal;
using FPBits = fputil::FPBits<double>;
FPBits xbits(x), ybits(y);
-
- bool x_sign = xbits.sign() == Sign::NEG;
- bool y_sign = ybits.sign() == Sign::NEG;
-
- FPBits x_abs = xbits.abs();
- FPBits y_abs = ybits.abs();
-
- uint64_t x_mant = xbits.get_mantissa();
- uint64_t y_mant = ybits.get_mantissa();
uint64_t x_u = xbits.uintval();
- uint64_t x_a = x_abs.uintval();
- uint64_t y_a = y_abs.uintval();
+ uint64_t y_u = ybits.uintval();
+ uint64_t y_a = ybits.abs().uintval();
- double e_x = static_cast<double>(xbits.get_exponent());
- uint64_t sign = 0;
-
- ///////// BEGIN - Check exceptional cases ////////////////////////////////////
- // If x or y is signaling NaN
- if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
+ if (LIBC_UNLIKELY((x_u & 0x0003'FFFF'FFFF'FFFF) == 0 ||
+ (y_u & 0x000F'FFFF'FFFF'FFFF) == 0)) {
+ if (auto r = check_special_inputs(x, y); LIBC_UNLIKELY(r.has_value()))
+ return r.value();
}
- // The double precision number that is closest to 1 is (1 - 2^-53), which has
- // log2(1 - 2^-53) ~ -1.715...p-53.
- // So if |y| > |1075 / log2(1 - 2^-53)|, and x is finite:
- // |y * log2(x)| = 0 or > 1075.
- // Hence x^y will either overflow or underflow if x is not zero.
- if (LIBC_UNLIKELY(y_mant == 0 || y_a > 0x43d7'4910'd52d'3052 ||
- x_u == FPBits::one().uintval() ||
+ double e_x = static_cast<double>(xbits.get_exponent());
+ uint64_t x_mant = xbits.get_mantissa();
+ bool is_neg = false;
+ double sign_d = 1.0;
+
+ if (LIBC_UNLIKELY(y_a <= Y_LOWER_BOUND || y_a >= Y_UPPER_BOUND ||
x_u >= FPBits::inf().uintval() ||
x_u < FPBits::min_normal().uintval())) {
- // Exceptional exponents.
- if (y == 0.0)
- return 1.0;
-
- switch (y_a) {
- case 0x3fe0'0000'0000'0000: { // y = +-0.5
- // TODO: speed up x^(-1/2) with rsqrt(x) when available.
- if (LIBC_UNLIKELY(
- (x == 0.0 || x_u == FPBits::inf(Sign::NEG).uintval()))) {
- // pow(-0, 1/2) = +0
- // pow(-inf, 1/2) = +inf
- // Make sure it works correctly for FTZ/DAZ.
- return y_sign ? 1.0 / (x * x) : (x * x);
- }
- return y_sign ? (1.0 / fputil::sqrt<double>(x)) : fputil::sqrt<double>(x);
- }
- case 0x3ff0'0000'0000'0000: // y = +-1.0
- return y_sign ? (1.0 / x) : x;
- case 0x4000'0000'0000'0000: // y = +-2.0;
- return y_sign ? (1.0 / (x * x)) : (x * x);
- }
-
- // |y| > |1075 / log2(1 - 2^-53)|.
- if (y_a > 0x43d7'4910'd52d'3052) {
- if (y_a >= 0x7ff0'0000'0000'0000) {
- // y is inf or nan
- if (y_mant != 0) {
- // y is NaN
- // pow(1, NaN) = 1
- // pow(x, NaN) = NaN
- return (x_u == FPBits::one().uintval()) ? 1.0 : y;
- }
-
- // Now y is +-Inf
- if (x_abs.is_nan()) {
- // pow(NaN, +-Inf) = NaN
- return x;
- }
-
- if (x_a == 0x3ff0'0000'0000'0000) {
- // pow(+-1, +-Inf) = 1.0
- return 1.0;
- }
-
- if (x == 0.0 && y_sign) {
- // pow(+-0, -Inf) = +inf and raise FE_DIVBYZERO
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_DIVBYZERO);
- return FPBits::inf().get_val();
- }
- // pow (|x| < 1, -inf) = +inf
- // pow (|x| < 1, +inf) = 0.0
- // pow (|x| > 1, -inf) = 0.0
- // pow (|x| > 1, +inf) = +inf
- return ((x_a < FPBits::one().uintval()) == y_sign)
- ? FPBits::inf().get_val()
- : 0.0;
- }
- // x^y will overflow / underflow in double precision. Set y to a
- // large enough exponent but not too large, so that the computations
- // won't overflow in double precision.
- y = y_sign ? -0x1.0p100 : 0x1.0p100;
- }
-
- // y is finite and non-zero.
-
- if (x_u == FPBits::one().uintval()) {
- // pow(1, y) = 1
- return 1.0;
- }
-
- // TODO: Speed things up with pow(2, y) = exp2(y) and pow(10, y) = exp10(y).
-
- if (x == 0.0) {
- bool out_is_neg = x_sign && is_odd_integer(y);
- if (y_sign) {
- // pow(0, negative number) = inf
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_DIVBYZERO);
- return FPBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
- }
- // pow(0, positive number) = 0
- return out_is_neg ? -0.0 : 0.0;
- }
-
- if (x_a == FPBits::inf().uintval()) {
- bool out_is_neg = x_sign && is_odd_integer(y);
- if (y_sign)
- return out_is_neg ? -0.0 : 0.0;
- return FPBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
- }
-
- if (x_a > FPBits::inf().uintval()) {
- // x is NaN.
- // pow (aNaN, 0) is already taken care above.
- return x;
- }
-
- // Normalize denormal inputs.
- if (x_a < FPBits::min_normal().uintval()) {
- FPBits x_norm(x * 0x1.0p64);
- e_x = static_cast<double>(x_norm.get_exponent()) - 64.0;
- x_mant = x_norm.get_mantissa();
- }
-
- // x is finite and negative, and y is a finite integer.
- if (x_sign) {
- if (is_integer(y)) {
- x = -x;
- if (is_odd_integer(y))
- // sign = -1.0;
- sign = 0x8000'0000'0000'0000;
- } else {
- // pow( negative, non-integer ) = NaN
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
- }
- }
+ if (auto r = check_exceptional_cases(x, y, e_x, x_mant, is_neg, sign_d);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
}
- ///////// END - Check exceptional cases //////////////////////////////////////
-
// x^y = 2^( y * log2(x) )
// = 2^( y * ( e_x + log2(m_x) ) )
// First we compute log2(x) = e_x + log2(m_x)
@@ -382,66 +182,90 @@ LIBC_INLINE double pow(double x, double y) {
// log2(m_x) = log2( (1 + dx) / r )
// = log2(1 + dx) - log2(r).
- // In order for the overall computations x^y = 2^(y * log2(x)) to have the
- // relative errors < 2^-52 (1ULP), we will need to evaluate the exponent part
- // y * log2(x) with absolute errors < 2^-52 (or better, 2^-53). Since the
- // whole exponent range for double precision is bounded by
- // |y * log2(x)| < 1076 ~ 2^10, we need to evaluate log2(x) with absolute
- // errors < 2^-53 * 2^-10 = 2^-63.
+ // As analyzed in the overview comment above, we evaluate the cubic part
+ // (dx * C1 + dx^2 * C2 + dx^3 * C3) accurately using DoubleDouble (n = 3)
+ // and approximate higher terms with dx^4 * P(dx) to ensure relative error
+ // < 2^-63 for this fast path.
- // With that requirement, we use the following degree-6 polynomial
- // approximation:
- // P(dx) ~ log2(1 + dx) / dx
+ // Degree-5 polynomial approximation for:
+ // P(dx) ~ (log2(1 + dx) - (dx - dx^2/2 + dx^3/3)/log(2)) / dx^4
// Generated by Sollya with:
- // > P = fpminimax(log2(1 + x)/x, 6, [|D...|], [-2^-8, 2^-7]); P;
- // > dirtyinfnorm(log2(1 + x) - x*P, [-2^-8, 2^-7]);
- // 0x1.d03cc...p-66
- constexpr double COEFFS[] = {0x1.71547652b82fep0, -0x1.71547652b82e7p-1,
- 0x1.ec709dc3b1fd5p-2, -0x1.7154766124215p-2,
- 0x1.2776bd90259d8p-2, -0x1.ec586c6f3d311p-3,
- 0x1.9c4775eccf524p-3};
- // Error: ulp(dx^2) <= (2^-7)^2 * 2^-52 = 2^-66
- // Extra errors from various computations and rounding directions, the overall
- // errors we can be bounded by 2^-65.
-
- DoubleDouble dx_c0;
-
- // Perform exact range reduction and exact product dx * c0.
+ // > P = fpminimax((log2(1 + x) - (x - x^2/2 + x^3/3)/log(2))/x^4, 5,
+ // [|D...|], [-2^-8, 2^-7]);
+ // > dirtyinfnorm((log2(1 + x) - (x - x^2/2 + x^3/3)/log(2))/x - x^3*P,
+ // [-2^-8, 2^-7]);
+ // 0x1.a643c...p-74
+ // > dirtyinfnorm((log2(1 + x) - (x - x^2/2 + x^3/3)/log(2))/x^4 - P,
+ // [-2^-8, 2^-7]);
+ // 0x1.b81c5...p-53
+ constexpr double COEFFS[] = {-0x1.71547652b82fdp-2, 0x1.2776c50ef8f9bp-2,
+ -0x1.ec709dc4f0fedp-3, 0x1.a617677f716dep-3,
+ -0x1.715423d54d5p-3, 0x1.44e6355fc4d03p-3};
+
+ // Constants for C_k = (-1)^(k-1) / (k * log(2)) in DoubleDouble:
+ // C1 = 1 / log(2)
+ constexpr DoubleDouble C1 = {0x1.777d0ffda0d24p-56, 0x1.71547652b82fep0};
+ // C2 = -1 / (2 * log(2))
+ constexpr DoubleDouble C2 = {-0x1.777d0ffda0d24p-57, -0x1.71547652b82fep-1};
+ // C3 = 1 / (3 * log(2))
+ constexpr DoubleDouble C3 = {0x1.b749fc15522bcp-50, 0x1.ec709dc3a03e2p-2};
+
+ // Perform exact range reduction.
#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
double dx = fputil::multiply_add(RD[idx_x], m_x.get_val(), -1.0); // Exact
- dx_c0 = fputil::exact_mult(COEFFS[0], dx);
#else
double c = FPBits(m_x.uintval() & 0x3fff'e000'0000'0000).get_val();
double dx =
fputil::multiply_add(RD[idx_x], m_x.get_val() - c, CD[idx_x]); // Exact
- dx_c0 = fputil::exact_mult<double, 28>(dx, COEFFS[0]); // Exact
#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // Evaluate the cubic part (dx * C1 + dx^2 * C2 + dx^3 * C3) and polynomial
+ // tail using a parallel Double-Double Estrin scheme.
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // Error-free transformation for r = C1.hi + dx * C2.hi:
+ double r_hi = fputil::multiply_add(dx, C2.hi, C1.hi);
+ double r_lo = fputil::multiply_add(dx, C2.hi, C1.hi - r_hi); // Exact error
+#else
+ DoubleDouble dx_c2 = fputil::exact_mult(dx, C2.hi);
+ DoubleDouble r_sum = fputil::exact_add(C1.hi, dx_c2.hi);
+ double r_hi = r_sum.hi;
+ double r_lo = r_sum.lo + dx_c2.lo;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ // Low parts polynomial evaluated in parallel:
+ // C1.lo + dx * (C2.lo + dx * C3.lo) + r_lo
+ double lo_tail = fputil::multiply_add(dx, C3.lo, C2.lo);
+ double lo_poly = fputil::multiply_add(dx, lo_tail, C1.lo + r_lo);
+
+ // Evaluate polynomial tail P(dx) using Estrin's scheme:
double dx2 = dx * dx;
- double c0 = fputil::multiply_add(dx, COEFFS[2], COEFFS[1]);
- double c1 = fputil::multiply_add(dx, COEFFS[4], COEFFS[3]);
- double c2 = fputil::multiply_add(dx, COEFFS[6], COEFFS[5]);
+ double c0 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
+ double c1 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
+ double c2 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);
+
+ double dx4 = dx2 * dx2;
+ double d0 = fputil::multiply_add(dx2, c1, c0);
+ double p = fputil::multiply_add(dx4, c2, d0);
+
+ // High part of cubic term: C3.hi + dx * P(dx)
+ double q = fputil::multiply_add(dx, p, C3.hi);
- double p = fputil::polyeval(dx2, c0, c1, c2);
+ // Combine dx^2 * q into lo_poly:
+ lo_poly = fputil::multiply_add(dx2, q, lo_poly);
- // s = e_x - log2(r) + dx * P(dx)
- // Absolute error bound:
- // |log2(x) - log2_x.hi - log2_x.lo| < 2^-65.
+ // Multiply by dx to get log2(1 + dx) in DoubleDouble:
+ DoubleDouble log2_1p = fputil::exact_mult(dx, r_hi);
+ log2_1p.lo = fputil::multiply_add(dx, lo_poly, log2_1p.lo);
- // Notice that e_x - log2(r).hi is exact, so we perform an exact sum of
- // e_x - log2(r).hi and the high part of the product dx * c0:
- // log2_x_hi.hi + log2_x_hi.lo = e_x - log2(r).hi + (dx * c0).hi
+ // Combine with e_x - log2(r):
DoubleDouble log2_x_hi =
- fputil::exact_add(e_x + LOG2_R_DD[idx_x].hi, dx_c0.hi);
- // The low part is dx^2 * p + low part of (dx * c0) + low part of -log2(r).
- double log2_x_lo =
- fputil::multiply_add(dx2, p, dx_c0.lo + LOG2_R_DD[idx_x].lo);
- // Perform accurate sums.
+ fputil::exact_add(e_x + LOG2_R_DD[idx_x].hi, log2_1p.hi);
+ double log2_x_lo = log2_1p.lo + LOG2_R_DD[idx_x].lo;
DoubleDouble log2_x = fputil::exact_add(log2_x_hi.hi, log2_x_lo);
log2_x.lo += log2_x_hi.lo;
// To compute 2^(y * log2(x)), we break the exponent into 3 parts:
- // y * log(2) = hi + mid + lo, where
+ // y * log2(x) = hi + mid + lo, where
// hi is an integer
// mid * 2^6 is an integer
// |lo| <= 2^-7
@@ -451,96 +275,190 @@ LIBC_INLINE double pow(double x, double y) {
// and 2^lo ~ 1 + lo * P(lo).
// Thus, we have:
// hi + mid = 2^-6 * round( 2^6 * y * log2(x) )
- // If we restrict the output such that |hi| < 150, (hi + mid) uses (8 + 6)
+ // If we restrict the output such that |hi| < 512, (hi + mid) uses (9 + 6)
// bits, hence, if we use double precision to perform
// round( 2^6 * y * log2(x))
- // the lo part is bounded by 2^-7 + 2^(-(52 - 14)) = 2^-7 + 2^-38
+ // the lo part is bounded by 2^-7 + 2^(-(52 - 15)) = 2^-7 + 2^-37
// In the following computations:
// y6 = 2^6 * y
// hm = 2^6 * (hi + mid) = round(2^6 * y * log2(x)) ~ round(y6 * s)
// lo6 = 2^6 * lo = 2^6 * (y - (hi + mid)) = y6 * log2(x) - hm.
- double y6 = y * 0x1.0p6; // Exact.
+ constexpr double SCALE = 0x1.0p6;
+ double y6 = y * SCALE; // Exact.
DoubleDouble y6_log2_x = fputil::exact_mult(y6, log2_x.hi);
y6_log2_x.lo = fputil::multiply_add(y6, log2_x.lo, y6_log2_x.lo);
// Check overflow/underflow.
double scale = 1.0;
+ bool is_denorm = false;
// |2^(hi + mid) - exp2_hi_mid| <= ulp(exp2_hi_mid) / 2
- // Clamp the exponent part into smaller range that fits double precision.
- // For those exponents that are out of range, the final conversion will round
- // them correctly to inf/max float or 0/min float accordingly.
- constexpr double UPPER_EXP_BOUND = 512.0 * 0x1.0p6;
+
+ // The fast computation for 2^hi below requires that:
+ // |hi| < 512, or equivalently, |hm| < 512 * 2^6.
+ // This guarantees that the biased exponent:
+ // exp_biased = (hm_i >> 6) + EXP_BIAS = hi + 1023
+ // is strictly within the normal double range:
+ // 1023 - 511 <= exp_biased <= 1023 + 511, or 512 <= exp_biased <= 1534.
+ // Hence, 2^hi is always a normal, non-zero, finite power of 2, and the
+ // multiplication upper * exp2_hi is exact and never underflows or overflows.
+ //
+ // From the edge case checks above:
+ // 2^(-54) / 1074 <= |y| <= 1075 * 2^53, and
+ // |log_2(1 - 2^(-53))| <= |log2(x)| <= 1074.
+ // So their product is bounded by:
+ // 2^-117 < |y * log2(x)| < 2^74.
+ //
+ // The meaningful range for y * log2(x) in double precision is:
+ // -1075 <= y * log2(x) <= 1024.
+ // Any value > 1024 overflows, and any value < -1075 underflows.
+ //
+ // When |y * log2(x)| >= 511, we shift the exponent by an offset S:
+ // y * log2(x) = (y * log2(x) - S) + S
+ // and multiply by scale = 2^S at the end.
+ // To ensure the shifted exponent (y * log2(x) - S) stays within [-511, 511]:
+ // - For positive range [511, 1024]:
+ // 511 - S > -511 ==> S < 1022
+ // 1024 - S < 511 ==> S > 513
+ // - For negative range [-1076, -511]:
+ // -511 + S < 511 ==> S <= 1022
+ // -1076 + S > -511 ==> S >= 565
+ // Combined, the shift must satisfy: 565 <= S <= 1022.
+ // We choose S = 600, which yields:
+ // y * log2(x) - 600 in [-89, 424] for the positive range, and
+ // y * log2(x) + 600 in [-476, 89] for the negative range,
+ // both fitting well within [-511, 511].
+ //
+ // For exponents that are completely out of range:
+ // y * log2(x) > 1025 or y * log2(x) < -1076,
+ // the product can be as large as 2^74, so subtracting or adding 600 * 64
+ // would still overflow a 32-bit int when computing hm_i.
+ // We return early with overflow or underflow in these cases.
+ //
+ // Alternatively, for less branching or in SIMD / vector implementations, one
+ // could clamp y6_log2_x.hi to:
+ // - UPPER_EXP_BOUND (511 * 64) for overflow, which when multiplied by
+ // scale = 2^600 yields 2^1111 and correctly overflows.
+ // - -500 * 64 for underflow, which when multiplied by scale = 2^-600 yields
+ // 2^-1100 and correctly underflows.
+ // However, in this correctly rounded version, such clamping can cause the
+ // Ziv accuracy test to fail on the clamped exponent, unnecessarily
+ // triggering the slow accurate paths for overflow or underflow cases.
+
+ constexpr double UPPER_EXP_BOUND = 511.0 * SCALE;
if (LIBC_UNLIKELY(FPBits(y6_log2_x.hi).abs().get_val() >= UPPER_EXP_BOUND)) {
if (FPBits(y6_log2_x.hi).sign() == Sign::POS) {
- scale = 0x1.0p512;
- y6_log2_x.hi -= 512.0 * 64.0;
- if (y6_log2_x.hi > 513.0 * 64.0)
- y6_log2_x.hi = 513.0 * 64.0;
+ if (y6_log2_x.hi > 1025.0 * SCALE)
+ return set_overflow(is_neg);
+ scale = 0x1.0p600;
+ y6_log2_x.hi -= 600.0 * SCALE;
} else {
- scale = 0x1.0p-512;
- y6_log2_x.hi += 512.0 * 64.0;
- if (y6_log2_x.hi < (-1076.0 + 512.0) * 64.0)
- y6_log2_x.hi = -564.0 * 64.0;
+ if (y6_log2_x.hi <= -1021.0 * SCALE) {
+ if (y6_log2_x.hi < -1076.0 * SCALE)
+ return set_underflow(is_neg);
+ is_denorm = true;
+ } else {
+ scale = 0x1.0p-600;
+ y6_log2_x.hi += 600.0 * SCALE;
+ }
}
}
double hm = fputil::nearest_integer(y6_log2_x.hi);
// lo6 = 2^6 * lo.
- double lo6_hi = y6_log2_x.hi - hm;
- double lo6 = lo6_hi + y6_log2_x.lo;
+ DoubleDouble lo6 = fputil::exact_add(y6_log2_x.hi - hm, y6_log2_x.lo);
int hm_i = static_cast<int>(hm);
unsigned idx_y = static_cast<unsigned>(hm_i) & 0x3f;
// 2^hi
- int64_t exp2_hi_i = static_cast<int64_t>(
- static_cast<uint64_t>(static_cast<int64_t>(hm_i >> 6))
- << FPBits::FRACTION_LEN);
+ int hi = hm_i >> 6;
+ double exp2_hi = 1.0;
+ if (LIBC_LIKELY(!is_denorm)) {
+ int64_t exp2_hi_i = static_cast<int64_t>(
+ static_cast<uint64_t>(hi + FPBits::EXP_BIAS) << FPBits::FRACTION_LEN);
+ exp2_hi = FPBits(static_cast<uint64_t>(exp2_hi_i)).get_val();
+ }
+
// 2^mid
- int64_t exp2_mid_hi_i =
- static_cast<int64_t>(FPBits(EXP2_MID1[idx_y].hi).uintval());
- int64_t exp2_mid_lo_i =
- static_cast<int64_t>(FPBits(EXP2_MID1[idx_y].mid).uintval());
- // (-1)^sign * 2^hi * 2^mid
- // Error <= 2^hi * 2^-53
- uint64_t exp2_hm_hi_i =
- static_cast<uint64_t>(exp2_hi_i + exp2_mid_hi_i) + sign;
- // The low part could be 0.
- uint64_t exp2_hm_lo_i =
- idx_y != 0 ? static_cast<uint64_t>(exp2_hi_i + exp2_mid_lo_i) + sign
- : sign;
- double exp2_hm_hi = FPBits(exp2_hm_hi_i).get_val();
- double exp2_hm_lo = FPBits(exp2_hm_lo_i).get_val();
-
- // Degree-5 polynomial approximation P(lo6) ~ 2^(lo6 / 2^6) = 2^(lo).
- // Generated by Sollya with:
- // > P = fpminimax(2^(x/64), 5, [|1, D...|], [-2^-1, 2^-1]);
- // > dirtyinfnorm(2^(x/64) - P, [-0.5, 0.5]);
- // 0x1.a2b77e618f5c4c176fd11b7659016cde5de83cb72p-60
- constexpr double EXP2_COEFFS[] = {0x1p0,
- 0x1.62e42fefa39efp-7,
- 0x1.ebfbdff82a23ap-15,
- 0x1.c6b08d7076268p-23,
- 0x1.3b2ad33f8b48bp-31,
- 0x1.5d870c4d84445p-40};
-
- double lo6_sqr = lo6 * lo6;
-
- double d0 = fputil::multiply_add(lo6, EXP2_COEFFS[2], EXP2_COEFFS[1]);
- double d1 = fputil::multiply_add(lo6, EXP2_COEFFS[4], EXP2_COEFFS[3]);
- double pp = fputil::polyeval(lo6_sqr, d0, d1, EXP2_COEFFS[5]);
-
- double r = fputil::multiply_add(exp2_hm_hi * lo6, pp, exp2_hm_lo);
- r += exp2_hm_hi;
-
- return r * scale;
+ DoubleDouble exp2_mid{EXP2_MID1[idx_y].mid * sign_d,
+ EXP2_MID1[idx_y].hi * sign_d};
+
+ // Polynomial expansion for 2^(lo6/64):
+ // 2^(lo6/64) ~ 1 + lo6 * (log(2)/64) + lo6^2 * P(lo6)
+ // The linear term is computed in DoubleDouble, and the degree-3 polynomial
+ // P(lo6) is evaluated with standard double precision.
+ //
+ // hi and lo parts of log(2)/64, generated by Sollya with:
+ // > a = D(log(2)/64);
+ // > b = D(log(2)/64 - a);
+ constexpr DoubleDouble LOG_2_OVER_64 = {0x1.abc9e3b39803fp-62,
+ 0x1.62e42fefa39efp-7};
+
+ // Degree-3 polynomial approximation for (2^(lo6/64) - 1 - lo6*log(2)/64) /
+ // lo6^2: Generated by Sollya with:
+ // > f = (2^(x/64) - 1 - x*log(2)/64) / x^2;
+ // > P = fpminimax(f, 5, [|D...|], [-0.5, 0.5]);
+ // > dirtyinfnorm((f - P) * x^2, [-0.5, 0.5]);
+ // 0x1.1778...p-73
+ constexpr double EXP2_COEFFS[] = {
+ 0x1.ebfbdff82c58fp-15, 0x1.c6b08d704a0cp-23, 0x1.3b2ab6fb4d08fp-31,
+ 0x1.5d87fe77a735dp-40, 0x1.430d835610044p-49, 0x1.ffe67c38112c3p-59};
+
+ DoubleDouble lo_log_2 = fputil::quick_mult(lo6, LOG_2_OVER_64);
+ DoubleDouble lo_log_2_p1 = fputil::exact_add(1.0, lo_log_2.hi);
+ lo_log_2_p1.lo += lo_log_2.lo;
+
+ double lo6_sq = lo6.hi * lo6.hi;
+ double e0 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[1], EXP2_COEFFS[0]);
+ double e1 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[3], EXP2_COEFFS[2]);
+ double e2 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[5], EXP2_COEFFS[4]);
+
+ double lo6_4 = lo6_sq * lo6_sq;
+ double f0 = fputil::multiply_add(lo6_sq, e0, lo_log_2_p1.lo);
+ double f1 = fputil::multiply_add(lo6_sq, e2, e1);
+
+ lo_log_2_p1.lo = fputil::multiply_add(lo6_4, f1, f0);
+ DoubleDouble r = fputil::quick_mult(exp2_mid, lo_log_2_p1);
+
+ // Absolute error bound for r:
+ // - Base error from exp2 stage is bounded by 0x1.0p-64.
+ // - Propagated error from log2(x):
+ // |y| * err(log2(x)) * (2^mid * log(2)) <= |y| * 2^-73.5
+ // Dynamic error bound adapting to exponent scale |y|:
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double err_r =
+ fputil::multiply_add(FPBits(y).abs().get_val(), 0x1.8p-73, 0x1.0p-64);
+#else
+ // Without FMA, intermediate roundings increase log2(x) and exp2 errors.
+ double err_r =
+ fputil::multiply_add(FPBits(y).abs().get_val(), 0x1.cp-72, 0x1.0p-63);
+#endif
+
+ if (LIBC_UNLIKELY(is_denorm)) {
+ if (auto r_denorm = ziv_test_denorm(hi, r.hi, r.lo, err_r, is_neg);
+ LIBC_LIKELY(r_denorm.has_value())) {
+ return r_denorm.value();
+ }
+ return pow_accurate(x, y, is_neg, static_cast<int>(e_x), idx_x, dx);
+ }
+
+ double upper = r.hi + (r.lo + err_r);
+ double lower = r.hi + (r.lo - err_r);
+ if (LIBC_LIKELY(upper == lower)) {
+ double tmp1 = upper * exp2_hi;
+ return tmp1 * scale;
+ }
+
+ return pow_accurate(x, y, is_neg, static_cast<int>(e_x), idx_x, dx);
}
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_H
diff --git a/libc/src/__support/math/pow_accurate_128.h b/libc/src/__support/math/pow_accurate_128.h
new file mode 100644
index 0000000000000..36d03e3c09e42
--- /dev/null
+++ b/libc/src/__support/math/pow_accurate_128.h
@@ -0,0 +1,375 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// 128-bit accurate path for double-precision pow(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_128_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_128_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.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/common.h"
+#include "src/__support/frac128.h"
+#include "src/__support/frac64.h"
+#include "src/__support/integer_literals.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/common_constants.h"
+#include "src/__support/math/pow_accurate_256.h"
+#include "src/__support/math/pow_utils.h"
+#include "src/__support/uint128.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace pow_internal {
+
+using DFloat128 = typename fputil::DyadicFloat<128>;
+using MantissaType = typename DFloat128::MantissaType;
+using LIBC_NAMESPACE::operator""_u128;
+
+// exp2_f128 shares EXP2_MID_FRAC256 from pow_accurate_256.h via
+// Frac256::to_frac128().
+
+// Polynomial approximation for log2(1 + dx):
+// For dx in [-2^-8, 2^-7], we approximate:
+// P(x) ~ log2(1 + x) / (2 * x)
+// then:
+// log2(1 + dx) = 2 * dx * P(dx).
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > prec = 256;
+// > P = fpminimax(log2(1 + x) / (2 * x), 15, [|127...|], [-2^(-8), 2^(-7)],
+// fixed);
+// > for i from 0 to 15 do {
+// c = coeff(P, i);
+// c_int = round(abs(c) * 2^127, 256, RN);
+// lo = c_int mod 2^64;
+// hi = (c_int - lo) / 2^64;
+// print("Frac128({", lo, "ULL,", hi, "ULL}), // a_", i);
+// };
+// > dirtyinfnorm(log2(1 + x) - 2 * x * P, [-2^(-8), 2^(-7)]);
+// 0x1.c6524...p-139 < 2^-138.
+//
+// Expanding P(x) in Horner form:
+// P(x) = a_0 - a_1 * x + a_2 * x^2 - a_3 * x^3 + ...
+// We store the absolute values of the coefficients |a_k| as Frac128.
+//
+// With y = |dx|,
+// - if dx >= 0:
+// P(y) = |a_0| - y * (|a_1| - y * (|a_2| - ...))
+// - if dx < 0:
+// P(-y) = |a_0| + y * (|a_1| + y * (|a_2| + ...))
+LIBC_INLINE_VAR constexpr Frac128 LOG2_POLY_128[16] = {
+ Frac128({0xdf43ff68348e9f44ULL, 0x5c551d94ae0bf85dULL}), // a_0
+ Frac128({0xefa1ffb41a474fa2ULL, 0x2e2a8eca5705fc2eULL}), // a_1
+ Frac128({0x9fc15522bc34caf6ULL, 0x1ec709dc3a03fd74ULL}), // a_2
+ Frac128({0x77d0ffda0d2838eaULL, 0x171547652b82fe17ULL}), // a_3
+ Frac128({0x2ca733033277ada0ULL, 0x12776c50ef9bfe79ULL}), // a_4
+ Frac128({0x4fe0aa8003f351fbULL, 0x0f6384ee1d01febaULL}), // a_5
+ Frac128({0xb2490fe2b526029dULL, 0x0d30bb153d6f6c9fULL}), // a_6
+ Frac128({0xbc03c43080705d70ULL, 0x0b8aa3b295c17f0bULL}), // a_7
+ Frac128({0x8db7017adcd9fd4dULL, 0x0a42589ebe015466ULL}), // a_8
+ Frac128({0x7a1820b5f15e8892ULL, 0x093bb62877cdff26ULL}), // a_9
+ Frac128({0xfb43ea327ddb209cULL, 0x0864d424ca0b0f9dULL}), // a_10
+ Frac128({0xca32556a02086730ULL, 0x07b1c2770e8b0fabULL}), // a_11
+ Frac128({0xb05fb6d59979dddeULL, 0x071a3d579502ad57ULL}), // a_12
+ Frac128({0x4abbb99c78464ef1ULL, 0x06985d880590b72bULL}), // a_13
+ Frac128({0x86477a0bb2fb6379ULL, 0x06282cce144a837eULL}), // a_14
+ Frac128({0x25363223e4195930ULL, 0x05c5acde28e4f0bfULL}), // a_15
+};
+
+// Accurate log2(x) in 128-bit precision reusing range reduction from fast pass:
+// x = 2^x_e * m_x, with 1 <= m_x < 2
+// r = RD[idx_x]
+// dx = r * m_x - 1, with -2^-8 <= dx < 2^-7
+// log2(x) = x_e + (-log2(r)) + log2(1 + dx).
+// The value -log2(r) is provided by LOG2_RD_FRAC256[idx_x].
+// log2(1 + dx) is evaluated with a degree-15 minimax polynomial using a 2-tier
+// Horner scheme (Frac64 and Frac128).
+LIBC_INLINE DFloat128 log2_f128(int x_e, unsigned idx_x, double dx) {
+ using FPBits = fputil::FPBits<double>;
+
+ if (dx == 0.0) {
+ Frac128 log2_m = LOG2_RD_FRAC256[idx_x].to_frac128();
+ DFloat128 m_df(Sign::POS, -127, log2_m);
+ return fputil::quick_add(DFloat128(static_cast<double>(x_e)), m_df);
+ }
+
+ double abs_dx = (dx < 0.0) ? -dx : dx;
+ FPBits bits(abs_dx);
+ // Append hidden bit.
+ uint64_t mant = bits.get_mantissa() | (1ULL << 52);
+ int shift = 127 + bits.get_exponent() - 52;
+ Frac128 y =
+ (shift >= 0) ? Frac128((UInt<128>(mant) << shift).val) : Frac128(0);
+
+ // Evaluate log2(1 + dx) using Horner scheme in 2 stages:
+ // - Degree 8-15: 64-bit precision evaluation.
+ // - Degree 0-7: 128-bit precision evaluation.
+ //
+ // Since y = |dx| <= 2^-7, truncation errors at degree k is bounded by:
+ // y^k <= 2^(-7k).
+ //
+ // Step 1: 64-bit evaluation, with truncation error at degree 8 is bounded by:
+ // 2^-63 * y^8 <= 2^-119.
+ Frac64 y64 = y.to_frac64();
+ Frac64 p64 = LOG2_POLY_128[15].to_frac64();
+ if (dx >= 0.0) {
+ for (int k = 14; k >= 8; --k)
+ p64 = LOG2_POLY_128[k].to_frac64() - ((p64 * y64) << 1);
+ } else {
+ for (int k = 14; k >= 8; --k)
+ p64 = LOG2_POLY_128[k].to_frac64() + ((p64 * y64) << 1);
+ }
+
+ // Step 2: 128-bit evaluation.
+ Frac128 p({0, p64.val[0]});
+ if (dx >= 0.0) {
+ for (int k = 7; k >= 0; --k)
+ p = LOG2_POLY_128[k] - ((p * y) << 1);
+ } else {
+ for (int k = 7; k >= 0; --k)
+ p = LOG2_POLY_128[k] + ((p * y) << 1);
+ }
+
+ // log2(1 + dx) = 2 * y * P(y):
+ Frac128 log2_1p = (y * p) << 2;
+ Frac128 log2_rd = LOG2_RD_FRAC256[idx_x].to_frac128();
+ Frac128 log2_m = (dx >= 0.0) ? (log2_rd + log2_1p) : (log2_rd - log2_1p);
+
+ DFloat128 m_df(Sign::POS, -127, log2_m);
+ return fputil::quick_add(DFloat128(static_cast<double>(x_e)), m_df);
+}
+
+// Polynomial approximation for (2^x - 1) / x:
+// For x in [-2^-8, 2^-7], we approximate:
+// P(x) ~ (2^x - 1) / x
+// then:
+// 2^x = 1 + x * P(x).
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > prec = 500;
+// > P = fpminimax((2^x - 1)/x, 11, [|127...|], [-2^(-8), 2^(-7)], fixed);
+// > for i from 0 to 11 do {
+// c = coeff(P, i);
+// c_int = round(abs(c) * 2^127, 256, RN);
+// lo = c_int mod 2^64;
+// hi = (c_int - lo) / 2^64;
+// print("Frac128({", lo, "ULL,", hi, "ULL}), // c_", i);
+// };
+// > dirtyinfnorm(2^x - (1 + x * P), [-2^(-8), 2^(-7)]);
+// 0x1.4c648...p-139 < 2^-138.
+//
+// Expanding P(x) in Horner form:
+// P(x) = c_0 + c_1 * x + c_2 * x^2 + ...
+// where c_k ~ (log(2))^(k+1) / (k+1)! are all positive.
+//
+// With u = |x|,
+// - if x >= 0:
+// P(u) = c_0 + u * (c_1 + u * (c_2 + ...))
+// - if x < 0:
+// P(-u) = c_0 - u * (c_1 - u * (c_2 - ...))
+LIBC_INLINE_VAR constexpr Frac128 EXP2_POLY_128[12] = {
+ Frac128({0xe4f1d9cc01f97b58ULL, 0x58b90bfbe8e7bcd5ULL}), // c_0
+ Frac128({0x6f16b06ec9735fcbULL, 0x1ebfbdff82c58ea8ULL}), // c_1
+ Frac128({0xcce9d8aeccad816eULL, 0x071ac235c1282fe2ULL}), // c_2
+ Frac128({0x9ccbbe0b53ecfa40ULL, 0x013b2ab6fba4e772ULL}), // c_3
+ Frac128({0x20e2fed5a256d4e6ULL, 0x002bb0ffcf14ce62ULL}), // c_4
+ Frac128({0xdbd2c2a4af40b4b6ULL, 0x00050c244be1b1e1ULL}), // c_5
+ Frac128({0x1a199575c57ddfaaULL, 0x00007ff2ff1622c3ULL}), // c_6
+ Frac128({0x11fd7e22defbeaf5ULL, 0x00000b160111d2e4ULL}), // c_7
+ Frac128({0x7681395eb574dbfbULL, 0x000000da929e9cafULL}), // c_8
+ Frac128({0x1814285a9da99308ULL, 0x0000000f267a8ac6ULL}), // c_9
+ Frac128({0xafddba4f901557d4ULL, 0x00000000f46563d0ULL}), // c_10
+ Frac128({0x7ae715abc468c5f0ULL, 0x000000000e1de423ULL}), // c_11
+};
+
+// Compute 2^z in 128-bit precision:
+// Range reduction:
+// k = round(z * 64)
+// hi = k >> 6
+// idx = k & 0x3f
+// lo = z - k * 2^-6, with |lo| <= 2^-7.
+// Then:
+// 2^z = 2^hi * EXP2_MID_FRAC256[idx] * 2^lo
+// = 2^hi * EXP2_MID_FRAC256[idx] * (1 + lo * P(lo)).
+LIBC_INLINE DFloat128 exp2_f128(const DFloat128 &z) {
+ double z_d = static_cast<double>(z);
+ double z_scaled = z_d * 64.0;
+ double kd = fputil::nearest_integer(z_scaled);
+ int k = static_cast<int>(kd);
+
+ int hi = k >> 6;
+ unsigned idx = static_cast<unsigned>(k & 0x3f);
+
+ DFloat128 kd_f128(kd * 0x1.0p-6);
+ DFloat128 lo = fputil::quick_add(z, -kd_f128);
+
+ Frac128 m = EXP2_MID_FRAC256[idx].to_frac128();
+ if (LIBC_UNLIKELY(lo.mantissa.is_zero()))
+ return DFloat128(Sign::POS, hi - 127, m);
+
+ bool lo_is_neg = (lo.sign == Sign::NEG);
+ int shift = -127 - lo.exponent;
+ Frac128 u = (shift < 128) ? Frac128((lo.mantissa >> shift).val) : Frac128(0);
+
+ // Evaluate 2^lo - 1 = lo * P(lo) using Horner scheme in 2 stages:
+ // - Degree 7-11: 64-bit precision evaluation.
+ // - Degree 0-6: 128-bit precision evaluation.
+ //
+ // Since u = |lo| <= 2^-7, truncation errors at degree i is bounded by:
+ // u^i <= 2^(-7i).
+ //
+ // Step 1: 64-bit evaluation, with truncation error at degree 7 is bounded by:
+ // 2^-63 * u^7 <= 2^-112.
+ Frac64 u64 = u.to_frac64();
+ Frac64 p64 = EXP2_POLY_128[11].to_frac64();
+ if (lo_is_neg) {
+ for (int i = 10; i >= 7; --i)
+ p64 = EXP2_POLY_128[i].to_frac64() - ((p64 * u64) << 1);
+ } else {
+ for (int i = 10; i >= 7; --i)
+ p64 = EXP2_POLY_128[i].to_frac64() + ((p64 * u64) << 1);
+ }
+
+ // Step 2: 128-bit evaluation.
+ Frac128 p({0, p64.val[0]});
+ if (lo_is_neg) {
+ for (int i = 6; i >= 0; --i)
+ p = EXP2_POLY_128[i] - ((p * u) << 1);
+ } else {
+ for (int i = 6; i >= 0; --i)
+ p = EXP2_POLY_128[i] + ((p * u) << 1);
+ }
+
+ // Reconstruction:
+ // t = u * P(u) ~ 2^u - 1 as Frac128
+ // mt = m * t ~ m * (2^u - 1) as Frac128
+ // m_final = m +- mt ~ m * (1 +- t) ~ m * 2^lo
+ Frac128 t = (u * p) << 1;
+ Frac128 mt = (m * t) << 1;
+ Frac128 m_final = lo_is_neg ? (m - mt) : (m + mt);
+
+ if ((m_final.val[1] & (1ULL << 63)) == 0) {
+ m_final = m_final << 1;
+ --hi;
+ }
+
+ return DFloat128(Sign::POS, hi - 127, m_final);
+}
+
+// Accurate pow(x, y) reusing range reduction parameters from the fast pass.
+LIBC_INLINE double pow_accurate(double x, double y, bool is_neg, int x_e,
+ unsigned idx_x, double dx) {
+ DFloat128 log2_x = log2_f128(x_e, idx_x, dx);
+ DFloat128 y_f128(y);
+ DFloat128 z = fputil::quick_mul(y_f128, log2_x);
+
+ // 2^1025 > max normal double.
+ double z_d = static_cast<double>(z);
+ if (LIBC_UNLIKELY(z_d >= 1025.0))
+ return set_overflow(is_neg);
+
+ // 2^-1076 < min subnormal double.
+ if (LIBC_UNLIKELY(z_d <= -1076.0))
+ return set_underflow(is_neg);
+
+ // For 0 < |z| <= 2^-55, x^y is between 1 - 2^-54 and 1 + 2^-53.
+ if (LIBC_UNLIKELY(!z.mantissa.is_zero() && z.exponent + 127 <= -55)) {
+ volatile double one = 1.0;
+ volatile double eps = (z.sign == Sign::NEG) ? -0x1.0p-100 : 0x1.0p-100;
+ double res = one + eps;
+ return is_neg ? -res : res;
+ }
+
+ DFloat128 r = exp2_f128(z);
+ if (is_neg)
+ r.sign = Sign::NEG;
+
+ int unbiased_exp = r.exponent + 127;
+ if (LIBC_UNLIKELY(unbiased_exp >= 1024))
+ return set_overflow(is_neg);
+
+ // Check if r is close to a 54-bit rounding boundary (either an exact 53-bit
+ // float or a midpoint).
+ // For normal numbers, the 54 leading bits occupy bits [127:74] of r.mantissa,
+ // leaving shift = 128 - 54 = 74 fractional bits below the boundary.
+ // For subnormal numbers (unbiased_exp < -1022), the boundary shifts right.
+ int shift = 74;
+ if (LIBC_UNLIKELY(unbiased_exp < -1022))
+ shift = 74 - (-1022 - unbiased_exp);
+
+ bool is_boundary_candidate = false;
+ // Lauter & Lefevre (2009) showed that if x^y is not an exact 54-bit number,
+ // the distance to the nearest 54-bit boundary is at least:
+ // |x^y - o_54(x^y)| / x^y >= 2^-114.
+ // In r.mantissa, this minimum distance is 2^-114 * 2^127 = 2^13.
+ //
+ // Since eps_bound >= 2^13, when shift < 14, the maximum possible distance to
+ // the nearest boundary, 2^(shift - 1) <= 2^12, is strictly less than
+ // eps_bound, making every value a candidate. Checking shift >= 14 also
+ // prevents undefined negative shifts for denormals.
+ if (shift >= 14) {
+ // Distance from r.mantissa to the nearest multiple of 2^shift:
+ MantissaType mask = (MantissaType(1) << shift) - 1;
+ MantissaType rem = r.mantissa & mask;
+ MantissaType half = MantissaType(1) << (shift - 1);
+ MantissaType dist =
+ (rem <= half) ? rem : ((MantissaType(1) << shift) - rem);
+
+ // Evaluation error in r.mantissa is bounded by:
+ // log(2) * |y| * AbsErr(log2_x) * 2^127 + AbsErr(exp2_z)
+ // <= log(2) * 32 * |y| + 2^13
+ // < 32 * |y| + 2^13.
+ // For |y| >= 2^59, 32 * |y| >= 2^64, so we cap eps_bound at 2^64 to avoid
+ // 64-bit integer overflow.
+ double abs_y = (y < 0.0) ? -y : y;
+ MantissaType eps_bound = MantissaType(1) << 13;
+ if (LIBC_LIKELY(abs_y < 0x1.0p59))
+ eps_bound += (MantissaType(static_cast<uint64_t>(abs_y)) << 5);
+ else
+ eps_bound = (MantissaType(1) << 64);
+ is_boundary_candidate = (dist < eps_bound);
+ } else {
+ is_boundary_candidate = true;
+ }
+
+ if (LIBC_UNLIKELY(is_boundary_candidate)) {
+ uint64_t exact_m = 0;
+ int exact_exp = 0;
+ if (is_exact_rounding_boundary(x, y, exact_m, exact_exp)) {
+ int l = 64 - cpp::countl_zero(exact_m);
+ DFloat128 exact_f128(r.sign, exact_exp + l - 128,
+ MantissaType(exact_m) << (128 - l));
+ exact_f128.normalize();
+
+ return static_cast<double>(exact_f128);
+ }
+
+ return pow_accurate_256(y, is_neg, x_e, idx_x, dx);
+ }
+
+ double res = static_cast<double>(r);
+ if (LIBC_UNLIKELY(fputil::FPBits<double>(res).is_inf()))
+ return set_overflow(is_neg);
+ return res;
+}
+
+} // namespace pow_internal
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_128_H
diff --git a/libc/src/__support/math/pow_accurate_256.h b/libc/src/__support/math/pow_accurate_256.h
new file mode 100644
index 0000000000000..bb2deae4c4fe8
--- /dev/null
+++ b/libc/src/__support/math/pow_accurate_256.h
@@ -0,0 +1,818 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// 256-bit accurate path for double-precision pow(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_256_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_256_H
+
+#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/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/common.h"
+#include "src/__support/frac256.h"
+#include "src/__support/frac64.h"
+#include "src/__support/integer_literals.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/pow_utils.h"
+#include "src/__support/uint128.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace pow_internal {
+
+using DFloat256 = typename fputil::DyadicFloat<256>;
+using Mantissa256 = typename DFloat256::MantissaType;
+
+// Polynomial approximation for log2(1 + dx):
+// For dx in [-2^-8, 2^-7], we approximate:
+// P(x) ~ log2(1 + x) / (2 * x)
+// then:
+// log2(1 + dx) = 2 * dx * P(dx).
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > prec = 512;
+// > P = fpminimax(log2(1 + x) / (2 * x), 30, [|255...|], [-2^(-8), 2^(-7)],
+// fixed);
+// > for i from 0 to 30 do {
+// c = coeff(P, i);
+// c_int = round(abs(c) * 2^255, 512, RN);
+// print(c_int);
+// };
+// > dirtyinfnorm(log2(1 + x) - 2 * x * P, [-2^(-8), 2^(-7)]);
+// 0x1.bb204...p-270 < 2^-266.
+//
+// Expanding P(x) in Horner form:
+// P(x) = a_0 - a_1 * x + a_2 * x^2 - a_3 * x^3 + ...
+// We store the absolute values of the coefficients |a_k| as Frac256.
+//
+// With y = |dx|,
+// - if dx >= 0:
+// P(y) = |a_0| - y * (|a_1| - y * (|a_2| - ...))
+// - if dx < 0:
+// P(-y) = |a_0| + y * (|a_1| + y * (|a_2| + ...))
+LIBC_INLINE_VAR constexpr Frac256 LOG2_POLY_256[31] = {
+ Frac256({0x45928b3668d09924ULL, 0x75abbd546eb4ad2cULL,
+ 0xdf43ff68348e9f44ULL, 0x5c551d94ae0bf85dULL}), // a_0
+ Frac256({0x22c9459b34684c8aULL, 0x3ad5deaa375a5696ULL,
+ 0xefa1ffb41a474fa2ULL, 0x2e2a8eca5705fc2eULL}), // a_1
+ Frac256({0x1730d91222eb39a3ULL, 0x27393f1c24e6e464ULL,
+ 0x9fc15522bc2f8a6cULL, 0x1ec709dc3a03fd74ULL}), // a_2
+ Frac256({0x1164a2cd9b640d52ULL, 0x1d6aef551bad2b4bULL,
+ 0x77d0ffda0d23a7d1ULL, 0x171547652b82fe17ULL}), // a_3
+ Frac256({0xa783b5e0b932720eULL, 0x1788bf77495755d5ULL,
+ 0x2ca73314d74fb974ULL, 0x12776c50ef9bfe79ULL}), // a_4
+ Frac256({0x0b986b9531f9d509ULL, 0x139c9f8e12737232ULL,
+ 0x4fe0aa915e17c536ULL, 0x0f6384ee1d01febaULL}), // a_5
+ Frac256({0x97a9d262490024e5ULL, 0x10cf6430a219cf98ULL,
+ 0xb22e490ee2efcd9cULL, 0x0d30bb153d6f6c9fULL}), // a_6
+ Frac256({0x8a48b0b4c37bacd5ULL, 0x8eb577aa8dd695a4ULL,
+ 0xbbe87fed0691d3e8ULL, 0x0b8aa3b295c17f0bULL}), // a_7
+ Frac256({0xb0583590fcb47a20ULL, 0x62686a5eb6f7bc5bULL,
+ 0x354071b63eba8379ULL, 0x0a42589ebe01547cULL}), // a_8
+ Frac256({0xed3f5ebc7c98a745ULL, 0x0bc45fbba4b78402ULL,
+ 0x9653998a6ba7dcbaULL, 0x093bb62877cdff3cULL}), // a_9
+ Frac256({0x7359984444b13a62ULL, 0x96556e4d1cd07a22ULL,
+ 0x2b91d166906a0e7aULL, 0x0864d424ca011694ULL}), // a_10
+ Frac256({0x33fa41c7bb5906d1ULL, 0x09ce4f8602a57acbULL,
+ 0x27f05548af0be29bULL, 0x07b1c2770e80ff5dULL}), // a_11
+ Frac256({0x3cd309be8234c315ULL, 0xf55cdd54ec97ae7fULL,
+ 0x388f13a58de39618ULL, 0x071a3d5a34c5d807ULL}), // a_12
+ Frac256({0xa481c98db0b1aff6ULL, 0x093722cc884edbbfULL,
+ 0xd91724877177e6ceULL, 0x06985d8a9eb7b64fULL}), // a_13
+ Frac256({0xc5ebfbd7f8f11467ULL, 0xba98bfd7dd3b3734ULL,
+ 0xb98d1106f26fe87aULL, 0x0627cec5a533ff7dULL}), // a_14
+ Frac256({0x8f5b01e8d6775018ULL, 0x761efcd442a1522bULL,
+ 0xddf43ff68348e853ULL, 0x05c551d94ae0bf85ULL}), // a_15
+ Frac256({0xaa2b7924b8634b56ULL, 0x0bb2cd446fa9eca9ULL,
+ 0xd0e5e1d8f4097fddULL, 0x056e6b26dd0fc350ULL}), // a_16
+ Frac256({0x5a3861ff74fe9726ULL, 0x67655566edd6ee5dULL,
+ 0x1aa038db217aa2f2ULL, 0x05212c4f5f00aa3eULL}), // a_17
+ Frac256({0x8e5b3868f6fd9a2bULL, 0x86ad9ad6b3808b9dULL,
+ 0x26b2bc9961df0451ULL, 0x04dc0f07d343ff99ULL}), // a_18
+ Frac256({0xd5aff98461aaaa10ULL, 0x6d86f9b7f2a0e47bULL,
+ 0x4b29cb04bafb62c9ULL, 0x049ddb143be6ff9eULL}), // a_19
+ Frac256({0x73e005b24f0b5921ULL, 0x79754c39f3650d18ULL,
+ 0x3b6462eeb984c6bfULL, 0x046593b1bf252435ULL}), // a_20
+ Frac256({0xbfe649c5ee1a2476ULL, 0xa979d5e9dd8792fcULL,
+ 0x16a5af95b3d52886ULL, 0x04326a1265008b4aULL}), // a_21
+ Frac256({0xb240c3c658e84513ULL, 0x5a0dc3a6d6dcfbcbULL,
+ 0x3250ee4b9aa7a10cULL, 0x0403b35f8200853cULL}), // a_22
+ Frac256({0xa79c9e1b001812d1ULL, 0xc371e02ecf7b2649ULL,
+ 0xbc02b8b45c8cbd91ULL, 0x03d8e13b87407f7cULL}), // a_23
+ Frac256({0x396cbdb266e04da2ULL, 0xf42521f706b70fb3ULL,
+ 0xc4949cec1b5b8eb0ULL, 0x03b17c102febc7f2ULL}), // a_24
+ Frac256({0x623e4d738fe3f137ULL, 0x6d666a8379a8770dULL,
+ 0xddbb37796ef31775ULL, 0x038d1ead1a5f0977ULL}), // a_25
+ Frac256({0x90aab3e1d0b7c642ULL, 0x65353252593e6d5bULL,
+ 0x0a50523f5de27e92ULL, 0x036b72dfa00d56f3ULL}), // a_26
+ Frac256({0x57b2859b2032fbedULL, 0xe13ac78494bc3804ULL,
+ 0xc717595704ee1587ULL, 0x034c2ec97b40b274ULL}), // a_27
+ Frac256({0xaacab0e905a586e3ULL, 0x4c5437b1d3769115ULL,
+ 0x3db148838720c9e9ULL, 0x032f124ccfb211fcULL}), // a_28
+ Frac256({0xb0f6f05d9338d95fULL, 0x43c399e43c54795fULL,
+ 0x44b5111834acb3d6ULL, 0x031357dab283c21cULL}), // a_29
+ Frac256({0x676b154931b3d988ULL, 0x8e0bff4ea0623d39ULL,
+ 0xe8a2fa21caed0f90ULL, 0x02d819fe5a43f238ULL}), // a_30
+};
+
+// Polynomial approximation for (2^x - 1) / x:
+// For x in [-2^-8, 2^-7], we approximate:
+// P(x) ~ (2^x - 1) / x
+// then:
+// 2^x = 1 + x * P(x).
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > prec = 512;
+// > P = fpminimax((2^x - 1) / x, 21, [|255...|], [-2^(-8), 2^(-7)], fixed);
+// > for i from 0 to 21 do {
+// c = coeff(P, i);
+// c_int = round(abs(c) * 2^255, 512, RN);
+// print(c_int);
+// };
+// > dirtyinfnorm(2^x - (1 + x * P), [-2^(-8), 2^(-7)]);
+// 0x1.1554a...p-269 < 2^-268.
+//
+// Expanding P(x) in Horner form:
+// P(x) = c_0 + c_1 * x + c_2 * x^2 + ...
+// where c_k ~ (log(2))^(k+1) / (k+1)! are all positive.
+//
+// With u = |x|,
+// - if x >= 0:
+// P(u) = c_0 + u * (c_1 + u * (c_2 + ...))
+// - if x < 0:
+// P(-u) = c_0 - u * (c_1 - u * (c_2 - ...))
+LIBC_INLINE_VAR constexpr Frac256 EXP2_POLY_256[22] = {
+ Frac256({0xc5068badc5d57d16ULL, 0xa079a193394c5b16ULL,
+ 0xe4f1d9cc01f97b57ULL, 0x58b90bfbe8e7bcd5ULL}), // c_0
+ Frac256({0xc2be93bbb1396e6eULL, 0xa3a2751c30ce69d4ULL,
+ 0x6f16b06ec9735fcaULL, 0x1ebfbdff82c58ea8ULL}), // c_1
+ Frac256({0x4753198ade221489ULL, 0xa7ae23a226d00887ULL,
+ 0xcce9d8aeccaf4b7bULL, 0x071ac235c1282fe2ULL}), // c_2
+ Frac256({0x699b709699e59d1cULL, 0x72478ea53e63911dULL,
+ 0x9ccbbe0b53eeac50ULL, 0x013b2ab6fba4e772ULL}), // c_3
+ Frac256({0xc8cdc3747ffed019ULL, 0x60aed94d2dce32e1ULL,
+ 0x20e2fed34a297d86ULL, 0x002bb0ffcf14ce62ULL}), // c_4
+ Frac256({0x549592ed6dd93585ULL, 0xcfb314ffccc47bb0ULL,
+ 0xdbd2c2a261ac8d07ULL, 0x00050c244be1b1e1ULL}), // c_5
+ Frac256({0x5d0fe48ff2c6f3f5ULL, 0x4ace1152e8810feeULL,
+ 0x1a1ac547321f639aULL, 0x00007ff2ff1622c3ULL}), // c_6
+ Frac256({0x7d81cc760a6cb8d0ULL, 0x2586e1a0f7107ab9ULL,
+ 0x11fec7ff3036d3beULL, 0x00000b160111d2e4ULL}), // c_7
+ Frac256({0x08c1c3fcab2a43d5ULL, 0x84518cb8caa9ba63ULL,
+ 0x3e1ed253872d27fcULL, 0x000000da929e9cafULL}), // c_8
+ Frac256({0x4cefc19f553e7564ULL, 0xb4ce2a0608a16541ULL,
+ 0xc764fb7ed0eca973ULL, 0x0000000f267a8ac5ULL}), // c_9
+ Frac256({0x561e9aeb9041ed25ULL, 0x90a4edc7b3b041c3ULL,
+ 0x8dd92607abccaf23ULL, 0x00000000f465639aULL}), // c_10
+ Frac256({0x277d61ed12e707bcULL, 0xd17730e76ad24fc6ULL,
+ 0x7e14c2f15ab43f0bULL, 0x000000000e1deb28ULL}), // c_11
+ Frac256({0xd72a83926cf51bcaULL, 0x7b8d8d1512799500ULL,
+ 0x8b3687cb140d6180ULL, 0x0000000000c0b0c9ULL}), // c_12
+ Frac256({0xd3a3f9891ddac9deULL, 0x2ae5760179b7f299ULL,
+ 0x26ac3c54b9f8a1b1ULL, 0x0000000000098a4bULL}), // c_13
+ Frac256({0xa41178da4233dfdbULL, 0xfc55ac7cb3287dbcULL,
+ 0xa10ec1008799ec55ULL, 0x00000000000070dbULL}), // c_14
+ Frac256({0xbd533bacbc8f3ab8ULL, 0x93b765826d13ae70ULL,
+ 0xa26b9e7e2ce48e3fULL, 0x00000000000004e3ULL}), // c_15
+ Frac256({0xe33450363b8c0d3bULL, 0xc09a8ea3f919a324ULL,
+ 0x088968384b4faac9ULL, 0x0000000000000033ULL}), // c_16
+ Frac256({0x2bf8d71ab53ade5eULL, 0x48f9567d04173c92ULL,
+ 0xf7176bdb43695d72ULL, 0x0000000000000001ULL}), // c_17
+ Frac256({0x52f008e29e649d62ULL, 0x79d54d84b16ad619ULL,
+ 0x125a7ecb835b4a46ULL, 0x0000000000000000ULL}), // c_18
+ Frac256({0x3e1c2787b4d48f1bULL, 0x3588e474c82d686bULL,
+ 0x00a2d6625a829f19ULL, 0x0000000000000000ULL}), // c_19
+ Frac256({0x5e8c8f1088f16459ULL, 0xc735fd5e15f67f2eULL,
+ 0x00055ff17b41b4b0ULL, 0x0000000000000000ULL}), // c_20
+ Frac256({0xfd5b04031fd9182dULL, 0x8d6a9e1233c035ddULL,
+ 0x00002b59f51f6ec7ULL, 0x0000000000000000ULL}), // c_21
+};
+
+// Lookup table for 2^(i / 64) for i = 0..63 as Frac256.
+// Bit 255 represents 2^0 = 1, bits 254..0 represent the fractional part.
+// Generated by Sollya with:
+// > prec = 512;
+// > for i from 0 to 63 do {
+// v = round(2^(i / 64) * 2^255, 256, RN);
+// print(v);
+// };
+LIBC_INLINE_VAR constexpr Frac256 EXP2_MID_FRAC256[64] = {
+ Frac256({0x0000000000000000ULL, 0x0000000000000000ULL,
+ 0x0000000000000000ULL, 0x8000000000000000ULL}),
+ Frac256({0xd08075ac1f200e4cULL, 0x9eb851655e2e5c4dULL,
+ 0x7be56527bd14def4ULL, 0x8164d1f3bc030773ULL}),
+ Frac256({0x2502f15067378a17ULL, 0x29f1a4afbefa5d7cULL,
+ 0x3e2a475b46520bffULL, 0x82cd8698ac2ba1d7ULL}),
+ Frac256({0x806bddad09d9c4a3ULL, 0x0d96b414ec4c9d06ULL,
+ 0x1af92eca13fd1582ULL, 0x843a28c3acde4046ULL}),
+ Frac256({0x5d42b362af1ee859ULL, 0x148a0459e7585151ULL,
+ 0xc5c95b8c2154c1b2ULL, 0x85aac367cc487b14ULL}),
+ Frac256({0x5229a7352c9b247bULL, 0x259ac58894f4fcb3ULL,
+ 0x3a1727c57b52a956ULL, 0x871f61969e8d1010ULL}),
+ Frac256({0x8bc3587fb118c94dULL, 0xe623d58b3772ba13ULL,
+ 0x5df8d76c98c67562ULL, 0x88980e8092da8527ULL}),
+ Frac256({0xe9c32d22e935007dULL, 0x259c4df53d76e910ULL,
+ 0x080ca1d92c3680c2ULL, 0x8a14d575496efd9aULL}),
+ Frac256({0x91e135ee84a3f734ULL, 0x1aa84ffbebac349fULL,
+ 0xfbe4628758a53c90ULL, 0x8b95c1e3ea8bd6e6ULL}),
+ Frac256({0x724a166325437476ULL, 0x183926ae7d718dc2ULL,
+ 0xb4c7b4968e41ad36ULL, 0x8d1adf5b7e5ba9e5ULL}),
+ Frac256({0xeb90ce3700bf59b6ULL, 0xa11037230b367828ULL,
+ 0x2dc0144c8783d4c5ULL, 0x8ea4398b45cd53c0ULL}),
+ Frac256({0x6f398dfe3f7903f1ULL, 0x43e90e15c2002132ULL,
+ 0x775814a8494e87e2ULL, 0x9031dc431466b1dcULL}),
+ Frac256({0xf1203caf65bfb9b9ULL, 0x1942b34816fb4f26ULL,
+ 0x0fd6d8e0ae5ac9d8ULL, 0x91c3d373ab11c336ULL}),
+ Frac256({0x583eab6852a22bb1ULL, 0x2748c36eeaffa273ULL,
+ 0xd339940e9d924ee7ULL, 0x935a2b2f13e6e92bULL}),
+ Frac256({0x035fb634c2e63a0fULL, 0x4856046901ff6c05ULL,
+ 0x2e8afad12551de54ULL, 0x94f4efa8fef70961ULL}),
+ Frac256({0x9a22b1526bb6a2e4ULL, 0xe0e68d9f200c5358ULL,
+ 0x48ea9b683a9c22c4ULL, 0x96942d3720185a00ULL}),
+ Frac256({0xd78b65cbefa7bb70ULL, 0x5e139a1b14fa8178ULL,
+ 0x46ad23182e42f6f6ULL, 0x9837f0518db8a96fULL}),
+ Frac256({0x1560e51a5df911dcULL, 0x8ac981ca9ceca6b3ULL,
+ 0xe43086cb34b5fcaeULL, 0x99e0459320b7fa64ULL}),
+ Frac256({0x9769d9b0a908a786ULL, 0x0928b5fce34cdf21ULL,
+ 0xa2a817a2a3cc3f1fULL, 0x9b8d39b9d54e5538ULL}),
+ Frac256({0x33a6fe2d4fd53e8aULL, 0x1ff17c29677589a0ULL,
+ 0xde494cf050e99b0bULL, 0x9d3ed9a72cffb750ULL}),
+ Frac256({0x21f977fe7c7fa118ULL, 0x65c15c122133e2a2ULL,
+ 0xa0911f09ebb9fdd1ULL, 0x9ef5326091a111adULL}),
+ Frac256({0x9f33f7bc78dc629fULL, 0x782a0735d02b1a20ULL,
+ 0x192dc79edb0fd9a9ULL, 0xa0b0510fb9714fc2ULL}),
+ Frac256({0x5a7a799221808de9ULL, 0x9da4384dbc2c8eaeULL,
+ 0x9b7a04ef80cfdea7ULL, 0xa27043030c496818ULL}),
+ Frac256({0x4c72418596cc5bd0ULL, 0xbae743abfbc07376ULL,
+ 0x0d1db4831781e1eeULL, 0xa43515ae09e6809eULL}),
+ Frac256({0x2589c98a8290d3f0ULL, 0x1dd170ace2bcfc17ULL,
+ 0x1cbd7f621710701bULL, 0xa5fed6a9b15138eaULL}),
+ Frac256({0xdd30939a1d1e929cULL, 0x01424bd194d3999eULL,
+ 0x9ec5b4d5039f72afULL, 0xa7cd93b4e9653569ULL}),
+ Frac256({0x325c9e2203504517ULL, 0x3951f214c02d824aULL,
+ 0x541e24ec3531fa73ULL, 0xa9a15ab4ea7c0ef8ULL}),
+ Frac256({0x967357d6b36df9f8ULL, 0x7ad59ec00ebe6393ULL,
+ 0x658023b2759e0079ULL, 0xab7a39b5a93ed337ULL}),
+ Frac256({0xb165f141833a67daULL, 0x6be409407034fdedULL,
+ 0x4980a8c8f59a2ec4ULL, 0xad583eea42a14ac6ULL}),
+ Frac256({0x5a8c73beaa946990ULL, 0xa4502c14f429ded9ULL,
+ 0xdf26101ccbb35032ULL, 0xaf3b78ad690a4374ULL}),
+ Frac256({0x97ced890d5b0b0c0ULL, 0x757cfb9913adc577ULL,
+ 0x87d037e96d215d8eULL, 0xb123f581d2ac258fULL}),
+ Frac256({0xba1e54cf684354dfULL, 0xfa6e051d6f8bc3ffULL,
+ 0x3ecf14dc798a519bULL, 0xb311c412a9112489ULL}),
+ Frac256({0xed17ac8583339915ULL, 0x1d6f60ba893ba84cULL,
+ 0x597d89b3754abe9fULL, 0xb504f333f9de6484ULL}),
+ Frac256({0x20850e774a86cd8fULL, 0xf88abbe777df360eULL,
+ 0x07165f0ddd541a59ULL, 0xb6fd91e328d17791ULL}),
+ Frac256({0x322d7893ed4da9a8ULL, 0xa5ab16cf451056edULL,
+ 0x1b879778566b65a1ULL, 0xb8fbaf4762fb9ee9ULL}),
+ Frac256({0x6c373a75c2828202ULL, 0x02f30d0bdcaa516dULL,
+ 0x74d519d24593838cULL, 0xbaff5ab2133e45fbULL}),
+ Frac256({0x0d9a4be023ece032ULL, 0x15b34bbcb0298f41ULL,
+ 0xa8811fb66d0faf7aULL, 0xbd08a39f580c36beULL}),
+ Frac256({0x83ea957596be426dULL, 0xa13fc7e6faf9c830ULL,
+ 0xe815d0abcbf0b850ULL, 0xbf1799b67a731082ULL}),
+ Frac256({0xdefefee72ae7a33dULL, 0x6b2e5dd607a9969cULL,
+ 0x7c457d59a50087b5ULL, 0xc12c4cca66709456ULL}),
+ Frac256({0x5b718d616c4fef19ULL, 0x6b9f89b7dabbcb2bULL,
+ 0x20ec856128b83a42ULL, 0xc346ccda24976407ULL}),
+ Frac256({0xc7686006e4e6c093ULL, 0x6b0f939998251a36ULL,
+ 0x3e2ad0c964dd9f37ULL, 0xc5672a115506daddULL}),
+ Frac256({0xcea65224bc9900d0ULL, 0x4da570a2c574a304ULL,
+ 0xc13a2e3976c0277eULL, 0xc78d74c8abb9b15cULL}),
+ Frac256({0xf4dd023ff93c7ffbULL, 0x257ac0db1f419377ULL,
+ 0x80e1f92a0511697eULL, 0xc9b9bd866e2f27a2ULL}),
+ Frac256({0x639aa6f940962626ULL, 0xeb8a25b7b40c0426ULL,
+ 0xf4907c8f45ebf6dcULL, 0xcbec14fef2727c5cULL}),
+ Frac256({0x2bbd398af35c079fULL, 0x6f28610b8c36485aULL,
+ 0xe235838f95f2c6edULL, 0xce248c151f8480e3ULL}),
+ Frac256({0x2a33269ab05c3e5dULL, 0x11546d3ea28976d6ULL,
+ 0xd6d45c6559a4d502ULL, 0xd06333daef2b2594ULL}),
+ Frac256({0xfa7663033f05357bULL, 0x52029c0b81f7be57ULL,
+ 0x12248e57c3de4028ULL, 0xd2a81d91f12ae45aULL}),
+ Frac256({0xfa628009459a2417ULL, 0xb8e7a32e5783da5cULL,
+ 0x5921deffa6262c5aULL, 0xd4f35aabcfedfa1fULL}),
+ Frac256({0xb5c13ada0e77829aULL, 0x1d733af522058b16ULL,
+ 0x39a68bb9902d3fdeULL, 0xd744fccad69d6af4ULL}),
+ Frac256({0xb70cfbb1bdf6eb5dULL, 0xc0edda4d891be43dULL,
+ 0xfe873deca3e12babULL, 0xd99d15c278afd7b5ULL}),
+ Frac256({0x613b0d1dbfa0d717ULL, 0x481e1ab725b12d56ULL,
+ 0x3d840d5a9e29aa64ULL, 0xdbfbb797daf23755ULL}),
+ Frac256({0xcc2490c8643ef6b4ULL, 0x01438495eacdf256ULL,
+ 0xdd07a2d9e8466859ULL, 0xde60f4825e0e9123ULL}),
+ Frac256({0x1cb99d3f1ff298a2ULL, 0x224b251b33092002ULL,
+ 0x065895048dd333caULL, 0xe0ccdeec2a94e111ULL}),
+ Frac256({0xfa8fcbb2e85b853fULL, 0xf358a8d368fceaeaULL,
+ 0x09bfe90795980eecULL, 0xe33f8972be8a5a51ULL}),
+ Frac256({0xcefcd5b62a14b818ULL, 0xaacd6065b6e9f6acULL,
+ 0x1e5e8f4a4edbb0ecULL, 0xe5b906e77c8348a8ULL}),
+ Frac256({0x3a1c6473409c261dULL, 0xfe312f84fa665204ULL,
+ 0x791790d0ac70c7ddULL, 0xe8396a503c4bdc68ULL}),
+ Frac256({0x17d8d1e8ca31880bULL, 0xc4faace043b7f91cULL,
+ 0xd02d75b3706e54faULL, 0xeac0c6e7dd24392eULL}),
+ Frac256({0xc8e7c95b06416e6dULL, 0x3787630a764ae4c9ULL,
+ 0x600d2db6a64bfb12ULL, 0xed4f301ed9942b84ULL}),
+ Frac256({0x9392870834f21a53ULL, 0xd4a277eaddaa925cULL,
+ 0x46561cf6948db912ULL, 0xefe4b99bdcdaf5cbULL}),
+ Frac256({0x7c43b0ea5d43228dULL, 0x2cf0b49df0bd70e9ULL,
+ 0xe8980a9cc8f47a4bULL, 0xf281773c59ffb139ULL}),
+ Frac256({0xbdd80329364aa2a0ULL, 0x6f510308677709f5ULL,
+ 0x7b9d0c7aed980fc3ULL, 0xf5257d152486cc2cULL}),
+ Frac256({0xef6797b5a11efb7cULL, 0xe914ffb4723793f1ULL,
+ 0xfe90d496d60fb6eaULL, 0xf7d0df730ad13bb8ULL}),
+ Frac256({0x4844b29bf4af18e8ULL, 0x8006fe21a95d14dcULL,
+ 0x7c25bb14315d7fccULL, 0xfa83b2db722a033aULL}),
+ Frac256({0x9d2285b6754edd61ULL, 0x061b7bb285a60791ULL,
+ 0x853f3a5931e0ee03ULL, 0xfd3e0c0cf486c174ULL}),
+};
+
+// Lookup table for -log2(RD[i]) for i = 0..127 as Frac256,
+// where RD[i] = 2^-8 * ceil(2^8 * (1 - 2^-8) / (1 + i * 2^-7)) is the argument
+// reduction constant from common_constants.h.
+// Generated by Sollya with:
+// > prec = 512;
+// > for i from 0 to 127 do {
+// rd = 2^(-8) * ceil(2^8 * (1 - 2^(-8)) / (1 + i * 2^(-7)));
+// v = round(-log2(rd) * 2^255, 256, RN);
+// print(v);
+// };
+LIBC_INLINE_VAR constexpr Frac256 LOG2_RD_FRAC256[128] = {
+ Frac256({0x0000000000000000ULL, 0x0000000000000000ULL,
+ 0x0000000000000000ULL, 0x0000000000000000ULL}),
+ Frac256({0x6e26f44a7ba5ed8aULL, 0xd4a6b5a62ff68790ULL,
+ 0xb5d184a2c615b70aULL, 0x0172c7ba20f73275ULL}),
+ Frac256({0x262f5ca40b7fa26cULL, 0x3ea8c6b85ced9bdaULL,
+ 0xca906c23ef817e0bULL, 0x02e87dd0c3e6aac6ULL}),
+ Frac256({0x95b7f4fcee92de52ULL, 0xe775fc0ab1c9f028ULL,
+ 0x48f836042de0dc32ULL, 0x04612e39315abe0aULL}),
+ Frac256({0xa3a1e195dfe2d23bULL, 0xfc9f0716c9c77bc6ULL,
+ 0x7970e03f821c75d5ULL, 0x05dce53276563557ULL}),
+ Frac256({0xf021a80d769faed4ULL, 0x1dc18a2c9aa3711eULL,
+ 0x155660710eb2a091ULL, 0x075baf47c7faff81ULL}),
+ Frac256({0xa1c8545273e2d13eULL, 0x2c217796c1fcdb43ULL,
+ 0x631514aef39ce630ULL, 0x08dd9953002a4e86ULL}),
+ Frac256({0x75a862230522a2eaULL, 0xc18cf448bef176a4ULL,
+ 0x050799beaaab2940ULL, 0x0a62b07f3457c407ULL}),
+ Frac256({0x9e9fded86733d3caULL, 0xbaad60b1bf6bcd42ULL,
+ 0x9da288fc615a727dULL, 0x0beb024b67dda633ULL}),
+ Frac256({0x754ee58441dbf839ULL, 0xe1725598f1deb4d6ULL,
+ 0xf22dbbaced44516cULL, 0x0cb0657cd5dbe4f6ULL}),
+ Frac256({0xef894f2363f702c4ULL, 0x4cea550b2568cf82ULL,
+ 0x0d939dceecdd9ce0ULL, 0x0e3da945b878e27dULL}),
+ Frac256({0xe96924afd238c8c9ULL, 0x57c3feaa2020e02eULL,
+ 0x99596a8e2e84c8f4ULL, 0x0fce4aee0e88b274ULL}),
+ Frac256({0xb86f9d5aa2c46088ULL, 0xcc91d77315a3145bULL,
+ 0xa487dfb264b2a99fULL, 0x1097e38ce606492bULL}),
+ Frac256({0x6ffb458008d61f66ULL, 0x7d13d9c8e8a9007cULL,
+ 0xd23af3271ce44c70ULL, 0x122dadc2ab3496d2ULL}),
+ Frac256({0xc5b06c09a5744f6fULL, 0x2d931f14daa7bc7cULL,
+ 0x644ac793db28c412ULL, 0x13c6fb650cde50a1ULL}),
+ Frac256({0xb70ec311ae623b5fULL, 0x47a2302a7c41bfa9ULL,
+ 0xa74a794230202b5bULL, 0x1494f863b8df34baULL}),
+ Frac256({0x319dd5425ceb6df9ULL, 0x2cace034400340f2ULL,
+ 0xa7d70047ddacbac0ULL, 0x1633a8bf437ce10aULL}),
+ Frac256({0xbd7cbee513e8b29cULL, 0x623a9d853e546763ULL,
+ 0x19cb9577a5aea7b3ULL, 0x17046031c79f84beULL}),
+ Frac256({0xa214ecf13179e463ULL, 0x639cad5bff474c1aULL,
+ 0x6a9b7e2df60d2bdcULL, 0x18a8980abfbd3266ULL}),
+ Frac256({0x549b2f738887ca0dULL, 0x13e96e51cbcddd86ULL,
+ 0xaa32d50b40cf8ce7ULL, 0x197c1cb13c7ec085ULL}),
+ Frac256({0x0fc30982c995dbc4ULL, 0x0c5d7f995f1384d0ULL,
+ 0xde69308bc912aa0fULL, 0x1b2602497d53458cULL}),
+ Frac256({0x1766d8107cc2463bULL, 0xadedce27f01160b5ULL,
+ 0xf02bdee0b9f5de06ULL, 0x1bfc67a7fff4cc06ULL}),
+ Frac256({0x911ae4dceeb91543ULL, 0xa0018a21469a1f11ULL,
+ 0x4574e09b954ede83ULL, 0x1dac22d3e441d2feULL}),
+ Frac256({0x45b6ce87115cb69eULL, 0x603bb472c90ec7e1ULL,
+ 0xe40c5e6d7829a1b2ULL, 0x1e857d3d361367bdULL}),
+ Frac256({0xc2ae18e36f1a943bULL, 0xc2d3510c6f5fd964ULL,
+ 0x44ca200650512c0aULL, 0x203b3779f4c3a8bbULL}),
+ Frac256({0x5063aa93fcadc0b9ULL, 0x4805be54f509252bULL,
+ 0x552203798e04bf52ULL, 0x21179c1b2bf46fd8ULL}),
+ Frac256({0x6f2dba8b95f5c7c9ULL, 0x2079020c8a004162ULL,
+ 0x6a1f00babcdb8b0aULL, 0x22d380a6c7e2b0e4ULL}),
+ Frac256({0xe7d3c5ed0de90057ULL, 0x49c1a6f5afa47381ULL,
+ 0xa60108dfb23650c7ULL, 0x23b30593aa4e106bULL}),
+ Frac256({0x40b4cbc998adb82dULL, 0x1e4a0e0dbeb3195dULL,
+ 0x083e072a57679e5aULL, 0x24939a56279ad89aULL}),
+ Frac256({0x34552da7c8bea27fULL, 0xf185be16b7fcc49cULL,
+ 0x4492f1bc6b3e5770ULL, 0x2657fdc6e1dcd0cbULL}),
+ Frac256({0x6d3834bd854918b2ULL, 0xbb77e46da7b566d9ULL,
+ 0x5697a3886ffcccdaULL, 0x273bd1c2ab3edefeULL}),
+ Frac256({0xb402daaee365795fULL, 0x06e4fbd9714b80fbULL,
+ 0xd394fe8cca7d9625ULL, 0x2820c02f87d9451cULL}),
+ Frac256({0x9d7b13b46dda1849ULL, 0x885dce30d5e9e813ULL,
+ 0xc9d0b5ca5a8e7bb5ULL, 0x29edf7659d8b30f1ULL}),
+ Frac256({0x121032713a5003c9ULL, 0xe0b08724e7f04c97ULL,
+ 0x3cd6715512f1784cULL, 0x2ad645cd6af1c939ULL}),
+ Frac256({0x9c9cd255e7b18d6eULL, 0x58b9b632cfed6a65ULL,
+ 0x5adb21d3765b875dULL, 0x2bbfb9e3dd5c1c88ULL}),
+ Frac256({0x362e0218efac2fa7ULL, 0x10fb850fefe496c5ULL,
+ 0x4840199e302970e4ULL, 0x2caa569330c4eed6ULL}),
+ Frac256({0xb9254dfdf461e98cULL, 0x9a4ea80721314ac5ULL,
+ 0xdb502402c94092ccULL, 0x2d961ed0cb91d406ULL}),
+ Frac256({0xbeeea0c45c970f7eULL, 0x3bd29380bd1f0b11ULL,
+ 0xe278bad54ec9e6cfULL, 0x2f713e059e555a63ULL}),
+ Frac256({0x6f816641264b17aeULL, 0xf3f42009108dab4eULL,
+ 0x7561496ab4e4b293ULL, 0x30609b21823fa654ULL}),
+ Frac256({0x4429d9e262f3c8c6ULL, 0xc7395ab7fe8e9835ULL,
+ 0xd536fc5bec1a57b8ULL, 0x315130157f7a64ccULL}),
+ Frac256({0xd2f50c802ff2764aULL, 0xac331ec8776891deULL,
+ 0xe2357ab8cc98c9eeULL, 0x3243001249ba76feULL}),
+ Frac256({0x13bf99bcdd15e820ULL, 0x0c891fcfa03fdd31ULL,
+ 0x51b7e816f77f7b63ULL, 0x33360e552d8d64deULL}),
+ Frac256({0xc20d37d96d02df04ULL, 0x6d860cae2ed327dfULL,
+ 0x2ffa76fafcba2917ULL, 0x351ff2e30214bc30ULL}),
+ Frac256({0x85fb7f8ad3538e0dULL, 0x03d2c5dcbd9ccc4aULL,
+ 0x1ec47c7145831457ULL, 0x3616cfe9e8d01feaULL}),
+ Frac256({0x829c6e703f0fc512ULL, 0x0031e528bbef9eadULL,
+ 0xc4ce7959dfb11374ULL, 0x370ef8af6360dfdfULL}),
+ Frac256({0xe83c947d78828170ULL, 0xac2ed3668bc0c3b9ULL,
+ 0xfa8ae31eec3ba722ULL, 0x380870b3c5fb66f6ULL}),
+ Frac256({0x0cd67db0d807d403ULL, 0x0a9caaa68f883accULL,
+ 0x2bb1588cc9e47f8eULL, 0x39033b85a8bfc871ULL}),
+ Frac256({0x16420cec83688abaULL, 0x7b9e3ba7068f7b86ULL,
+ 0x1c5a0baeb329e73eULL, 0x39ff5cc235a256c5ULL}),
+ Frac256({0x65af196f4ce5b140ULL, 0x6a2513c4f89aa3a6ULL,
+ 0xa96b573a7ed69eedULL, 0x3afcd815786af187ULL}),
+ Frac256({0x319c80ace7918092ULL, 0x2ca495197c5993eaULL,
+ 0x8c8946414c6a14b1ULL, 0x3bfbb13ab0dc5613ULL}),
+ Frac256({0xe91dc50aab435aa9ULL, 0xe06bcb114a0f1719ULL,
+ 0xf6a2b59276887aabULL, 0x3cfbebfca715669dULL}),
+ Frac256({0x63d58c2be0e5b043ULL, 0x99a78444f0d00323ULL,
+ 0x930f8ba9f0570f47ULL, 0x3dfd8c36023f0ab6ULL}),
+ Frac256({0xd66a1a98fa3886e2ULL, 0x55226b72b74410b6ULL,
+ 0xaf2e6fea614b834dULL, 0x3f0095d1a19a0331ULL}),
+ Frac256({0x88c9482df3a84033ULL, 0x25aa4f7373799779ULL,
+ 0x670197a0e8f3ba74ULL, 0x40050ccaf800ca8cULL}),
+ Frac256({0x60ad16b1d436b2cfULL, 0x28336d5fee3ef522ULL,
+ 0xcd9cfff75e149b95ULL, 0x410af52e69f26263ULL}),
+ Frac256({0x14359d48fe9047deULL, 0x1991256093b9dc1cULL,
+ 0xc3fcaf8df7db7c03ULL, 0x42125319ae3bbf05ULL}),
+ Frac256({0x16fadce756865ed5ULL, 0xbae450ac5754344aULL,
+ 0x5cc3da171dd99950ULL, 0x431b2abc31565be7ULL}),
+ Frac256({0x8492c2612d5a8759ULL, 0x34ee848b0b781887ULL,
+ 0x89cd3dd41df9689bULL, 0x442580577b936762ULL}),
+ Frac256({0x3ad4311182915175ULL, 0x60c67a245f78bb52ULL,
+ 0x8283ccdf555594a0ULL, 0x4531583f9a2be203ULL}),
+ Frac256({0x5b3a3aef25dbcdefULL, 0xfdeb5563207896e0ULL,
+ 0x45eba230bf4dbea8ULL, 0x463eb6db8b4f066dULL}),
+ Frac256({0x32087d5975d8fe6bULL, 0xea99e677177c285cULL,
+ 0x02356a22199e7587ULL, 0x474da0a5ad495303ULL}),
+ Frac256({0x04a7b4a6562b4343ULL, 0xddd02924cf83def7ULL,
+ 0x77a639bfdd27aeb2ULL, 0x485e1a2c30df9ea9ULL}),
+ Frac256({0xe61b71e474e99f75ULL, 0x75abfd08dec73401ULL,
+ 0xd7220e04ebb0e2a4ULL, 0x497028118efabeb7ULL}),
+ Frac256({0x2a2437241ee48e55ULL, 0x8e16aa8188fd490aULL,
+ 0xb71b554e74851c3cULL, 0x4a83cf0d01c16e3cULL}),
+ Frac256({0xb80082d9fb082e0aULL, 0x67ff158637af6596ULL,
+ 0x077e50d0c2749c04ULL, 0x4b9913eb013f5ec5ULL}),
+ Frac256({0xb80082d9fb082e0aULL, 0x67ff158637af6596ULL,
+ 0x077e50d0c2749c04ULL, 0x4b9913eb013f5ec5ULL}),
+ Frac256({0x68aa5b4f917d44ffULL, 0xe30b7c2d6ff98938ULL,
+ 0x8925e378d67caee1ULL, 0x4caffb8dc3b9a196ULL}),
+ Frac256({0x642224ca9e7cc367ULL, 0xd122ba0a2e1a73faULL,
+ 0x9a95f528f2c754f3ULL, 0x4dc88aedc1d1ee96ULL}),
+ Frac256({0x313b1a081bc4ea73ULL, 0x09c491c06681cc48ULL,
+ 0x9336b66e4ac8a9deULL, 0x4ee2c71a3e9bb4b6ULL}),
+ Frac256({0x63d1c43831c10bdaULL, 0xf4272036db2b8f69ULL,
+ 0xa293ec16410a6ee4ULL, 0x4ffeb539d3c7579aULL}),
+ Frac256({0xd9740fe9e9c5253fULL, 0x1b73dad61ee48894ULL,
+ 0x202655dbb6b0071eULL, 0x511c5a8b02098837ULL}),
+ Frac256({0xd9740fe9e9c5253fULL, 0x1b73dad61ee48894ULL,
+ 0x202655dbb6b0071eULL, 0x511c5a8b02098837ULL}),
+ Frac256({0x2e147034099ccfcfULL, 0xb03b7d7e6bd0493fULL,
+ 0xe55be97611f87779ULL, 0x523bbc64c5e64350ULL}),
+ Frac256({0x58e4c41655df5f19ULL, 0x740e9519cc12e4a0ULL,
+ 0xbb0e246ec2cef169ULL, 0x535ce0373108b235ULL}),
+ Frac256({0x0834df65b40b7c12ULL, 0xb61832b8bb4c9f7fULL,
+ 0xbfe9dbebf2e8a45dULL, 0x547fcb8c0852f0c0ULL}),
+ Frac256({0x8b32b1fbdcbd3f0cULL, 0x4408786d49566334ULL,
+ 0x613e33c06c95a688ULL, 0x55a4840766d29904ULL}),
+ Frac256({0xdc92a6fefa30f4c6ULL, 0x4d2754039098a562ULL,
+ 0x6da8120164a04966ULL, 0x56cb0f6865c8ea03ULL}),
+ Frac256({0xdc92a6fefa30f4c6ULL, 0x4d2754039098a562ULL,
+ 0x6da8120164a04966ULL, 0x56cb0f6865c8ea03ULL}),
+ Frac256({0x313af26502f8a6cdULL, 0x8dff0ebab8d36942ULL,
+ 0x9a1977b5b995b421ULL, 0x57f37389c9f76d14ULL}),
+ Frac256({0x2d7c1749f6f423c0ULL, 0xdf64485cf695a7a3ULL,
+ 0xdd9926d3f02373c8ULL, 0x591db662b664264cULL}),
+ Frac256({0xaa7293ec76e59d5bULL, 0xc74c7f4741bb1f6fULL,
+ 0xd90b84e721864711ULL, 0x5a49de0764caa121ULL}),
+ Frac256({0xaa7293ec76e59d5bULL, 0xc74c7f4741bb1f6fULL,
+ 0xd90b84e721864711ULL, 0x5a49de0764caa121ULL}),
+ Frac256({0xf662658135c18184ULL, 0x5f0bcac4e6cfec7bULL,
+ 0x748d68b767f88088ULL, 0x5b77f0a9e3f18cfbULL}),
+ Frac256({0x80348894da13c701ULL, 0x49c65b1a5a602129ULL,
+ 0xe718f240e6bcbf3cULL, 0x5ca7f49adc1f1f5aULL}),
+ Frac256({0xdce91e02b40b35d8ULL, 0xc46327805ec0076cULL,
+ 0xed1f4b0d4b62c07cULL, 0x5dd9f04a59e91469ULL}),
+ Frac256({0xdce91e02b40b35d8ULL, 0xc46327805ec0076cULL,
+ 0xed1f4b0d4b62c07cULL, 0x5dd9f04a59e91469ULL}),
+ Frac256({0x5f884b8ddadcf74dULL, 0xf5e3dadf04bd0ff3ULL,
+ 0xf9cb2cc55748a4ccULL, 0x5f0dea489f9fed21ULL}),
+ Frac256({0xa609c77144ded04cULL, 0xb30a91252b56f94aULL,
+ 0x572667587b10ca0dULL, 0x6043e946fd97f5dcULL}),
+ Frac256({0xa609c77144ded04cULL, 0xb30a91252b56f94aULL,
+ 0x572667587b10ca0dULL, 0x6043e946fd97f5dcULL}),
+ Frac256({0xf93ce8eff413527cULL, 0x8554c2c17f944fa1ULL,
+ 0x360c2ae2103c7c0dULL, 0x617bf418b195b338ULL}),
+ Frac256({0x7b3285d76164c890ULL, 0x07d4b4b5500472a5ULL,
+ 0x0b4a9afdc5fabbe4ULL, 0x62b611b3cda69037ULL}),
+ Frac256({0x7b3285d76164c890ULL, 0x07d4b4b5500472a5ULL,
+ 0x0b4a9afdc5fabbe4ULL, 0x62b611b3cda69037ULL}),
+ Frac256({0x599f2a65c6ef7dd4ULL, 0x3edbea325add2accULL,
+ 0x1d9267663010bca1ULL, 0x63f2493226b211bfULL}),
+ Frac256({0x5509d9aded608bccULL, 0xdf7d4115ebc0ee87ULL,
+ 0x1ee1343fe7c9cb4aULL, 0x6530a1d24b136c10ULL}),
+ Frac256({0x5509d9aded608bccULL, 0xdf7d4115ebc0ee87ULL,
+ 0x1ee1343fe7c9cb4aULL, 0x6530a1d24b136c10ULL}),
+ Frac256({0x063711bbcff6a7caULL, 0x34bf67662d61c015ULL,
+ 0x05317356e8d480d0ULL, 0x667122f8818f20fdULL}),
+ Frac256({0xdb09e5e814c87f43ULL, 0x829344b50240232cULL,
+ 0x2ddb71189c56a8f0ULL, 0x67b3d42fd0fc4d02ULL}),
+ Frac256({0xdb09e5e814c87f43ULL, 0x829344b50240232cULL,
+ 0x2ddb71189c56a8f0ULL, 0x67b3d42fd0fc4d02ULL}),
+ Frac256({0x4c8dec7f0307142cULL, 0xca3b2dcc7941a5dfULL,
+ 0x3fe30528818495d6ULL, 0x68f8bd2b10fd80d6ULL}),
+ Frac256({0x841a6fb2da05be09ULL, 0xdb0c195c5da64fbfULL,
+ 0x5ff4edf5f974522eULL, 0x6a3fe5c604297860ULL}),
+ Frac256({0x841a6fb2da05be09ULL, 0xdb0c195c5da64fbfULL,
+ 0x5ff4edf5f974522eULL, 0x6a3fe5c604297860ULL}),
+ Frac256({0x6f295eb94b4ff3f4ULL, 0x524c383b17757feeULL,
+ 0xc716be9bc093ec11ULL, 0x6b8956067c08b2ceULL}),
+ Frac256({0xd50c7b85bba151bcULL, 0x6ac0ed3b39116e44ULL,
+ 0xae0d3f8a58b459b2ULL, 0x6cd5161d8751e5e0ULL}),
+ Frac256({0xd50c7b85bba151bcULL, 0x6ac0ed3b39116e44ULL,
+ 0xae0d3f8a58b459b2ULL, 0x6cd5161d8751e5e0ULL}),
+ Frac256({0xcee3b58a450ab308ULL, 0x7822b754be5b62abULL,
+ 0x5babcf87c69ea8a5ULL, 0x6e232e68aad484a1ULL}),
+ Frac256({0xcee3b58a450ab308ULL, 0x7822b754be5b62abULL,
+ 0x5babcf87c69ea8a5ULL, 0x6e232e68aad484a1ULL}),
+ Frac256({0x7ea793f02baad929ULL, 0xb0c4015f8fdff17dULL,
+ 0xd843902f5aad7542ULL, 0x6f73a77325861c69ULL}),
+ Frac256({0xe19972fa315f9621ULL, 0xd5e23a4c7b4973b0ULL,
+ 0xa1251311eb06fd8aULL, 0x70c689f7402d26f1ULL}),
+ Frac256({0xe19972fa315f9621ULL, 0xd5e23a4c7b4973b0ULL,
+ 0xa1251311eb06fd8aULL, 0x70c689f7402d26f1ULL}),
+ Frac256({0xab2afce4184639aeULL, 0x4df1d7bf78e23ef9ULL,
+ 0x269d2c8d7342a3c3ULL, 0x721bdedfa92a22ceULL}),
+ Frac256({0xab2afce4184639aeULL, 0x4df1d7bf78e23ef9ULL,
+ 0x269d2c8d7342a3c3ULL, 0x721bdedfa92a22ceULL}),
+ Frac256({0x261c1d2c4d56a905ULL, 0xb94dd8c8f0834794ULL,
+ 0xc6e6db59262e2e6fULL, 0x7373af48dce652a8ULL}),
+ Frac256({0xdb6ddbdb00d73945ULL, 0x1ad7c182a716c033ULL,
+ 0x99d63ecf5dd4529eULL, 0x74ce04829b7674c1ULL}),
+ Frac256({0xdb6ddbdb00d73945ULL, 0x1ad7c182a716c033ULL,
+ 0x99d63ecf5dd4529eULL, 0x74ce04829b7674c1ULL}),
+ Frac256({0x22ba4e8b413991d3ULL, 0x95b97a0e1d121d02ULL,
+ 0xfd9776f25acec4acULL, 0x762ae8116c071e93ULL}),
+ Frac256({0x22ba4e8b413991d3ULL, 0x95b97a0e1d121d02ULL,
+ 0xfd9776f25acec4acULL, 0x762ae8116c071e93ULL}),
+ Frac256({0x7420ca3f82a950a2ULL, 0xa3757861c1116ee5ULL,
+ 0x1845a2a3336f47ccULL, 0x778a63b02eb032a6ULL}),
+ Frac256({0x7420ca3f82a950a2ULL, 0xa3757861c1116ee5ULL,
+ 0x1845a2a3336f47ccULL, 0x778a63b02eb032a6ULL}),
+ Frac256({0x481eb4627658b4afULL, 0x16a73e812a9e4565ULL,
+ 0xc1c1e586711df5eaULL, 0x78ec8151bd552842ULL}),
+ Frac256({0x481eb4627658b4afULL, 0x16a73e812a9e4565ULL,
+ 0xc1c1e586711df5eaULL, 0x78ec8151bd552842ULL}),
+ Frac256({0xfce168eaef943079ULL, 0xce4c86d28e4be331ULL,
+ 0xb27e43da520fbdb7ULL, 0x7a514b229c409e33ULL}),
+ Frac256({0xfce168eaef943079ULL, 0xce4c86d28e4be331ULL,
+ 0xb27e43da520fbdb7ULL, 0x7a514b229c409e33ULL}),
+ Frac256({0x59576b46a119898aULL, 0x254008aeb4167a33ULL,
+ 0x9faebec15b2e2b43ULL, 0x7bb8cb8abb32fb44ULL}),
+ Frac256({0x59576b46a119898aULL, 0x254008aeb4167a33ULL,
+ 0x9faebec15b2e2b43ULL, 0x7bb8cb8abb32fb44ULL}),
+ Frac256({0x10e7d4a6c2ef9745ULL, 0x3ead121a489569ffULL,
+ 0xb23b03bdcfdea0d7ULL, 0x7d230d2f47a5baceULL}),
+ Frac256({0x10e7d4a6c2ef9745ULL, 0x3ead121a489569ffULL,
+ 0xb23b03bdcfdea0d7ULL, 0x7d230d2f47a5baceULL}),
+ Frac256({0xa828a9bde1ec7e79ULL, 0xe33209b70d9a5be1ULL,
+ 0x071c84ffe86b0bbbULL, 0x7e901af4910f7ae8ULL}),
+ Frac256({0x0000000000000000ULL, 0x0000000000000000ULL,
+ 0x0000000000000000ULL, 0x8000000000000000ULL}),
+};
+
+// Accurate log2(x) in 256-bit precision reusing range reduction from fast pass:
+// x = 2^x_e * m_x, with 1 <= m_x < 2
+// r = RD[idx_x]
+// dx = r * m_x - 1, with -2^-8 <= dx < 2^-7
+// log2(x) = x_e + (-log2(r)) + log2(1 + dx).
+// The value -log2(r) is provided by LOG2_RD_FRAC256[idx_x].
+// log2(1 + dx) is evaluated with a degree-30 minimax polynomial using a 3-tier
+// Horner scheme (Frac64, Frac128, and Frac256).
+LIBC_INLINE DFloat256 log2_f256(int x_e, unsigned idx_x, double dx) {
+ using FPBits = fputil::FPBits<double>;
+
+ if (dx == 0.0) {
+ Frac256 log2_m = LOG2_RD_FRAC256[idx_x];
+ DFloat256 m_df(Sign::POS, -255, log2_m);
+ return fputil::quick_add(DFloat256(static_cast<double>(x_e)), m_df);
+ }
+
+ double abs_dx = (dx < 0.0) ? -dx : dx;
+ FPBits bits(abs_dx);
+ // Append hidden bit.
+ uint64_t mant = bits.get_mantissa() | (1ULL << 52);
+ int shift = 255 + bits.get_exponent() - 52;
+ Frac256 y =
+ (shift >= 0) ? Frac256((UInt<256>(mant) << shift).val) : Frac256(0);
+
+ // Evaluate log2(1 + dx) using Horner scheme in 3 stages:
+ // - Degree 18-30: 64-bit precision evaluation.
+ // - Degree 8-17: 128-bit precision evaluation.
+ // - Degree 0-7: 256-bit precision evaluation.
+ //
+ // Since y = |dx| <= 2^-7, truncation errors at degree k is bounded by:
+ // y^k <= 2^(-7k).
+ //
+ // Step 1: 64-bit evaluation, truncation error at degree 18 is bounded by:
+ // 2^-63 * y^18 <= 2^-189.
+ Frac64 y64 = y.to_frac64();
+ Frac64 p64 = LOG2_POLY_256[30].to_frac64();
+ if (dx >= 0.0) {
+ for (int k = 29; k >= 18; --k)
+ p64 = LOG2_POLY_256[k].to_frac64() - ((p64 * y64) << 1);
+ } else {
+ for (int k = 29; k >= 18; --k)
+ p64 = LOG2_POLY_256[k].to_frac64() + ((p64 * y64) << 1);
+ }
+
+ // Step 2: 128-bit evaluation, truncation error at degree 8 is bounded by:
+ // 2^-127 * y^8 <= 2^-183.
+ Frac128 y128 = y.to_frac128();
+ Frac128 p128({0, p64.val[0]});
+ if (dx >= 0.0) {
+ for (int k = 17; k >= 8; --k)
+ p128 = LOG2_POLY_256[k].to_frac128() - ((p128 * y128) << 1);
+ } else {
+ for (int k = 17; k >= 8; --k)
+ p128 = LOG2_POLY_256[k].to_frac128() + ((p128 * y128) << 1);
+ }
+
+ // Step 3: 256-bit evaluation.
+ Frac256 p({0, 0, p128.val[0], p128.val[1]});
+ if (dx >= 0.0) {
+ for (int k = 7; k >= 0; --k)
+ p = LOG2_POLY_256[k] - ((p * y) << 1);
+ } else {
+ for (int k = 7; k >= 0; --k)
+ p = LOG2_POLY_256[k] + ((p * y) << 1);
+ }
+
+ // log2(1 + dx) = 2 * y * P(y):
+ Frac256 log2_1p = (y * p) << 2;
+ Frac256 log2_rd = LOG2_RD_FRAC256[idx_x];
+ Frac256 log2_m = (dx >= 0.0) ? (log2_rd + log2_1p) : (log2_rd - log2_1p);
+
+ DFloat256 m_df(Sign::POS, -255, log2_m);
+ return fputil::quick_add(DFloat256(static_cast<double>(x_e)), m_df);
+}
+
+// Range reduction:
+// k = round(z * 64)
+// hi = k >> 6
+// idx = k & 0x3f
+// lo = z - k * 2^-6, with |lo| <= 2^-7.
+// Then:
+// 2^z = 2^hi * EXP2_MID_FRAC256[idx] * 2^lo
+// = 2^hi * EXP2_MID_FRAC256[idx] * (1 + lo * P(lo)).
+LIBC_INLINE DFloat256 exp2_f256(const DFloat256 &z) {
+ double z_d = static_cast<double>(z);
+ double z_scaled = z_d * 64.0;
+ double kd = fputil::nearest_integer(z_scaled);
+ int k = static_cast<int>(kd);
+
+ int hi = k >> 6;
+ unsigned idx = static_cast<unsigned>(k & 0x3f);
+
+ DFloat256 kd_f256(kd * 0x1.0p-6);
+ DFloat256 lo = fputil::quick_add(z, -kd_f256);
+
+ Frac256 m = EXP2_MID_FRAC256[idx];
+ if (LIBC_UNLIKELY(lo.mantissa.is_zero()))
+ return DFloat256(Sign::POS, hi - 255, m);
+
+ bool lo_is_neg = (lo.sign == Sign::NEG);
+ int shift = -255 - lo.exponent;
+ Frac256 u = (shift < 256) ? Frac256((lo.mantissa >> shift).val) : Frac256(0);
+
+ // Evaluate 2^lo - 1 = lo * P(lo) using Horner scheme in 3 stages:
+ // - Degree 12-21: 64-bit precision evaluation.
+ // - Degree 6-11: 128-bit precision evaluation.
+ // - Degree 0-5: 256-bit precision evaluation.
+ //
+ // Since u = |lo| <= 2^-7, truncation errors at degree i is bounded by:
+ // u^i <= 2^(-7i).
+ //
+ // Step 1: 64-bit evaluation, truncation error at degree 12 is bounded by:
+ // 2^-63 * u^12 <= 2^-147.
+ Frac64 u64 = u.to_frac64();
+ Frac64 p64 = EXP2_POLY_256[21].to_frac64();
+ if (lo_is_neg) {
+ for (int i = 20; i >= 12; --i)
+ p64 = EXP2_POLY_256[i].to_frac64() - ((p64 * u64) << 1);
+ } else {
+ for (int i = 20; i >= 12; --i)
+ p64 = EXP2_POLY_256[i].to_frac64() + ((p64 * u64) << 1);
+ }
+
+ // Step 2: 128-bit evaluation, truncation error at degree 6 is bounded by:
+ // 2^-127 * u^6 <= 2^-169.
+ Frac128 u128 = u.to_frac128();
+ Frac128 p128({0, p64.val[0]});
+ if (lo_is_neg) {
+ for (int i = 11; i >= 6; --i)
+ p128 = EXP2_POLY_256[i].to_frac128() - ((p128 * u128) << 1);
+ } else {
+ for (int i = 11; i >= 6; --i)
+ p128 = EXP2_POLY_256[i].to_frac128() + ((p128 * u128) << 1);
+ }
+
+ // Step 3: 256-bit evaluation.
+ Frac256 p({0, 0, p128.val[0], p128.val[1]});
+ if (lo_is_neg) {
+ for (int i = 5; i >= 0; --i)
+ p = EXP2_POLY_256[i] - ((p * u) << 1);
+ } else {
+ for (int i = 5; i >= 0; --i)
+ p = EXP2_POLY_256[i] + ((p * u) << 1);
+ }
+
+ // Reconstruction:
+ // t = u * P(u) ~ 2^u - 1 as Frac256
+ // mt = m * t ~ m * (2^u - 1) as Frac256
+ // m_final = m +- mt ~ m * (1 +- t) ~ m * 2^lo
+ Frac256 t = (u * p) << 1;
+ Frac256 mt = (m * t) << 1;
+ Frac256 m_final = lo_is_neg ? (m - mt) : (m + mt);
+
+ if ((m_final.val[3] & (1ULL << 63)) == 0) {
+ m_final = m_final << 1;
+ --hi;
+ }
+
+ return DFloat256(Sign::POS, hi - 255, m_final);
+}
+
+// Accurate pow(x, y) reusing range reduction parameters in 256-bit precision.
+// Lauter & Lefevre (2009) showed that if x^y is not an exact 54-bit number, the
+// distance to the nearest 54-bit boundary is at least:
+// |x^y - o_54(x^y)| / x^y >= 2^-114.
+// Since exact 54-bit boundaries are already handled in pow_accurate, the
+// 256-bit evaluation with error < 2^-250 is sufficient to correctly round all
+// remaining cases.
+LIBC_INLINE double pow_accurate_256(double y, bool is_neg, int x_e,
+ unsigned idx_x, double dx) {
+ DFloat256 log2_x = log2_f256(x_e, idx_x, dx);
+ DFloat256 y_f256(y);
+ DFloat256 z = fputil::quick_mul(y_f256, log2_x);
+
+ // For 0 < |z| <= 2^-55, x^y is between 1 - 2^-54 and 1 + 2^-53.
+ if (LIBC_UNLIKELY(!z.mantissa.is_zero() && z.exponent + 255 <= -55)) {
+ volatile double one = 1.0;
+ volatile double eps = (z.sign == Sign::NEG) ? -0x1.0p-100 : 0x1.0p-100;
+ double res = one + eps;
+ return is_neg ? -res : res;
+ }
+
+ DFloat256 r = exp2_f256(z);
+ if (is_neg)
+ r.sign = Sign::NEG;
+
+ int unbiased_exp = r.exponent + 255;
+ if (LIBC_UNLIKELY(unbiased_exp >= 1024))
+ return set_overflow(is_neg);
+
+ double res = static_cast<double>(r);
+ if (LIBC_UNLIKELY(fputil::FPBits<double>(res).is_inf()))
+ return set_overflow(is_neg);
+
+ return res;
+}
+
+} // namespace pow_internal
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_256_H
diff --git a/libc/src/__support/math/pow_fast.h b/libc/src/__support/math/pow_fast.h
new file mode 100644
index 0000000000000..5960bea60642b
--- /dev/null
+++ b/libc/src/__support/math/pow_fast.h
@@ -0,0 +1,304 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Fast implementation header for double-precision pow(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_FAST_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_FAST_H
+
+#include "src/__support/FPUtil/FPBits.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/common.h"
+#include "src/__support/macros/config.h"
+#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"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+LIBC_INLINE double pow(double x, double y) {
+ using namespace pow_internal;
+ using FPBits = fputil::FPBits<double>;
+
+ FPBits xbits(x), ybits(y);
+ uint64_t x_u = xbits.uintval();
+ uint64_t y_u = ybits.uintval();
+ uint64_t y_a = ybits.abs().uintval();
+
+ if (LIBC_UNLIKELY((x_u & 0x0003'FFFF'FFFF'FFFF) == 0 ||
+ (y_u & 0x000F'FFFF'FFFF'FFFF) == 0)) {
+ if (auto r = check_special_inputs(x, y); LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ }
+
+ double e_x = static_cast<double>(xbits.get_exponent());
+ uint64_t x_mant = xbits.get_mantissa();
+ bool is_neg = false;
+ double sign_d = 1.0;
+
+ if (LIBC_UNLIKELY(y_a <= Y_LOWER_BOUND || y_a >= Y_UPPER_BOUND ||
+ x_u >= FPBits::inf().uintval() ||
+ x_u < FPBits::min_normal().uintval())) {
+ if (auto r = check_exceptional_cases(x, y, e_x, x_mant, is_neg, sign_d);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ }
+
+ // x^y = 2^( y * log2(x) )
+ // = 2^( y * ( e_x + log2(m_x) ) )
+ // First we compute log2(x) = e_x + log2(m_x)
+
+ // Extract exponent field of x.
+
+ // Use the highest 7 fractional bits of m_x as the index for look up tables.
+ unsigned idx_x = static_cast<unsigned>(x_mant >> (FPBits::FRACTION_LEN - 7));
+ // Add the hidden bit to the mantissa.
+ // 1 <= m_x < 2
+ FPBits m_x = FPBits(x_mant | 0x3ff0'0000'0000'0000);
+
+ // Reduced argument for log2(m_x):
+ // dx = r * m_x - 1.
+ // The computation is exact, and -2^-8 <= dx < 2^-7.
+ // Then m_x = (1 + dx) / r, and
+ // log2(m_x) = log2( (1 + dx) / r )
+ // = log2(1 + dx) - log2(r).
+
+ // See the overview comment in pow.h for the detailed error analysis.
+ // To ensure the relative error is strictly < 1 ULP across the entire double
+ // precision range, we evaluate the prefix log2(e) * (dx - dx^2/2) accurately
+ // using DoubleDouble (n = 2), and approximate higher terms with dx^3 * P(dx).
+
+ // Constants for C_k = (-1)^(k-1) / (k * log(2)) in DoubleDouble:
+ // C1 = 1 / log(2) = log2(e)
+ constexpr DoubleDouble C1 = {0x1.777d0ffda0d24p-56, 0x1.71547652b82fep0};
+ // C2 = -1 / (2 * log(2)) = -log2(e) / 2
+ constexpr DoubleDouble C2 = {-0x1.777d0ffda0d24p-57, -0x1.71547652b82fep-1};
+
+ // Degree-5 polynomial approximation for:
+ // P(dx) ~ (log2(1 + dx) - (dx - dx^2/2) / log(2)) / dx^3
+ // Generated by Sollya with:
+ // > P = fpminimax((log2(1 + x) - (x - x^2/2)/log(2))/x^3, 5, [|D...|],
+ // [-2^-8, 2^-7]);
+ // > dirtyinfnorm((log2(1 + x) - (x - x^2/2)/log(2))/x - x^2*P,
+ // [-2^-8, 2^-7]);
+ // 0x1.cfc09...p-67
+ constexpr double COEFFS[] = {
+ 0x1.ec709dc3a03fcp-2, -0x1.71547652b7546p-2, 0x1.2776c50fc020bp-2,
+ -0x1.ec70a327cae39p-3, 0x1.a616f87d961abp-3, -0x1.6d89893e348ddp-3,
+ };
+
+ // Perform exact range reduction.
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double dx = fputil::multiply_add(RD[idx_x], m_x.get_val(), -1.0); // Exact
+#else
+ double c = FPBits(m_x.uintval() & 0x3fff'e000'0000'0000).get_val();
+ double dx =
+ fputil::multiply_add(RD[idx_x], m_x.get_val() - c, CD[idx_x]); // Exact
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ // Error-free transformation for r = C1 + dx * C2:
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double r_hi = fputil::multiply_add(dx, C2.hi, C1.hi);
+ double r_lo = fputil::multiply_add(dx, C2.hi, C1.hi - r_hi);
+#else
+ DoubleDouble dx_c2 = fputil::exact_mult(dx, C2.hi);
+ DoubleDouble r_sum = fputil::exact_add(C1.hi, dx_c2.hi);
+ double r_hi = r_sum.hi;
+ double r_lo = r_sum.lo + dx_c2.lo;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ r_lo = fputil::multiply_add(dx, C2.lo, r_lo + C1.lo);
+
+ // dx * (r_hi + r_lo) = dx * C1 + dx^2 * C2
+ DoubleDouble dx_r = fputil::exact_mult(dx, r_hi);
+ dx_r.lo = fputil::multiply_add(dx, r_lo, dx_r.lo);
+
+ // Evaluate polynomial tail P(dx) using Estrin's scheme:
+ double dx2 = dx * dx;
+ double c0 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
+ double c1 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
+ double c2 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);
+
+ double dx4 = dx2 * dx2;
+ double d0 = fputil::multiply_add(dx2, c1, c0);
+ double p = fputil::multiply_add(dx4, c2, d0);
+
+ double dx3 = dx2 * dx;
+ double log2_x_lo =
+ fputil::multiply_add(dx3, p, dx_r.lo + LOG2_R_DD[idx_x].lo);
+
+ DoubleDouble log2_x_hi =
+ fputil::exact_add(e_x + LOG2_R_DD[idx_x].hi, dx_r.hi);
+
+ DoubleDouble log2_x = fputil::exact_add(log2_x_hi.hi, log2_x_lo);
+ log2_x.lo += log2_x_hi.lo;
+
+ // To compute 2^(y * log2(x)), we break the exponent into 3 parts:
+ // y * log2(x) = hi + mid + lo, where
+ // hi is an integer
+ // mid * 2^6 is an integer
+ // |lo| <= 2^-7
+ // Then:
+ // x^y = 2^(y * log2(x)) = 2^hi * 2^mid * 2^lo,
+ // In which 2^mid is obtained from a look-up table of size 2^6 = 64 elements,
+ // and 2^lo ~ 1 + lo * P(lo).
+ // Thus, we have:
+ // hi + mid = 2^-6 * round( 2^6 * y * log2(x) )
+ // If we restrict the output such that |hi| < 512, (hi + mid) uses (9 + 6)
+ // bits, hence, if we use double precision to perform
+ // round( 2^6 * y * log2(x))
+ // the lo part is bounded by 2^-7 + 2^(-(52 - 15)) = 2^-7 + 2^-37
+
+ // In the following computations:
+ // y6 = 2^6 * y
+ // hm = 2^6 * (hi + mid) = round(2^6 * y * log2(x)) ~ round(y6 * s)
+ // lo6 = 2^6 * lo = 2^6 * (y - (hi + mid)) = y6 * log2(x) - hm.
+ constexpr double SCALE = 0x1.0p6;
+ double y6 = y * SCALE; // Exact.
+
+ DoubleDouble y6_log2_x = fputil::exact_mult(y6, log2_x.hi);
+ y6_log2_x.lo = fputil::multiply_add(y6, log2_x.lo, y6_log2_x.lo);
+
+ // Check overflow/underflow.
+ double scale = 1.0;
+
+ // |2^(hi + mid) - exp2_hi_mid| <= ulp(exp2_hi_mid) / 2
+
+ // The fast computation for 2^hi below requires that:
+ // |hi| < 512, or equivalently, |hm| < 512 * 2^6.
+ // This guarantees that the biased exponent:
+ // exp_biased = (hm_i >> 6) + EXP_BIAS = hi + 1023
+ // is strictly within the normal double range:
+ // 1023 - 511 <= exp_biased <= 1023 + 511, or 512 <= exp_biased <= 1534.
+ // Hence, 2^hi is always a normal, non-zero, finite power of 2, and the
+ // multiplication upper * exp2_hi is exact and never underflows or overflows.
+ //
+ // From the edge case checks above:
+ // 2^(-54) / 1074 <= |y| <= 1075 * 2^53, and
+ // |log_2(1 - 2^(-53))| <= |log2(x)| <= 1074.
+ // So their product is bounded by:
+ // 2^-117 < |y * log2(x)| < 2^74.
+ //
+ // The meaningful range for y * log2(x) in double precision is:
+ // -1075 <= y * log2(x) <= 1024.
+ // Any value > 1024 overflows, and any value < -1075 underflows.
+ //
+ // When |y * log2(x)| >= 511, we shift the exponent by an offset S:
+ // y * log2(x) = (y * log2(x) - S) + S
+ // and multiply by scale = 2^S at the end.
+ // To ensure the shifted exponent (y * log2(x) - S) stays within [-511, 511]:
+ // - For positive range [511, 1024]:
+ // 511 - S > -511 ==> S < 1022
+ // 1024 - S < 511 ==> S > 513
+ // - For negative range [-1076, -511]:
+ // -511 + S < 511 ==> S <= 1022
+ // -1076 + S > -511 ==> S >= 565
+ // Combined, the shift must satisfy: 565 <= S <= 1022.
+ // We choose S = 600, which yields:
+ // y * log2(x) - 600 in [-89, 424] for the positive range, and
+ // y * log2(x) + 600 in [-476, 89] for the negative range,
+ // both fitting well within [-511, 511].
+ //
+ // For exponents that are completely out of range:
+ // y * log2(x) > 1024 or y * log2(x) < -1076,
+ // the product can be as large as 2^74, so subtracting or adding 600 * 64
+ // would still overflow a 32-bit int when computing hm_i.
+ // We return early with overflow or underflow in these cases.
+ //
+ // Alternatively, for less branching or in SIMD / vector implementations, one
+ // could clamp y6_log2_x.hi to:
+ // - UPPER_EXP_BOUND (511 * 64) for overflow, which when multiplied by
+ // scale = 2^600 yields 2^1111 and correctly overflows.
+ // - -500 * 64 for underflow, which when multiplied by scale = 2^-600 yields
+ // 2^-1100 and correctly underflows.
+ // Note that in a correctly rounded version, such clamping can cause the
+ // Ziv accuracy test to fail on the clamped exponent, unnecessarily
+ // triggering the slow accurate paths for overflow or underflow cases.
+
+ constexpr double UPPER_EXP_BOUND = 511.0 * SCALE;
+ if (LIBC_UNLIKELY(FPBits(y6_log2_x.hi).abs().get_val() >= UPPER_EXP_BOUND)) {
+ if (FPBits(y6_log2_x.hi).sign() == Sign::POS) {
+ if (y6_log2_x.hi > 1024.0 * SCALE)
+ return set_overflow(is_neg);
+ scale = 0x1.0p600;
+ y6_log2_x.hi -= 600.0 * SCALE;
+ } else {
+ if (y6_log2_x.hi < -1076.0 * SCALE)
+ return set_underflow(is_neg);
+ scale = 0x1.0p-600;
+ y6_log2_x.hi += 600.0 * SCALE;
+ }
+ }
+
+ double hm = fputil::nearest_integer(y6_log2_x.hi);
+
+ // lo6 = 2^6 * lo.
+ DoubleDouble lo6 = fputil::exact_add(y6_log2_x.hi - hm, y6_log2_x.lo);
+
+ int hm_i = static_cast<int>(hm);
+ unsigned idx_y = static_cast<unsigned>(hm_i) & 0x3f;
+
+ // 2^hi
+ int64_t exp2_hi_i =
+ static_cast<int64_t>(static_cast<uint64_t>((hm_i >> 6) + FPBits::EXP_BIAS)
+ << FPBits::FRACTION_LEN);
+ double exp2_hi = FPBits(static_cast<uint64_t>(exp2_hi_i)).get_val();
+
+ // 2^mid
+ DoubleDouble exp2_mid{EXP2_MID1[idx_y].mid * sign_d,
+ EXP2_MID1[idx_y].hi * sign_d};
+
+ // Polynomial expansion for 2^(lo6/64):
+ // 2^(lo6/64) ~ 1 + lo6 * (log(2)/64) + lo6^2 * P(lo6)
+ // The linear term is computed in DoubleDouble, and the degree-3 polynomial
+ // P(lo6) is evaluated with standard double precision.
+ //
+ // hi and lo parts of log(2)/64, generated by Sollya with:
+ // > a = D(log(2)/64);
+ // > b = D(log(2)/64 - a);
+ constexpr DoubleDouble LOG_2_OVER_64 = {0x1.abc9e3b39803fp-62,
+ 0x1.62e42fefa39efp-7};
+
+ // Degree-3 polynomial approximation for (2^(lo6/64) - 1 - lo6*log(2)/64) /
+ // lo6^2: Generated by Sollya with:
+ // > f = (2^(x/64) - 1 - x*log(2)/64) / x^2;
+ // > P = fpminimax(f, 3, [|D...|], [-0.5, 0.5]);
+ // > dirtyinfnorm((f - P) * x^2, [-0.5, 0.5]);
+ // 0x1.439a...p-58
+ constexpr double EXP2_COEFFS[] = {
+ 0x1.ebfbdff82bb76p-15, 0x1.c6b08d704acbep-23, 0x1.3b2acb2c56647p-31,
+ 0x1.5d8803cdafb8ep-40};
+
+ DoubleDouble lo_log_2 = fputil::quick_mult(lo6, LOG_2_OVER_64);
+ DoubleDouble lo_log_2_p1 = fputil::exact_add(1.0, lo_log_2.hi);
+ lo_log_2_p1.lo += lo_log_2.lo;
+
+ double lo6_sq = lo6.hi * lo6.hi;
+ double e0 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[1], EXP2_COEFFS[0]);
+ double e1 = fputil::multiply_add(lo6.hi, EXP2_COEFFS[3], EXP2_COEFFS[2]);
+ double pp = fputil::multiply_add(lo6_sq, e1, e0);
+
+ lo_log_2_p1.lo = fputil::multiply_add(lo6_sq, pp, lo_log_2_p1.lo);
+ DoubleDouble r = fputil::quick_mult(exp2_mid, lo_log_2_p1);
+
+ double tmp1 = (r.hi + r.lo) * exp2_hi;
+ return tmp1 * scale;
+}
+
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_FAST_H
diff --git a/libc/src/__support/math/pow_utils.h b/libc/src/__support/math/pow_utils.h
new file mode 100644
index 0000000000000..a5d026737168e
--- /dev/null
+++ b/libc/src/__support/math/pow_utils.h
@@ -0,0 +1,671 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Common utilities and tables for double-precision pow(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_UTILS_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#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/double_double.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/common_constants.h"
+#include "src/__support/math/exp10.h"
+#include "src/__support/math/exp2.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},
+};
+
+// Check if x is an odd integer: the lowest set bit must be at the unit
+// position:
+// x_e + lsb == UNIT_EXPONENT.
+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);
+}
+
+// Check if x is an integer: the lowest set bit must be at or above the unit
+// position:
+// x_e + lsb >= UNIT_EXPONENT.
+LIBC_INLINE bool is_integer(double x) {
+ if (x == 0.0)
+ return true;
+ 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 double set_overflow(bool is_neg) {
+ fputil::set_errno_if_required(ERANGE);
+#if defined(LIBC_TARGET_CPU_HAS_FPU_DOUBLE)
+ volatile double x = is_neg ? -0x1.0p1023 : 0x1.0p1023;
+ return x * 2.0;
+#else
+ using FPBits = fputil::FPBits<double>;
+ fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_TOWARDZERO)
+ return is_neg ? -FPBits::max_normal().get_val()
+ : FPBits::max_normal().get_val();
+ if (rounding == FE_DOWNWARD)
+ return is_neg ? -FPBits::inf().get_val() : FPBits::max_normal().get_val();
+ if (rounding == FE_UPWARD)
+ return is_neg ? -FPBits::max_normal().get_val() : FPBits::inf().get_val();
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return is_neg ? -FPBits::inf().get_val() : FPBits::inf().get_val();
+#endif // LIBC_TARGET_CPU_HAS_FPU_DOUBLE
+}
+
+LIBC_INLINE double set_underflow(bool is_neg) {
+ fputil::set_errno_if_required(ERANGE);
+#if defined(LIBC_TARGET_CPU_HAS_FPU_DOUBLE)
+ volatile double x = is_neg ? -0x1.0p-1022 : 0x1.0p-1022;
+ return x * 0x1.0p-100;
+#else
+ using FPBits = fputil::FPBits<double>;
+ fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_UPWARD && !is_neg)
+ return FPBits::min_subnormal().get_val();
+ if (rounding == FE_DOWNWARD && is_neg)
+ return -FPBits::min_subnormal().get_val();
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return is_neg ? -0.0 : 0.0;
+#endif // LIBC_TARGET_CPU_HAS_FPU_DOUBLE
+}
+
+// Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We
+// assume further that
+// 1 <= mid < 2,
+// mid + lo < 2,
+// and |lo| << mid.
+//
+// Notice that, if 0 < x < 2^-1022:
+// double( 2^-1022 + x ) - 2^-1022 = double(x),
+// and if -2^-1022 < x < 0:
+// double(-2^-1022 + x ) - -2^-1022 = double(x).
+//
+// So if we scale |x| up by 2^1022, we can use the "+1 trick" (or "-1 trick"
+// when x < 0):
+// double( 1.0 + 2^1022 * x) - 1.0, for x > 0, and
+// double(-1.0 + 2^1022 * x) - -1.0, for x < 0,
+// to test how x is rounded in the denormal range. By adding +-1.0, the
+// results will have similar rounding points as denormal outputs, aligning the
+// least significant bit of the mantissa with 2^-52 * 2^-1022 = 2^-1074.
+//
+// Finally, subtracting 1023 in the exponent field converts the rounded normal
+// value into the subnormal value.
+template <bool SKIP_ZIV_TEST = false>
+LIBC_INLINE constexpr cpp::optional<double>
+ziv_test_denorm(int hi, double mid, double lo, double err, bool is_neg) {
+ using FPBits = typename fputil::FPBits<double>;
+
+ uint64_t scale_u = static_cast<uint64_t>(hi + 2045) << FPBits::FRACTION_LEN;
+ double exp_scale = cpp::bit_cast<double>(scale_u);
+
+ double mid_hi = mid * exp_scale;
+ double lo_scaled = lo * exp_scale;
+
+ double extra_factor = 0.0;
+ uint64_t scale_down = 0x3FE0'0000'0000'0000ULL; // 1022 in the exponent field.
+
+ // Result is denormal if |mid_hi + lo_scaled| < 1.0.
+ // In that case, add +-1.0 to mimic denormal rounding points.
+ if (is_neg) {
+ if ((-1.0 - mid_hi) < lo_scaled) {
+ // Extra rounding step is needed, which adds more rounding errors.
+ extra_factor = -1.0;
+ scale_down = 0x3FF0'0000'0000'0000ULL; // 1023 in the exponent field.
+ err += 0x1.0p-52;
+ }
+ } else {
+ if ((1.0 - mid_hi) > lo_scaled) {
+ // Extra rounding step is needed, which adds more rounding errors.
+ extra_factor = 1.0;
+ scale_down = 0x3FF0'0000'0000'0000ULL; // 1023 in the exponent field.
+ err += 0x1.0p-52;
+ }
+ }
+
+ // By adding +-1.0, the hardware adder performs rounding at the exact same
+ // bit position as denormal rounding.
+ if constexpr (SKIP_ZIV_TEST) {
+ double r = extra_factor + (mid_hi + lo_scaled);
+ return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(r) - scale_down);
+ } else {
+ double err_scaled = err * exp_scale;
+
+ double lo_u = lo_scaled + err_scaled;
+ double lo_l = lo_scaled - err_scaled;
+
+ double upper = extra_factor + (mid_hi + lo_u);
+ double lower = extra_factor + (mid_hi + lo_l);
+
+ if (LIBC_LIKELY(upper == lower)) {
+ return cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(upper) - scale_down);
+ }
+
+ return cpp::nullopt;
+ }
+}
+
+// Upper bound for y = 1075 / |log2(1 - 2^(-53)|, generated by Sollya:
+// > y = round(-1075 / log2(1 - 2^(-53)), D, RU);
+// > y;
+// 0x1.74910d52d3052p62
+// > printdouble(y);
+// 0x43d74910d52d3052
+constexpr uint64_t Y_UPPER_BOUND = 0x43d7'4910'd52d'3052;
+// Lower bound for y = 2^(-54) / 1074), generated by Sollya:
+// > y = round(2^(-54) / 1074), D, RD);
+// > y;
+// 0x1.e829f39aef509p-65
+// > printdouble(y);
+// 0x3bee829f39ae'f509
+constexpr uint64_t Y_LOWER_BOUND = 0x3bee'829f'39ae'f509;
+
+// Fast checks for special inputs:
+// x = 0, +-1, 2^k, 10, +- inf
+// y = 0, +-1, 2, 0.5
+LIBC_ALWAYS_INLINE cpp::optional<double> check_special_inputs(double x,
+ double y) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x), ybits(y);
+
+ bool x_sign = xbits.sign() == Sign::NEG;
+ bool y_sign = ybits.sign() == Sign::NEG;
+
+ FPBits x_abs = xbits.abs();
+ FPBits y_abs = ybits.abs();
+
+ uint64_t x_u = xbits.uintval();
+ uint64_t y_u = ybits.uintval();
+ uint64_t y_a = y_abs.uintval();
+
+ if (LIBC_UNLIKELY((x_u & 0x0003'FFFF'FFFF'FFFF) == 0) ||
+ ((y_u & 0x000F'FFFF'FFFF'FFFF) == 0)) {
+ // If x or y is signaling NaN
+ if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ if (x == 1.0 || y == 0.0)
+ return 1.0;
+
+ if (x == 0.0) {
+ if (y_abs.is_nan())
+ return y;
+ if (y_abs.is_inf())
+ return y_sign ? FPBits::inf().get_val() : 0.0;
+ bool out_is_neg = x_sign && is_odd_integer(y);
+ if (y_sign) {
+ // pow(0, negative number) = inf
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return FPBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
+ }
+ // pow(0, positive number) = 0
+ return out_is_neg ? -0.0 : 0.0;
+ }
+
+ if (y == 1.0)
+ return x;
+
+ if (y == 2.0)
+ return x * x;
+
+ if (y == 0.5 && !x_sign)
+ return fputil::sqrt<double>(x);
+
+ // TODO: Add special case y = -0.5 when rsqrt is available.
+
+ if (x == 2.0)
+ return math::exp2(y);
+
+ if (x == 10.0)
+ return math::exp10(y);
+
+ // For x = 2^(+- 2^n):
+ // x^y = 2^(+- 2^n * y).
+ if (!x_sign && x_abs.is_normal() && xbits.get_mantissa() == 0 &&
+ y_a > Y_LOWER_BOUND && y_a < Y_UPPER_BOUND) {
+ int e_x = xbits.get_exponent();
+ uint32_t abs_e = static_cast<uint32_t>(e_x > 0 ? e_x : -e_x);
+ if (cpp::has_single_bit(abs_e)) {
+ double hi = static_cast<double>(e_x) * y;
+ if (hi >= 1024.0)
+ return set_overflow(false);
+ if (hi < -1075.0)
+ return set_underflow(false);
+ if (hi >= -1022.0)
+ return math::exp2(hi);
+ }
+ }
+ }
+
+ return cpp::nullopt;
+}
+
+// Filters out extreme input ranges, infinities, NaNs, normalizes denormal
+// inputs, and handles negative bases. Returns a value if the result is
+// determined, or nullopt if regular evaluation should proceed (in which case x,
+// y, e_x, x_mant, is_neg, sign_d may be updated).
+LIBC_ALWAYS_INLINE cpp::optional<double>
+check_exceptional_cases(double &x, double &y, double &e_x, uint64_t &x_mant,
+ bool &is_neg, double &sign_d) {
+ using FPBits = fputil::FPBits<double>;
+
+ FPBits xbits(x), ybits(y);
+ bool x_sign = xbits.sign() == Sign::NEG;
+ bool y_sign = ybits.sign() == Sign::NEG;
+
+ FPBits x_abs = xbits.abs();
+ FPBits y_abs = ybits.abs();
+
+ uint64_t x_u = xbits.uintval();
+ uint64_t y_mant = ybits.get_mantissa();
+ uint64_t x_a = x_abs.uintval();
+ uint64_t y_a = y_abs.uintval();
+
+ if (LIBC_UNLIKELY(y_a <= Y_LOWER_BOUND || y_a >= Y_UPPER_BOUND ||
+ x_u >= FPBits::inf().uintval() ||
+ x_u < FPBits::min_normal().uintval())) {
+ // If x or y is signaling NaN
+ if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // |y| > |1075 / log2(1 - 2^-53)|.
+ if (y_a > 0x43d7'4910'd52d'3052) {
+ if (y_a >= 0x7ff0'0000'0000'0000) {
+ // y is inf or nan
+ if (y_mant != 0) {
+ // y is NaN
+ // pow(1, NaN) = 1 is already dealt with.
+ // pow(x, NaN) = NaN
+ return y;
+ }
+
+ // Now y is +-Inf
+ if (x_abs.is_nan()) {
+ // pow(NaN, +-Inf) = NaN
+ return x;
+ }
+
+ if (x_a == 0x3ff0'0000'0000'0000) {
+ // pow(+-1, +-Inf) = 1.0
+ return 1.0;
+ }
+
+ // pow (|x| < 1, -inf) = +inf
+ // pow (|x| < 1, +inf) = 0.0
+ // pow (|x| > 1, -inf) = 0.0
+ // pow (|x| > 1, +inf) = +inf
+ return ((x_a < FPBits::one().uintval()) == y_sign)
+ ? FPBits::inf().get_val()
+ : 0.0;
+ }
+ // x^y will overflow / underflow in double precision. Set y to a
+ // large enough exponent but not too large, so that the computations
+ // won't overflow in double precision.
+ y = y_sign ? -0x1.0p100 : 0x1.0p100;
+ }
+
+ // y is finite and non-zero.
+
+ if (x_a == FPBits::inf().uintval()) {
+ bool out_is_neg = x_sign && is_odd_integer(y);
+ Sign out_sign = out_is_neg ? Sign::NEG : Sign::POS;
+ return y_sign ? FPBits::zero(out_sign).get_val()
+ : FPBits::inf(out_sign).get_val();
+ }
+
+ if (x_a > FPBits::inf().uintval()) {
+ // x is NaN.
+ // pow (aNaN, 0) is already taken care above.
+ return x;
+ }
+
+ // Normalize denormal inputs.
+ if (x_a < FPBits::min_normal().uintval()) {
+ FPBits x_norm(x * 0x1.0p64);
+ e_x = static_cast<double>(x_norm.get_exponent()) - 64.0;
+ x_mant = x_norm.get_mantissa();
+ }
+
+ // x is finite and negative, and y is a finite integer.
+ if (x_sign) {
+ if (is_integer(y)) {
+ x = -x;
+ if (is_odd_integer(y)) {
+ is_neg = true;
+ sign_d = -1.0;
+ }
+ } else {
+ // pow( negative, non-integer ) = NaN
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ }
+
+ if (y_a <= Y_LOWER_BOUND) {
+ volatile double one = 1.0;
+ volatile double eps = ((x_a < FPBits::one().uintval()) == !y_sign)
+ ? -0x1.0p-100
+ : 0x1.0p-100;
+ return one + eps;
+ }
+ }
+
+ return cpp::nullopt;
+}
+
+// Check if x^y is an exact rounding boundary case (a 54-bit dyadic float,
+// which is either an exact 53-bit float or an exact halfway midpoint).
+//
+// Reference:
+// Lauter, C. and Lefevre, V., "Rounding Boundary Cases for Values of the
+// Exponential and Power Functions," IEEE Trans. Comput. 58(8):1063-1074.
+//
+// For a normalized 128-bit dyadic float with MSB at bit 127, the top 54 bits
+// occupy bits [127:74]. For any 54-bit number, bits [73:0] are zero.
+//
+// If x^y is not in F_54:
+// |x^y - o_54(x^y)| / x^y >= 2^-114.
+// With an evaluation error < 2^-115, the distance from bits [73:0] to the
+// nearest 54-bit boundary must satisfy:
+// dist >= (2^-114 - 2^-115) * 2^127 = 2^12.
+// We use a threshold of 2^-114 (dist < 2^13) to capture all boundary
+// candidates:
+//
+// 1. x = 2^e: x^y = 2^(e * y) is in F_54 iff e * y is an integer.
+// 2. x = m * 2^e, y = n * 2^f with m, n odd integers:
+// x^y in F_54 implies:
+// - 0 <= y <= 35, n <= 35, and f >= -5.
+// - e * y is an integer.
+// - When f < 0 (i.e. y = n / 2^(-f)): m^(2^f) is an integer, verified
+// by taking repeated square roots -f times.
+// - The resulting significand (m^(2^f))^n <= 2^54.
+//
+// Returns true and sets exact_m and exact_exp if x^y is an exact boundary.
+LIBC_INLINE bool is_exact_rounding_boundary(double x, double y,
+ uint64_t &exact_m, int &exact_exp) {
+ using FPBits = fputil::FPBits<double>;
+
+ FPBits xbits(x);
+ int x_e = 0;
+ uint64_t x_mant = xbits.get_mantissa();
+ if (LIBC_UNLIKELY(xbits.get_biased_exponent() == 0)) {
+ if (x_mant != 0) {
+ int shift = cpp::countl_zero(x_mant) - 11;
+ x_mant = (x_mant << shift) & FPBits::FRACTION_MASK;
+ x_e = -1022 - shift;
+ }
+ } else {
+ x_e = xbits.get_exponent();
+ }
+
+ // Case 1: x = 2^e.
+ if (x_mant == 0) {
+ double e = static_cast<double>(x_e);
+ double ey = e * y;
+ if (is_integer(ey)) {
+ exact_m = 1;
+ exact_exp = static_cast<int>(ey);
+ return true;
+ }
+ return false;
+ }
+
+ // Case 2: x is not a power of 2.
+ if (y < 0.0 || y > 35.0)
+ return false;
+
+ // Decompose y = n * 2^f with n odd.
+ FPBits ybits(y);
+ uint64_t y_mant = ybits.get_mantissa() | (1ULL << FPBits::FRACTION_LEN);
+ int y_exp = ybits.get_exponent() - static_cast<int>(FPBits::FRACTION_LEN);
+ int tz_y = cpp::countr_zero(y_mant);
+ uint64_t n = y_mant >> tz_y;
+ int f = y_exp + tz_y;
+
+ if (n > 35 || f < -5)
+ return false;
+
+ // Decompose x = m * 2^e with m odd.
+ uint64_t full_x_mant = x_mant | (1ULL << FPBits::FRACTION_LEN);
+ int tz_x = cpp::countr_zero(full_x_mant);
+ uint64_t m = full_x_mant >> tz_x;
+ int e = x_e - static_cast<int>(FPBits::FRACTION_LEN) + tz_x;
+
+ if (f < 0) {
+ // Non-integer exponent: y = n / 2^(-f) with -f in {1, ..., 5}.
+ double ey = static_cast<double>(e) * y;
+ if (!is_integer(ey))
+ return false;
+
+ // Check if m^(2^f) is an integer by taking repeated square roots -f times.
+ int count = -f;
+ uint64_t cur = m;
+ for (int i = 0; i < count; ++i) {
+ uint64_t s =
+ static_cast<uint64_t>(fputil::sqrt<double>(static_cast<double>(cur)));
+ if (s * s != cur)
+ return false;
+ cur = s;
+ }
+
+ // Compute res = cur^n and ensure res <= 2^54.
+ uint64_t res = 1;
+ for (uint64_t i = 0; i < n; ++i) {
+ if (res > (1ULL << 54) / cur)
+ return false;
+ res *= cur;
+ }
+ exact_m = res;
+ exact_exp = static_cast<int>(ey);
+ return true;
+ } else {
+ // Integer exponent: f >= 0, so y is an integer.
+ int y_int = static_cast<int>(y);
+ uint64_t res = 1;
+ for (int i = 0; i < y_int; ++i) {
+ if (res > (1ULL << 54) / m)
+ return false;
+ res *= m;
+ }
+ exact_m = res;
+ exact_exp = e * y_int;
+ return true;
+ }
+}
+
+} // namespace pow_internal
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_UTILS_H
diff --git a/libc/test/src/math/pow_test.cpp b/libc/test/src/math/pow_test.cpp
index 20e3ddfc8fbe5..d29451c13f39e 100644
--- a/libc/test/src/math/pow_test.cpp
+++ b/libc/test/src/math/pow_test.cpp
@@ -1,16 +1,28 @@
-//===-- Unittests for pow -------------------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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
+/// Unittests for pow.
+///
+//===----------------------------------------------------------------------===//
+#include "src/__support/macros/optimization.h"
#include "src/math/pow.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#define TOLERANCE 1
+#else
+#define TOLERANCE 0
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
using LlvmLibcPowTest = LIBC_NAMESPACE::testing::FPTest<double>;
using LIBC_NAMESPACE::testing::tlog;
@@ -25,14 +37,18 @@ TEST_F(LlvmLibcPowTest, TrickyInputs) {
{0x1.ffffffffffffcp-1, 0x1.fffffffffffffp-2},
{0x1.f558a88a8aadep-1, 0x1.88ap+12},
{0x1.e84d32731e593p-1, 0x1.2cb8p+13},
- {0x1.ffffffffffffcp-1, 0x1.fffffffffffffp-2},
+ {0x1.fd98b527935a2p-1, 0x1.71bp+13},
+ {0x1.01ff9791c39f9p+0, 0x1.62cc9e927e14ap+16},
+ {0x1.01ff3ca3af38bp+0, -0x1.605d9e6a8dce2p+16},
+ {0x1.01ff42c5cad93p+0, 0x1.63b5bb6d49873p+16},
+ {0x1.01e6666666669p+0, 0x1.fffffffffffffp+15},
};
for (auto input : INPUTS) {
double x = input.x;
double y = input.y;
EXPECT_MPFR_MATCH(mpfr::Operation::Pow, input, LIBC_NAMESPACE::pow(x, y),
- 1.5);
+ TOLERANCE + 0.5);
}
}
@@ -77,7 +93,8 @@ TEST_F(LlvmLibcPowTest, InFloatRange) {
mpfr::BinaryInput<double> inputs{x, y};
if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Pow, inputs,
- result, 1.5, rounding_mode)) {
+ result, TOLERANCE + 0.5,
+ rounding_mode)) {
++fails;
while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
mpfr::Operation::Pow, inputs, result, tol, rounding_mode)) {
@@ -100,13 +117,15 @@ TEST_F(LlvmLibcPowTest, InFloatRange) {
}
if (fails) {
mpfr::BinaryInput<double> inputs{mx, my};
- EXPECT_MPFR_MATCH(mpfr::Operation::Pow, inputs, mr, 1.5, rounding_mode);
+ EXPECT_MPFR_MATCH(mpfr::Operation::Pow, inputs, mr, TOLERANCE + 0.5,
+ rounding_mode);
}
};
tlog << " Test Rounding To Nearest...\n";
test(mpfr::RoundingMode::Nearest);
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
tlog << " Test Rounding Downward...\n";
test(mpfr::RoundingMode::Downward);
@@ -115,4 +134,5 @@ TEST_F(LlvmLibcPowTest, InFloatRange) {
tlog << " Test Rounding Toward Zero...\n";
test(mpfr::RoundingMode::TowardZero);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}
diff --git a/libc/test/src/math/smoke/pow_test.cpp b/libc/test/src/math/smoke/pow_test.cpp
index ab8a0cf75294b..9e7599d9546d9 100644
--- a/libc/test/src/math/smoke/pow_test.cpp
+++ b/libc/test/src/math/smoke/pow_test.cpp
@@ -1,10 +1,15 @@
-//===-- Unittests for pow -------------------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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
+/// Smoke tests for pow.
+///
+//===----------------------------------------------------------------------===//
#include "hdr/fenv_macros.h"
#include "src/math/pow.h"
@@ -69,8 +74,7 @@ TEST_F(LlvmLibcPowTest, SpecialNumbers) {
EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::pow(zero, zero));
EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::pow(zero, neg_zero));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::pow(zero, inf));
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::pow(zero, neg_inf),
- FE_DIVBYZERO);
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::pow(zero, neg_inf));
EXPECT_FP_IS_NAN(LIBC_NAMESPACE::pow(zero, aNaN));
// pow( -0.0, exponent )
@@ -89,8 +93,7 @@ TEST_F(LlvmLibcPowTest, SpecialNumbers) {
EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::pow(neg_zero, zero));
EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::pow(neg_zero, neg_zero));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::pow(neg_zero, inf));
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::pow(neg_zero, neg_inf),
- FE_DIVBYZERO);
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::pow(neg_zero, neg_inf));
EXPECT_FP_IS_NAN(LIBC_NAMESPACE::pow(neg_zero, aNaN));
// pow( 1.0, exponent )
@@ -195,26 +198,24 @@ TEST_F(LlvmLibcPowTest, SpecialNumbers) {
EXPECT_FP_EQ(zero, LIBC_NAMESPACE::pow(-1.1, neg_inf));
// Exact powers of 2:
- // TODO: Enable these tests when we use exp2.
- // EXPECT_FP_EQ(0x1.0p15, LIBC_NAMESPACE::pow(2.0, 15.0));
- // EXPECT_FP_EQ(0x1.0p126, LIBC_NAMESPACE::pow(2.0, 126.0));
- // EXPECT_FP_EQ(0x1.0p-45, LIBC_NAMESPACE::pow(2.0, -45.0));
- // EXPECT_FP_EQ(0x1.0p-126, LIBC_NAMESPACE::pow(2.0, -126.0));
- // EXPECT_FP_EQ(0x1.0p-149, LIBC_NAMESPACE::pow(2.0, -149.0));
+ EXPECT_FP_EQ(0x1.0p15, LIBC_NAMESPACE::pow(2.0, 15.0));
+ EXPECT_FP_EQ(0x1.0p126, LIBC_NAMESPACE::pow(2.0, 126.0));
+ EXPECT_FP_EQ(0x1.0p-45, LIBC_NAMESPACE::pow(2.0, -45.0));
+ EXPECT_FP_EQ(0x1.0p-126, LIBC_NAMESPACE::pow(2.0, -126.0));
+ EXPECT_FP_EQ(0x1.0p-149, LIBC_NAMESPACE::pow(2.0, -149.0));
// Exact powers of 10:
- // TODO: Enable these tests when we use exp10
- // EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::pow(10.0, 0.0));
- // EXPECT_FP_EQ(10.0, LIBC_NAMESPACE::pow(10.0, 1.0));
- // EXPECT_FP_EQ(100.0, LIBC_NAMESPACE::pow(10.0, 2.0));
- // EXPECT_FP_EQ(1000.0, LIBC_NAMESPACE::pow(10.0, 3.0));
- // EXPECT_FP_EQ(10000.0, LIBC_NAMESPACE::pow(10.0, 4.0));
- // EXPECT_FP_EQ(100000.0, LIBC_NAMESPACE::pow(10.0, 5.0));
- // EXPECT_FP_EQ(1000000.0, LIBC_NAMESPACE::pow(10.0, 6.0));
- // EXPECT_FP_EQ(10000000.0, LIBC_NAMESPACE::pow(10.0, 7.0));
- // EXPECT_FP_EQ(100000000.0, LIBC_NAMESPACE::pow(10.0, 8.0));
- // EXPECT_FP_EQ(1000000000.0, LIBC_NAMESPACE::pow(10.0, 9.0));
- // EXPECT_FP_EQ(10000000000.0, LIBC_NAMESPACE::pow(10.0, 10.0));
+ EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::pow(10.0, 0.0));
+ EXPECT_FP_EQ(10.0, LIBC_NAMESPACE::pow(10.0, 1.0));
+ EXPECT_FP_EQ(100.0, LIBC_NAMESPACE::pow(10.0, 2.0));
+ EXPECT_FP_EQ(1000.0, LIBC_NAMESPACE::pow(10.0, 3.0));
+ EXPECT_FP_EQ(10000.0, LIBC_NAMESPACE::pow(10.0, 4.0));
+ EXPECT_FP_EQ(100000.0, LIBC_NAMESPACE::pow(10.0, 5.0));
+ EXPECT_FP_EQ(1000000.0, LIBC_NAMESPACE::pow(10.0, 6.0));
+ EXPECT_FP_EQ(10000000.0, LIBC_NAMESPACE::pow(10.0, 7.0));
+ EXPECT_FP_EQ(100000000.0, LIBC_NAMESPACE::pow(10.0, 8.0));
+ EXPECT_FP_EQ(1000000000.0, LIBC_NAMESPACE::pow(10.0, 9.0));
+ EXPECT_FP_EQ(10000000000.0, LIBC_NAMESPACE::pow(10.0, 10.0));
// Overflow / Underflow:
if (ROUNDING_MODES[i] != RoundingMode::Downward &&
>From cf636897161ed1f74255f480126f08e5ef0e9e6f Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 11 Sep 2026 02:32:01 +0000
Subject: [PATCH 2/2] Move FPBits definition to where it was used inside
pow_utils.h.
---
libc/src/__support/math/pow_utils.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/__support/math/pow_utils.h b/libc/src/__support/math/pow_utils.h
index a5d026737168e..900e9a39205b8 100644
--- a/libc/src/__support/math/pow_utils.h
+++ b/libc/src/__support/math/pow_utils.h
@@ -242,9 +242,9 @@ LIBC_INLINE double set_underflow(bool is_neg) {
volatile double x = is_neg ? -0x1.0p-1022 : 0x1.0p-1022;
return x * 0x1.0p-100;
#else
- using FPBits = fputil::FPBits<double>;
fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ using FPBits = fputil::FPBits<double>;
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD && !is_neg)
return FPBits::min_subnormal().get_val();
More information about the libc-commits
mailing list