[libc-commits] [libc] [llvm] [libc][math] Fix powf exact cases, improve its performance, and add code-size optimized float-only version. (PR #226588)
via libc-commits
libc-commits at lists.llvm.org
Fri Sep 25 13:50:56 PDT 2026
https://github.com/lntue created https://github.com/llvm/llvm-project/pull/226588
Powf improvement:
- Fix exact case detection
- Improve performance
- Add code-size optimized float-only option
- Reorganize the implementation selection and tests.
Assisted-by: Gemini is used for code-size and performance analysis, and for refactoring.
>From 24fda5df518a4f78077695172404038e8b7a384b Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 25 Sep 2026 20:44:58 +0000
Subject: [PATCH] [libc][math] Fix powf exact cases, improve its performance,
and add code-size optimized float-only version.
---
libc/src/__support/FPUtil/double_double.h | 16 +-
libc/src/__support/math/CMakeLists.txt | 70 +-
libc/src/__support/math/powf.h | 1047 +----------------
libc/src/__support/math/powf_double_eval.h | 548 +++++++++
libc/src/__support/math/powf_float_eval.h | 263 +++++
libc/src/__support/math/powf_small_tables.h | 120 --
libc/src/__support/math/powf_utils.h | 718 +++++++++++
libc/test/src/math/CMakeLists.txt | 2 +
libc/test/src/math/powf_test.cpp | 239 ++--
libc/test/src/math/smoke/CMakeLists.txt | 2 +
libc/test/src/math/smoke/powf_test.cpp | 492 ++++----
.../llvm-project-overlay/libc/BUILD.bazel | 71 +-
12 files changed, 2085 insertions(+), 1503 deletions(-)
create mode 100644 libc/src/__support/math/powf_double_eval.h
create mode 100644 libc/src/__support/math/powf_float_eval.h
delete mode 100644 libc/src/__support/math/powf_small_tables.h
create mode 100644 libc/src/__support/math/powf_utils.h
diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h
index 87b3887716b0fd..6e41bb9b2a4bed 100644
--- a/libc/src/__support/FPUtil/double_double.h
+++ b/libc/src/__support/FPUtil/double_double.h
@@ -71,19 +71,21 @@ LIBC_INLINE constexpr NumberPair<T> exact_add(T a, T b) {
// If `lsb(a.hi) >= ulp(b.hi)`, then the errors:
// err = (a.hi + a.lo + b.hi + b.lo) - (r.hi - r.lo) is bounded by:
// |err| < 2^(-2p) * ufp(r.hi) + 2^(-p + 2) * ufp(a.lo + b.lo).
-template <typename T>
+template <bool FAST2SUM = true, typename T>
LIBC_INLINE constexpr NumberPair<T> add(const NumberPair<T> &a,
const NumberPair<T> &b) {
- NumberPair<T> r = exact_add(a.hi, b.hi);
+ NumberPair<T> r = exact_add<FAST2SUM>(a.hi, b.hi);
T lo = a.lo + b.lo;
- return exact_add(r.hi, r.lo + lo);
+ T r_lo = r.lo + lo;
+ return exact_add<FAST2SUM>(r.hi, r_lo);
}
-// Assumption: |a.hi| >= |b|
-template <typename T>
+// Assumption: when FAST2SUM = true, |a.hi| >= |b|
+template <bool FAST2SUM = true, typename T>
LIBC_INLINE constexpr NumberPair<T> add(const NumberPair<T> &a, T b) {
- NumberPair<T> r = exact_add<false>(a.hi, b);
- return exact_add(r.hi, r.lo + a.lo);
+ NumberPair<T> r = exact_add<FAST2SUM>(a.hi, b);
+ T r_lo = r.lo + a.lo;
+ return exact_add<FAST2SUM>(r.hi, r_lo);
}
// Veltkamp's Splitting for double precision.
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 8d6a1495e1c110..7f6067fc7a2ff4 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -5402,24 +5402,84 @@ add_header_library(
)
add_header_library(
- powf
+ powf_utils
HDRS
- powf.h
- powf_small_tables.h
+ powf_utils.h
DEPENDS
.common_constants
.exp10f
.exp2f
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.CPP.bit
+ libc.src.__support.CPP.optional
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.FPUtil.rounding_mode
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.FPUtil.triple_double
+ libc.src.__support.common
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
+add_header_library(
+ powf_double_eval
+ HDRS
+ powf_double_eval.h
+ DEPENDS
+ .common_constants
+ .exp_constants
+ .powf_utils
libc.src.__support.CPP.algorithm
libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.manipulation_functions
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
libc.src.__support.FPUtil.polyeval
- libc.src.__support.FPUtil.sqrt
- libc.src.__support.FPUtil.triple_double
+ libc.src.__support.common
+ libc.src.__support.macros.config
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+)
+
+add_header_library(
+ powf_float_eval
+ HDRS
+ powf_float_eval.h
+ DEPENDS
+ .exp2f_float_utils
+ .powf_utils
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.double_double
+ 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.rounding_mode
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+)
+
+add_header_library(
+ powf
+ HDRS
+ powf.h
+ DEPENDS
+ .powf_double_eval
+ .powf_float_eval
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
)
add_header_library(
diff --git a/libc/src/__support/math/powf.h b/libc/src/__support/math/powf.h
index 944b477ac44612..2671215d7d7c76 100644
--- a/libc/src/__support/math/powf.h
+++ b/libc/src/__support/math/powf.h
@@ -14,1051 +14,32 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_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"
#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \
- defined(LIBC_MATH_HAS_SMALL_TABLES)
-
-#include "src/__support/math/powf_small_tables.h"
+ (defined(LIBC_MATH_HAS_SMALL_TABLES) || \
+ defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT))
-#else
+#include "src/__support/math/powf_float_eval.h"
+#define LIBC_MATH_POWF_IMPL float_eval
-#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
-#include "exp10f.h" // Speedup for powf(10, y) = exp10f(y)
-#include "exp2f.h" // Speedup for powf(2, y) = exp2f(y)
-#include "exp_constants.h"
+#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#include "src/__support/math/powf_double_eval.h"
+#define LIBC_MATH_POWF_IMPL double_eval
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS && LIBC_MATH_HAS_SMALL_TABLES
-
-#include "src/__support/CPP/algorithm.h"
-#include "src/__support/CPP/bit.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 powf(x, 1/2) = sqrtf(x)
-#include "src/__support/FPUtil/triple_double.h"
-#include "src/__support/common.h"
-#include "src/__support/macros/config.h"
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
namespace LIBC_NAMESPACE_DECL {
-
namespace math {
-namespace powf_internal {
-
-using fputil::DoubleDouble;
-using fputil::TripleDouble;
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-alignas(16) LIBC_INLINE_VAR constexpr DoubleDouble LOG2_R_DD[128] = {
- {0.0, 0.0},
- {-0x1.177c23362928cp-25, 0x1.72c8p-7},
- {-0x1.179e0caa9c9abp-22, 0x1.744p-6},
- {-0x1.c6cea541f5b7p-23, 0x1.184cp-5},
- {-0x1.66c4d4e554434p-22, 0x1.773ap-5},
- {-0x1.70700a00fdd55p-24, 0x1.d6ecp-5},
- {0x1.53002a4e86631p-23, 0x1.1bb3p-4},
- {0x1.fcd15f101c142p-25, 0x1.4c56p-4},
- {0x1.25b3eed319cedp-22, 0x1.7d6p-4},
- {-0x1.4195120d8486fp-22, 0x1.960dp-4},
- {0x1.45b878e27d0d9p-23, 0x1.c7b5p-4},
- {0x1.770744593a4cbp-22, 0x1.f9c9p-4},
- {0x1.c673032495d24p-22, 0x1.097ep-3},
- {-0x1.1eaa65b49696ep-22, 0x1.22dbp-3},
- {0x1.b2866f2850b22p-22, 0x1.3c6f8p-3},
- {0x1.8ee37cd2ea9d3p-25, 0x1.494f8p-3},
- {0x1.7e86f9c2154fbp-24, 0x1.633a8p-3},
- {0x1.8e3cfc25f0ce6p-26, 0x1.7046p-3},
- {0x1.57f7a64ccd537p-28, 0x1.8a898p-3},
- {-0x1.a761c09fbd2aep-22, 0x1.97c2p-3},
- {0x1.24bea9a2c66f3p-22, 0x1.b26p-3},
- {-0x1.60002ccfe43f5p-25, 0x1.bfc68p-3},
- {0x1.69f220e97f22cp-22, 0x1.dac2p-3},
- {-0x1.6164f64c210ep-22, 0x1.e858p-3},
- {-0x1.0c1678ae89767p-24, 0x1.01d9cp-2},
- {-0x1.f26a05c813d57p-22, 0x1.08bdp-2},
- {0x1.4d8fc561c8d44p-24, 0x1.169cp-2},
- {-0x1.362ad8f7ca2dp-22, 0x1.1d984p-2},
- {0x1.2b13cd6c4d042p-22, 0x1.249ccp-2},
- {-0x1.1c8f11979a5dbp-22, 0x1.32cp-2},
- {0x1.c2ab3edefe569p-23, 0x1.39de8p-2},
- {0x1.7c3eca28e69cap-26, 0x1.4106p-2},
- {-0x1.34c4e99e1c6c6p-24, 0x1.4f6fcp-2},
- {-0x1.194a871b63619p-22, 0x1.56b24p-2},
- {0x1.e3dd5c1c885aep-23, 0x1.5dfdcp-2},
- {-0x1.6ccf3b1129b7cp-23, 0x1.6552cp-2},
- {-0x1.2f346e2bf924bp-23, 0x1.6cb1p-2},
- {-0x1.fa61aaa59c1d8p-23, 0x1.7b8ap-2},
- {0x1.90c11fd32a3abp-22, 0x1.8304cp-2},
- {0x1.57f7a64ccd537p-27, 0x1.8a898p-2},
- {0x1.249ba76fee235p-27, 0x1.9218p-2},
- {-0x1.aad2729b21ae5p-23, 0x1.99b08p-2},
- {0x1.71810a5e1818p-22, 0x1.a8ff8p-2},
- {-0x1.6172fe015e13cp-27, 0x1.b0b68p-2},
- {0x1.5ec6c1bfbf89ap-24, 0x1.b877cp-2},
- {0x1.678bf6cdedf51p-24, 0x1.c0438p-2},
- {0x1.c2d45fe43895ep-22, 0x1.c819cp-2},
- {-0x1.9ee52ed49d71dp-22, 0x1.cffbp-2},
- {0x1.5786af187a96bp-27, 0x1.d7e6cp-2},
- {0x1.3ab0dc56138c9p-23, 0x1.dfdd8p-2},
- {0x1.fe538ab34efb5p-22, 0x1.e7df4p-2},
- {-0x1.e4fee07aa4b68p-22, 0x1.efec8p-2},
- {-0x1.172f32fe67287p-22, 0x1.f804cp-2},
- {-0x1.9a83ff9ab9cc8p-22, 0x1.00144p-1},
- {-0x1.68cb06cece193p-22, 0x1.042bep-1},
- {0x1.8cd71ddf82e2p-22, 0x1.08494p-1},
- {0x1.5e18ab2df3ae6p-22, 0x1.0c6cap-1},
- {0x1.5dee4d9d8a273p-25, 0x1.1096p-1},
- {0x1.fcd15f101c142p-26, 0x1.14c56p-1},
- {-0x1.2474b0f992ba1p-23, 0x1.18faep-1},
- {0x1.4b5a92a606047p-24, 0x1.1d368p-1},
- {0x1.16186fcf54bbdp-22, 0x1.21786p-1},
- {0x1.18efabeb7d722p-27, 0x1.25c0ap-1},
- {-0x1.e5fc7d238691dp-24, 0x1.2a0f4p-1},
- {0x1.f5809faf6283cp-22, 0x1.2e644p-1},
- {0x1.f5809faf6283cp-22, 0x1.2e644p-1},
- {0x1.c6e1dcd0cb449p-22, 0x1.32bfep-1},
- {0x1.76e0e8f74b4d5p-22, 0x1.37222p-1},
- {-0x1.cb82c89692d99p-24, 0x1.3b8b2p-1},
- {-0x1.63161c5432aebp-22, 0x1.3ffaep-1},
- {0x1.458104c41b901p-22, 0x1.44716p-1},
- {0x1.458104c41b901p-22, 0x1.44716p-1},
- {-0x1.cd9d0cde578d5p-22, 0x1.48efp-1},
- {0x1.b9884591add87p-26, 0x1.4d738p-1},
- {0x1.c6042978605ffp-22, 0x1.51ff2p-1},
- {-0x1.fc4c96b37dcf6p-22, 0x1.56922p-1},
- {-0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
- {-0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
- {0x1.c4e4fbb68a4d1p-22, 0x1.5fcdcp-1},
- {-0x1.9d499bd9b3226p-23, 0x1.6476ep-1},
- {-0x1.f89b355ede26fp-23, 0x1.69278p-1},
- {-0x1.f89b355ede26fp-23, 0x1.69278p-1},
- {0x1.53c7e319f6e92p-24, 0x1.6ddfcp-1},
- {-0x1.b291f070528c7p-22, 0x1.729fep-1},
- {0x1.2967a451a7b48p-25, 0x1.7767cp-1},
- {0x1.2967a451a7b48p-25, 0x1.7767cp-1},
- {0x1.244fcff690fcep-22, 0x1.7c37ap-1},
- {0x1.46fd97f5dc572p-23, 0x1.810fap-1},
- {0x1.46fd97f5dc572p-23, 0x1.810fap-1},
- {-0x1.f3a7352663e5p-22, 0x1.85efep-1},
- {0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
- {0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
- {0x1.3226b211bf1d9p-23, 0x1.8fc92p-1},
- {0x1.d24b136c101eep-23, 0x1.94c28p-1},
- {0x1.d24b136c101eep-23, 0x1.94c28p-1},
- {0x1.7c40c7907e82ap-22, 0x1.99c48p-1},
- {-0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
- {-0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
- {-0x1.6a77813f94e01p-22, 0x1.a3e3p-1},
- {-0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
- {-0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
- {-0x1.f983f74d3138fp-23, 0x1.ae256p-1},
- {-0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
- {-0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
- {-0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
- {-0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
- {-0x1.19b4f3c72c4f8p-24, 0x1.bdceap-1},
- {0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
- {0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
- {-0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
- {-0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
- {-0x1.6e46335aae723p-24, 0x1.cdcecp-1},
- {-0x1.beb244c59f331p-22, 0x1.d3382p-1},
- {-0x1.beb244c59f331p-22, 0x1.d3382p-1},
- {0x1.16c071e93fd97p-27, 0x1.d8abap-1},
- {0x1.16c071e93fd97p-27, 0x1.d8abap-1},
- {0x1.d8175819530c2p-22, 0x1.de298p-1},
- {0x1.d8175819530c2p-22, 0x1.de298p-1},
- {0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
- {0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
- {0x1.914e204f19d94p-22, 0x1.e9452p-1},
- {0x1.914e204f19d94p-22, 0x1.e9452p-1},
- {0x1.c55d997da24fdp-22, 0x1.eee32p-1},
- {0x1.c55d997da24fdp-22, 0x1.eee32p-1},
- {-0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
- {-0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
- {0x1.7a4887bd74039p-22, 0x1.fa406p-1},
- {0.0, 1.0},
-};
-
-#else
-
-#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-LIBC_INLINE_VAR constexpr uint64_t ERR = 64;
-#else
-LIBC_INLINE_VAR constexpr uint64_t ERR = 128;
-#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-
-// We choose the precision of the high part to be 53 - 24 - 8, so that when
-// y * (e_x + LOG2_R_DD[i].hi) is exact.
-// Generated by Sollya with:
-// > for i from 0 to 127 do {
-// r = 2^-8 * ceil(2^8 * (1 - 2^-8) / (1 + i * 2^-7) );
-// a = -log2(r);
-// b = round(1 + a, 53 - 24 - 8, RN) - 1;
-// c = round(a - b, D, RN);
-// d = round(a - b - c, D, RN);
-// print("{", d, ",", c, ", ", b, "},");
-// };
-LIBC_INLINE_VAR constexpr TripleDouble LOG2_R_TD[128] = {
- {0.0, 0.0, 0.0},
- {0x1.84a2c615b70adp-79, -0x1.177c23362928cp-25, 0x1.72c8p-7},
- {-0x1.f27b820fd03eap-76, -0x1.179e0caa9c9abp-22, 0x1.744p-6},
- {-0x1.f27ef487c8f34p-77, -0x1.c6cea541f5b7p-23, 0x1.184cp-5},
- {-0x1.e3f80fbc71454p-76, -0x1.66c4d4e554434p-22, 0x1.773ap-5},
- {-0x1.9f8ef14d5f6eep-79, -0x1.70700a00fdd55p-24, 0x1.d6ecp-5},
- {0x1.452bbce7398c1p-77, 0x1.53002a4e86631p-23, 0x1.1bb3p-4},
- {-0x1.990555535afdp-81, 0x1.fcd15f101c142p-25, 0x1.4c56p-4},
- {0x1.447e30ad393eep-78, 0x1.25b3eed319cedp-22, 0x1.7d6p-4},
- {0x1.b7759da88a2dap-76, -0x1.4195120d8486fp-22, 0x1.960dp-4},
- {0x1.cee7766ece702p-78, 0x1.45b878e27d0d9p-23, 0x1.c7b5p-4},
- {-0x1.a55c745ecdc2fp-77, 0x1.770744593a4cbp-22, 0x1.f9c9p-4},
- {0x1.f7ec992caa67fp-77, 0x1.c673032495d24p-22, 0x1.097ep-3},
- {-0x1.433638c6ece3ep-77, -0x1.1eaa65b49696ep-22, 0x1.22dbp-3},
- {0x1.58f27b6518824p-76, 0x1.b2866f2850b22p-22, 0x1.3c6f8p-3},
- {-0x1.86bdcfdfd4a4cp-79, 0x1.8ee37cd2ea9d3p-25, 0x1.494f8p-3},
- {-0x1.ff7044a68a7fap-80, 0x1.7e86f9c2154fbp-24, 0x1.633a8p-3},
- {-0x1.aa21694561327p-81, 0x1.8e3cfc25f0ce6p-26, 0x1.7046p-3},
- {-0x1.d209f2d4239c6p-87, 0x1.57f7a64ccd537p-28, 0x1.8a898p-3},
- {-0x1.a55e97e60e632p-76, -0x1.a761c09fbd2aep-22, 0x1.97c2p-3},
- {0x1.261179225541ep-76, 0x1.24bea9a2c66f3p-22, 0x1.b26p-3},
- {-0x1.08fa30510fca9p-82, -0x1.60002ccfe43f5p-25, 0x1.bfc68p-3},
- {-0x1.63ec8d56242f9p-76, 0x1.69f220e97f22cp-22, 0x1.dac2p-3},
- {0x1.8bcdaf0534365p-76, -0x1.6164f64c210ep-22, 0x1.e858p-3},
- {0x1.1003282896056p-78, -0x1.0c1678ae89767p-24, 0x1.01d9cp-2},
- {0x1.01bcc7025fa92p-78, -0x1.f26a05c813d57p-22, 0x1.08bdp-2},
- {-0x1.fe8a8648e9ebcp-80, 0x1.4d8fc561c8d44p-24, 0x1.169cp-2},
- {0x1.08dfb23650c75p-79, -0x1.362ad8f7ca2dp-22, 0x1.1d984p-2},
- {-0x1.f8d5a89861a5ep-79, 0x1.2b13cd6c4d042p-22, 0x1.249ccp-2},
- {-0x1.a1c872983511ep-76, -0x1.1c8f11979a5dbp-22, 0x1.32cp-2},
- {0x1.e8e21bff3336bp-77, 0x1.c2ab3edefe569p-23, 0x1.39de8p-2},
- {0x1.fd1994fb2c4a1p-80, 0x1.7c3eca28e69cap-26, 0x1.4106p-2},
- {0x1.6b94b51cf76b1p-80, -0x1.34c4e99e1c6c6p-24, 0x1.4f6fcp-2},
- {-0x1.31d55da1d0f66p-76, -0x1.194a871b63619p-22, 0x1.56b24p-2},
- {-0x1.378b22691e28bp-77, 0x1.e3dd5c1c885aep-23, 0x1.5dfdcp-2},
- {0x1.99e302970e411p-83, -0x1.6ccf3b1129b7cp-23, 0x1.6552cp-2},
- {0x1.20164a049664dp-82, -0x1.2f346e2bf924bp-23, 0x1.6cb1p-2},
- {-0x1.d14aac4d864c3p-77, -0x1.fa61aaa59c1d8p-23, 0x1.7b8ap-2},
- {0x1.496ab4e4b293fp-79, 0x1.90c11fd32a3abp-22, 0x1.8304cp-2},
- {-0x1.d209f2d4239c6p-86, 0x1.57f7a64ccd537p-27, 0x1.8a898p-2},
- {0x1.eae3326327babp-81, 0x1.249ba76fee235p-27, 0x1.9218p-2},
- {0x1.fa05bddfded8cp-77, -0x1.aad2729b21ae5p-23, 0x1.99b08p-2},
- {-0x1.624140d175ba2p-77, 0x1.71810a5e1818p-22, 0x1.a8ff8p-2},
- {0x1.f1c5160c515c1p-81, -0x1.6172fe015e13cp-27, 0x1.b0b68p-2},
- {-0x1.86a6204eec8cp-79, 0x1.5ec6c1bfbf89ap-24, 0x1.b877cp-2},
- {0x1.718f761dd3915p-78, 0x1.678bf6cdedf51p-24, 0x1.c0438p-2},
- {-0x1.d4ee66c3700e4p-76, 0x1.c2d45fe43895ep-22, 0x1.c819cp-2},
- {-0x1.7d14533586306p-77, -0x1.9ee52ed49d71dp-22, 0x1.cffbp-2},
- {0x1.5ce9fb5a7bb5bp-81, 0x1.5786af187a96bp-27, 0x1.d7e6cp-2},
- {-0x1.ae6face57ad3bp-77, 0x1.3ab0dc56138c9p-23, 0x1.dfdd8p-2},
- {0x1.5ac93b443d55fp-78, 0x1.fe538ab34efb5p-22, 0x1.e7df4p-2},
- {0x1.f1753e0ae1e8fp-76, -0x1.e4fee07aa4b68p-22, 0x1.efec8p-2},
- {0x1.cdfd4c297069bp-76, -0x1.172f32fe67287p-22, 0x1.f804cp-2},
- {0x1.97a0e8f3ba742p-79, -0x1.9a83ff9ab9cc8p-22, 0x1.00144p-1},
- {-0x1.800450f5b2357p-78, -0x1.68cb06cece193p-22, 0x1.042bep-1},
- {-0x1.a839041241fe7p-78, 0x1.8cd71ddf82e2p-22, 0x1.08494p-1},
- {0x1.ed0b8eeccca86p-78, 0x1.5e18ab2df3ae6p-22, 0x1.0c6cap-1},
- {0x1.3dd41df9689b3p-79, 0x1.5dee4d9d8a273p-25, 0x1.1096p-1},
- {-0x1.990555535afdp-82, 0x1.fcd15f101c142p-26, 0x1.14c56p-1},
- {-0x1.1773d02c9055cp-77, -0x1.2474b0f992ba1p-23, 0x1.18faep-1},
- {-0x1.4aeef330c53c1p-78, 0x1.4b5a92a606047p-24, 0x1.1d368p-1},
- {0x1.8e6ff749ebacbp-77, 0x1.16186fcf54bbdp-22, 0x1.21786p-1},
- {0x1.c09d761c548ebp-84, 0x1.18efabeb7d722p-27, 0x1.25c0ap-1},
- {0x1.aaa73a428e1e4p-78, -0x1.e5fc7d238691dp-24, 0x1.2a0f4p-1},
- {-0x1.af2f3d8b63fbap-79, 0x1.f5809faf6283cp-22, 0x1.2e644p-1},
- {-0x1.af2f3d8b63fbap-79, 0x1.f5809faf6283cp-22, 0x1.2e644p-1},
- {0x1.78de359f2bb88p-77, 0x1.c6e1dcd0cb449p-22, 0x1.32bfep-1},
- {-0x1.415ae1a715618p-76, 0x1.76e0e8f74b4d5p-22, 0x1.37222p-1},
- {-0x1.4991b5375621fp-79, -0x1.cb82c89692d99p-24, 0x1.3b8b2p-1},
- {-0x1.827d37deb2236p-76, -0x1.63161c5432aebp-22, 0x1.3ffaep-1},
- {0x1.9576edac01c78p-77, 0x1.458104c41b901p-22, 0x1.44716p-1},
- {0x1.9576edac01c78p-77, 0x1.458104c41b901p-22, 0x1.44716p-1},
- {-0x1.05a27b81e2219p-77, -0x1.cd9d0cde578d5p-22, 0x1.48efp-1},
- {0x1.237616778b4bap-82, 0x1.b9884591add87p-26, 0x1.4d738p-1},
- {0x1.3b7d7e5d148bbp-76, 0x1.c6042978605ffp-22, 0x1.51ff2p-1},
- {-0x1.cc3f936a5977cp-79, -0x1.fc4c96b37dcf6p-22, 0x1.56922p-1},
- {0x1.20164a049664dp-83, -0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
- {0x1.20164a049664dp-83, -0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
- {-0x1.a212919a92f7ap-77, 0x1.c4e4fbb68a4d1p-22, 0x1.5fcdcp-1},
- {-0x1.b64b03f7230ddp-77, -0x1.9d499bd9b3226p-23, 0x1.6476ep-1},
- {-0x1.1ec6379e6e3b9p-77, -0x1.f89b355ede26fp-23, 0x1.69278p-1},
- {-0x1.1ec6379e6e3b9p-77, -0x1.f89b355ede26fp-23, 0x1.69278p-1},
- {-0x1.4ba44c03bfbbdp-78, 0x1.53c7e319f6e92p-24, 0x1.6ddfcp-1},
- {-0x1.c36fc650d030fp-77, -0x1.b291f070528c7p-22, 0x1.729fep-1},
- {-0x1.69e5693a7f067p-80, 0x1.2967a451a7b48p-25, 0x1.7767cp-1},
- {-0x1.69e5693a7f067p-80, 0x1.2967a451a7b48p-25, 0x1.7767cp-1},
- {0x1.6598aae91499ap-76, 0x1.244fcff690fcep-22, 0x1.7c37ap-1},
- {0x1.99d61ec432837p-77, 0x1.46fd97f5dc572p-23, 0x1.810fap-1},
- {0x1.99d61ec432837p-77, 0x1.46fd97f5dc572p-23, 0x1.810fap-1},
- {0x1.855c42078f81bp-76, -0x1.f3a7352663e5p-22, 0x1.85efep-1},
- {-0x1.59408e815107p-77, 0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
- {-0x1.59408e815107p-77, 0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
- {0x1.33b318085e50ap-78, 0x1.3226b211bf1d9p-23, 0x1.8fc92p-1},
- {0x1.343fe7c9cb4aep-79, 0x1.d24b136c101eep-23, 0x1.94c28p-1},
- {0x1.343fe7c9cb4aep-79, 0x1.d24b136c101eep-23, 0x1.94c28p-1},
- {-0x1.d19522e56fe6p-76, 0x1.7c40c7907e82ap-22, 0x1.99c48p-1},
- {-0x1.23b9d8ea55c3ep-77, -0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
- {-0x1.23b9d8ea55c3ep-77, -0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
- {0x1.829440c24aeb6p-78, -0x1.6a77813f94e01p-22, 0x1.a3e3p-1},
- {-0x1.624140d175ba2p-76, -0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
- {-0x1.624140d175ba2p-76, -0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
- {0x1.afa6f024fb045p-77, -0x1.f983f74d3138fp-23, 0x1.ae256p-1},
- {-0x1.603ad3a5d326dp-78, -0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
- {-0x1.603ad3a5d326dp-78, -0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
- {-0x1.0c1e0e5855d6ap-77, -0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
- {-0x1.0c1e0e5855d6ap-77, -0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
- {0x1.c817ad56baa16p-78, -0x1.19b4f3c72c4f8p-24, 0x1.bdceap-1},
- {0x1.44c47ac1bf62bp-77, 0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
- {0x1.44c47ac1bf62bp-77, 0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
- {-0x1.69b9465eae1e6p-78, -0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
- {-0x1.69b9465eae1e6p-78, -0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
- {-0x1.24a6d9d1d1904p-79, -0x1.6e46335aae723p-24, 0x1.cdcecp-1},
- {-0x1.3826144575ac4p-76, -0x1.beb244c59f331p-22, 0x1.d3382p-1},
- {-0x1.3826144575ac4p-76, -0x1.beb244c59f331p-22, 0x1.d3382p-1},
- {0x1.dbc96b3b12b25p-81, 0x1.16c071e93fd97p-27, 0x1.d8abap-1},
- {0x1.dbc96b3b12b25p-81, 0x1.16c071e93fd97p-27, 0x1.d8abap-1},
- {0x1.68a8ccdbd1f33p-77, 0x1.d8175819530c2p-22, 0x1.de298p-1},
- {0x1.68a8ccdbd1f33p-77, 0x1.d8175819530c2p-22, 0x1.de298p-1},
- {0x1.e586711df5ea1p-79, 0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
- {0x1.e586711df5ea1p-79, 0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
- {-0x1.bc25adf042483p-79, 0x1.914e204f19d94p-22, 0x1.e9452p-1},
- {-0x1.bc25adf042483p-79, 0x1.914e204f19d94p-22, 0x1.e9452p-1},
- {0x1.d7d82b65c5686p-76, 0x1.c55d997da24fdp-22, 0x1.eee32p-1},
- {0x1.d7d82b65c5686p-76, 0x1.c55d997da24fdp-22, 0x1.eee32p-1},
- {-0x1.3f108c0857ca3p-77, -0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
- {-0x1.3f108c0857ca3p-77, -0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
- {-0x1.bd800bca7a221p-78, 0x1.7a4887bd74039p-22, 0x1.fa406p-1},
- {0.0, 0.0, 1.0},
-};
-
-// Look up table for the second range reduction step:
-// Generated by Sollya with:
-// > for i from -64 to 128 do {
-// r = 2^-16 * nearestint(2^16 / (1 + i * 2^-14) );
-// a = -log2(r);
-// b = round(a, D, RN);
-// c = round(a - b, D, RN);
-// print("{", c, ", ", b, "},");
-// };
-LIBC_INLINE_VAR constexpr DoubleDouble LOG2_R2_DD[] = {
- {0x1.ff25180953e64p-62, -0x1.720c2ab2312a9p-8},
- {-0x1.15ffd79560d8fp-62, -0x1.6c4c92b1478ffp-8},
- {0x1.b8d6d6f2e3579p-62, -0x1.668ce3c873549p-8},
- {-0x1.5bfc3f0d5ef71p-62, -0x1.60cd1df6fde91p-8},
- {-0x1.d1f7a8777984ap-64, -0x1.5b0d413c30b5ep-8},
- {0x1.8e858515b8343p-66, -0x1.554d4d97551abp-8},
- {0x1.e165c4014c1f2p-62, -0x1.4f8d4307b46ecp-8},
- {0x1.0f84b2cc14c7ep-63, -0x1.49cd218c9800bp-8},
- {0x1.de618ed0db9a6p-62, -0x1.440ce9254916cp-8},
- {-0x1.f6b8587e64f22p-62, -0x1.3e4c99d110ee7p-8},
- {-0x1.7f793c84cfa63p-64, -0x1.388c338f38bdp-8},
- {-0x1.7d7ecf6258c9ap-65, -0x1.32cbb65f09aeep-8},
- {-0x1.810bc5ac188f5p-62, -0x1.2d0b223fcce81p-8},
- {-0x1.950035fc5b67cp-62, -0x1.274a7730cb841p-8},
- {0x1.4f47f3048cdadp-62, -0x1.2189b5314e95dp-8},
- {0x1.269519861e298p-68, -0x1.1bc8dc409f279p-8},
- {-0x1.5c2b0a46a7e2fp-62, -0x1.1607ec5e063b3p-8},
- {0x1.5001ac8f0bda8p-63, -0x1.1046e588cccap-8},
- {0x1.106f246af5d41p-62, -0x1.0a85c7c03bc4ap-8},
- {0x1.82a00583b34bap-66, -0x1.0354423e3c666p-8},
- {0x1.b6f37deb3137p-65, -0x1.fb25e19f11aecp-9},
- {-0x1.44a2140444811p-63, -0x1.efa310d6550ecp-9},
- {0x1.f5e68a763133fp-63, -0x1.e4201220d4858p-9},
- {0x1.692083115f0b9p-63, -0x1.d89ce57d219a6p-9},
- {0x1.144bb17b9ac9cp-63, -0x1.cd198ae9cdc3dp-9},
- {0x1.ee7f086d32c05p-63, -0x1.c19602656a671p-9},
- {-0x1.d4f1167538dbep-63, -0x1.b6124bee88d82p-9},
- {0x1.7df8d226c67ep-63, -0x1.aa8e6783ba5a2p-9},
- {0x1.60545d61b9512p-63, -0x1.9f0a5523901ebp-9},
- {0x1.54c99c291702p-63, -0x1.938614cc9b468p-9},
- {-0x1.a7e678d7280dep-64, -0x1.8801a67d6ce1p-9},
- {-0x1.6d419bbeb223ap-64, -0x1.7c7d0a3495ec9p-9},
- {0x1.ce2b9892e27e9p-64, -0x1.70f83ff0a7565p-9},
- {-0x1.a4db4eff7bd61p-63, -0x1.657347b031fa2p-9},
- {0x1.5bb04682fab82p-63, -0x1.59ee2171c6a2fp-9},
- {-0x1.78b8bfe6a3adep-64, -0x1.4e68cd33f60a3p-9},
- {0x1.574c3ce9b89b1p-63, -0x1.42e34af550d87p-9},
- {0x1.08fb216647b7bp-63, -0x1.375d9ab467a4dp-9},
- {0x1.ed5a50e7b919cp-66, -0x1.2bd7bc6fcaf56p-9},
- {0x1.91ad7a23f86fep-63, -0x1.2051b0260b3fp-9},
- {0x1.3ab2c932b8b0ap-64, -0x1.14cb75d5b8e54p-9},
- {-0x1.c63bcdf120f7ap-63, -0x1.09450d7d643a9p-9},
- {0x1.8af8c4ab4e82dp-64, -0x1.fb7cee373b008p-10},
- {0x1.a52c2ca9d8b9bp-65, -0x1.e46f655de9cc6p-10},
- {-0x1.460b177a58742p-64, -0x1.cd61806bf5166p-10},
- {0x1.611089de8d12ap-66, -0x1.b6533f5e7cf9bp-10},
- {-0x1.4209853cee70cp-69, -0x1.9f44a232a16eep-10},
- {0x1.964e032541a28p-64, -0x1.8835a8e5824c3p-10},
- {-0x1.fa9f94392637bp-66, -0x1.712653743f454p-10},
- {-0x1.3293693721a53p-64, -0x1.5a16a1dbf7eb6p-10},
- {-0x1.6e2af03c83c6ep-68, -0x1.43069419cbad5p-10},
- {-0x1.b5f05b9d5bd29p-65, -0x1.2bf62a2ad9d74p-10},
- {0x1.3db883c072f72p-64, -0x1.14e5640c4193p-10},
- {-0x1.a675a1c045304p-68, -0x1.fba8837643cf6p-11},
- {0x1.3b9c2aeb00068p-66, -0x1.cd85866933743p-11},
- {-0x1.2911a381901ebp-66, -0x1.9f61d0eb8f98bp-11},
- {-0x1.5ea75a74def03p-68, -0x1.713d62f7957c3p-11},
- {-0x1.305b92f93ffep-67, -0x1.43183c878218dp-11},
- {0x1.b7c8c8dd40d35p-68, -0x1.14f25d959223ap-11},
- {0x1.dc915d58a62f6p-66, -0x1.cd978c3804191p-12},
- {0x1.c7bc3fe53cd94p-66, -0x1.7148ec2a1bfc9p-12},
- {-0x1.427ce595cc53cp-67, -0x1.14f8daf5e3bcfp-12},
- {-0x1.d523885ac824cp-67, -0x1.714eb11fa5363p-13},
- {-0x1.945957f63330ap-69, -0x1.715193b17d35dp-14},
- {0, 0},
- {-0x1.88fb2ea8bf9eap-70, 0x1.7157590356aeep-14},
- {-0x1.5aeaee345d04ep-68, 0x1.715a3bc3593d5p-13},
- {-0x1.7fce430230132p-66, 0x1.1505d6ee104c5p-12},
- {-0x1.9a480f204ff09p-70, 0x1.716001718cb2bp-12},
- {-0x1.00e7233f2d8bdp-68, 0x1.cdbb9d77ae5a8p-12},
- {0x1.09d379fa18c5dp-67, 0x1.150c5586012b8p-11},
- {0x1.b6b9d90a104d3p-65, 0x1.433b951d0b231p-11},
- {0x1.4d9a3ea651885p-65, 0x1.716b8d86bc285p-11},
- {-0x1.7590b3a76f0f9p-67, 0x1.9f9c3ec8db94fp-11},
- {0x1.f183ca5b21bfep-65, 0x1.cdcda8e93107fp-11},
- {-0x1.a7e3465ba127p-66, 0x1.fbffcbed8465fp-11},
- {-0x1.7821f738d1221p-64, 0x1.151953edceec6p-10},
- {0x1.3bb4c0fb95359p-65, 0x1.2c331e5ca2e7dp-10},
- {0x1.236028e962f8p-64, 0x1.434d4546227fcp-10},
- {0x1.aaaa64d30f184p-66, 0x1.5a67c8ad32315p-10},
- {-0x1.a821b7cc57a7ap-64, 0x1.7182a894b69c6p-10},
- {-0x1.13d9d78aace21p-64, 0x1.889de4ff94838p-10},
- {-0x1.2f249a6b923ap-64, 0x1.9fb97df0b0cc2p-10},
- {-0x1.d47dc3664be7ap-68, 0x1.b6d5736af07e6p-10},
- {0x1.bd1522c6418fbp-64, 0x1.cdf1c57138c53p-10},
- {-0x1.bacdbb22d2163p-64, 0x1.e50e74066eee6p-10},
- {-0x1.ca7604812d77bp-64, 0x1.fc2b7f2d786a5p-10},
- {-0x1.2b6832f8830bfp-63, 0x1.09a473749d663p-9},
- {0x1.4e712033d0457p-65, 0x1.1533559e4de55p-9},
- {-0x1.473dd044017b5p-66, 0x1.20c26615409f1p-9},
- {-0x1.e033bcac726d3p-63, 0x1.2c51a4dae8915p-9},
- {-0x1.4a47a2b18a0fap-63, 0x1.37e111f0b8cb5p-9},
- {0x1.6f3615771c17bp-66, 0x1.4370ad58246ddp-9},
- {0x1.c0ee6c32d6236p-65, 0x1.4f0077129eabp-9},
- {0x1.fa94c99761b8fp-64, 0x1.5a906f219ac67p-9},
- {-0x1.979e6b473fbf8p-64, 0x1.662095868c153p-9},
- {0x1.30edde8d24c7bp-64, 0x1.71b0ea42e5fdap-9},
- {-0x1.d01594fe1421cp-64, 0x1.7d416d581bf7cp-9},
- {0x1.50bf7b995b49ap-63, 0x1.88d21ec7a18cdp-9},
- {-0x1.28ea2bcec5018p-63, 0x1.9462fe92ea57cp-9},
- {0x1.ed6add489c30bp-65, 0x1.9ff40cbb6a04bp-9},
- {0x1.201d5c3bbeb69p-64, 0x1.ab85494294517p-9},
- {-0x1.a05d0d4461ea9p-64, 0x1.b716b429dd0d3p-9},
- {-0x1.7c974c8a392fdp-63, 0x1.c2a84d72b8189p-9},
- {-0x1.f068238451bdep-64, 0x1.ce3a151e9965bp-9},
- {-0x1.5e4d95c6259c3p-66, 0x1.d9cc0b2ef4f83p-9},
- {-0x1.1fc262efaad6cp-63, 0x1.e55e2fa53ee53p-9},
- {0x1.49eee7abc7716p-63, 0x1.f0f08282eb533p-9},
- {-0x1.903de284d2782p-65, 0x1.fc8303c96e7a6p-9},
- {-0x1.ec564845134cbp-63, 0x1.040ad9bd1e522p-8},
- {-0x1.7692b7791cf1fp-66, 0x1.0861eadabc3dcp-8},
- {-0x1.37829afb11c1p-62, 0x1.0e2b6b51e4f7ep-8},
- {0x1.6706b91c3b0bap-62, 0x1.13f5030033459p-8},
- {-0x1.7558ccd710756p-62, 0x1.19beb1e6616c9p-8},
- {0x1.79f72a5bbe9dep-62, 0x1.1f88780529bb1p-8},
- {-0x1.e1297c110b25p-62, 0x1.2552555d46886p-8},
- {0x1.29930d567ca26p-62, 0x1.2b1c49ef72343p-8},
- {0x1.a08cbd7592a17p-65, 0x1.30e655bc67275p-8},
- {0x1.e4f9d4ac5db83p-62, 0x1.36b078c4dfd31p-8},
- {-0x1.ed1b0aafd30c2p-62, 0x1.3c7ab30996b1cp-8},
- {0x1.e78f0aa014b32p-62, 0x1.4245048b46462p-8},
- {0x1.8594548038a0fp-69, 0x1.480f6d4aa91c2p-8},
- {0x1.3df498168a333p-63, 0x1.4dd9ed4879c82p-8},
- {0x1.b1c502544f82ap-62, 0x1.53a4848572e77p-8},
- {-0x1.dc50552fe0da9p-63, 0x1.596f33024f203p-8},
- {-0x1.671d85c357d5ep-62, 0x1.5f39f8bfc9212p-8},
- {0x1.1c670cabccefap-64, 0x1.6504d5be9ba1ep-8},
- {-0x1.9983a9e98f318p-62, 0x1.6acfc9ff8162fp-8},
- {0x1.ae1a26af3eebep-62, 0x1.709ad583352d6p-8},
- {0x1.655eb510bfda3p-62, 0x1.7665f84a71d35p-8},
- {-0x1.e287bc0192e15p-64, 0x1.7c313255f22f8p-8},
- {0x1.cc4944139ccbfp-63, 0x1.81fc83a671257p-8},
- {0x1.4e09b4cb8645bp-62, 0x1.87c7ec3ca9a19p-8},
- {-0x1.5becc991e3a5fp-64, 0x1.8d936c1956991p-8},
- {-0x1.ddfa3f1e15ba8p-62, 0x1.935f033d3309ep-8},
- {-0x1.b7b06ea3fb362p-62, 0x1.992ab1a8f9facp-8},
- {0x1.32d614904e46cp-62, 0x1.9ef6775d667b4p-8},
- {-0x1.7186892b5bfaep-64, 0x1.a4c2545b33a3ep-8},
- {-0x1.d4de10b28dfd8p-62, 0x1.aa8e48a31c95cp-8},
- {0x1.4bb4b3bdc8175p-62, 0x1.b05a5435dc7adp-8},
- {0x1.9cedbd1d7fba5p-62, 0x1.b62677142e86p-8},
- {-0x1.0ed3379beaffdp-66, 0x1.bbf2b13ecdf2fp-8},
- {0x1.6e86a125567a6p-62, 0x1.c1bf02b67606p-8},
- {-0x1.35038e0c0a52cp-62, 0x1.c6184f1b326d9p-8},
- {0x1.05ef8bf5adf5ep-67, 0x1.cbe4c95b6c5abp-8},
- {-0x1.b7338b99a6b26p-65, 0x1.d1b15aeab217cp-8},
- {0x1.9e901c30c427ep-63, 0x1.d77e03c9bf0a4p-8},
- {-0x1.1f28a9c0b3d47p-62, 0x1.dd4ac3f94ea0ap-8},
- {-0x1.140ef760d3b63p-62, 0x1.e3179b7a1c52p-8},
- {-0x1.ab65b1037f517p-63, 0x1.e8e48a4ce39e7p-8},
- {-0x1.76940c457ce6dp-63, 0x1.eeb19072600edp-8},
- {0x1.da3ae65a605cfp-64, 0x1.f47eadeb4d34dp-8},
- {0x1.b15d0bce2ede6p-62, 0x1.fa4be2b866abp-8},
- {0x1.e02aa1fa9dc57p-61, 0x1.000c976d340a6p-7},
- {0x1.6be971a5565b9p-62, 0x1.02f34929068f3p-7},
- {-0x1.8a9319a6ed164p-64, 0x1.05da069008be7p-7},
- {0x1.825079f1e0ec5p-62, 0x1.08c0cfa298771p-7},
- {0x1.60d5749321466p-63, 0x1.0ba7a461139c8p-7},
- {-0x1.5b8f4c479e2ep-61, 0x1.0e8e84cbd8169p-7},
- {-0x1.e3e1248004e29p-62, 0x1.117570e343d17p-7},
- {0x1.9ac06487c375p-63, 0x1.145c68a7b4bddp-7},
- {0x1.f657ea5c03ea4p-62, 0x1.17436c1988d0dp-7},
- {-0x1.5a965659a05e2p-61, 0x1.1a2a7b391e04p-7},
- {-0x1.21ce9b9bfc512p-61, 0x1.1d119606d2554p-7},
- {-0x1.30fda247ad0e1p-61, 0x1.1ff8bc8303c7p-7},
- {-0x1.382c78a45cdeap-62, 0x1.22dfeeae10601p-7},
- {0x1.46ae4a64073d4p-61, 0x1.250d5bf952374p-7},
- {-0x1.dcad2cec3b84bp-62, 0x1.27f4a29740a2fp-7},
- {-0x1.413fbeb0b0635p-61, 0x1.2adbf4e50cdf9p-7},
- {0x1.f28e6a48bcb9p-61, 0x1.2dc352e315049p-7},
- {-0x1.96f286e1eb086p-61, 0x1.30aabc91b72ep-7},
- {-0x1.f88c04206dfa1p-61, 0x1.339231f1517c1p-7},
- {-0x1.11ea20e195841p-61, 0x1.3679b30242139p-7},
- {-0x1.d6e71452b674ap-63, 0x1.39613fc4e71dcp-7},
- {-0x1.57c578233b1b3p-61, 0x1.3c48d8399ec85p-7},
- {-0x1.ec430f03b76ep-63, 0x1.3f307c60c7455p-7},
- {0x1.e00dd1902ffb9p-61, 0x1.42182c3abecb5p-7},
- {-0x1.f22bcd96afe38p-61, 0x1.44ffe7c7e3957p-7},
- {0x1.08fd90f841d3p-61, 0x1.47e7af0893e2fp-7},
- {0x1.09594c5552bccp-62, 0x1.4acf81fd2df7ep-7},
- {-0x1.01a8a652e5602p-61, 0x1.4db760a6101c9p-7},
- {-0x1.826168febb3dp-64, 0x1.509f4b03989dcp-7},
- {-0x1.7eb21a35021e3p-62, 0x1.5387411625cccp-7},
- {-0x1.66cbc818e175p-61, 0x1.566f42de15ff4p-7},
- {0x1.9b784dd6cebdap-64, 0x1.5957505bc78f6p-7},
- {0x1.2b121ab482456p-61, 0x1.5b8562298c65bp-7},
- {-0x1.5d29869dd8233p-62, 0x1.5e6d842633702p-7},
- {-0x1.572a1b6cd63cfp-61, 0x1.6155b1d99f672p-7},
- {-0x1.a1f355360e877p-62, 0x1.643deb442eb59p-7},
- {-0x1.b6f1cd2e1c03fp-61, 0x1.672630663fcadp-7},
- {-0x1.2aaa11ccddcaep-61, 0x1.6a0e8140311aap-7},
- {0x1.3d979ddf4746cp-61, 0x1.6cf6ddd2611d4p-7},
- {-0x1.dc930484501f8p-63, 0x1.6fdf461d2e4f8p-7},
-};
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-LIBC_INLINE bool is_odd_integer(float x) {
- using FPBits = typename fputil::FPBits<float>;
- uint32_t x_u = cpp::bit_cast<uint32_t>(x);
- int32_t x_e =
- static_cast<int32_t>((x_u & FPBits::EXP_MASK) >> FPBits::FRACTION_LEN);
- int32_t lsb = cpp::countr_zero(x_u | FPBits::EXP_MASK);
- constexpr int32_t UNIT_EXPONENT =
- FPBits::EXP_BIAS + static_cast<int32_t>(FPBits::FRACTION_LEN);
- return (x_e + lsb == UNIT_EXPONENT);
-}
-
-LIBC_INLINE bool is_integer(float x) {
- using FPBits = typename fputil::FPBits<float>;
- uint32_t x_u = cpp::bit_cast<uint32_t>(x);
- int32_t x_e =
- static_cast<int32_t>((x_u & FPBits::EXP_MASK) >> FPBits::FRACTION_LEN);
- int32_t lsb = cpp::countr_zero(x_u | FPBits::EXP_MASK);
- constexpr int32_t UNIT_EXPONENT =
- FPBits::EXP_BIAS + static_cast<int32_t>(FPBits::FRACTION_LEN);
- return (x_e + lsb >= UNIT_EXPONENT);
-}
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-LIBC_INLINE bool larger_exponent(double a, double b) {
- using DoubleBits = typename fputil::FPBits<double>;
- return DoubleBits(a).get_biased_exponent() >=
- DoubleBits(b).get_biased_exponent();
-}
-
-// Calculate 2^(y * log2(x)) in double-double precision.
-// At this point we can reuse the following values:
-// idx_x: index for extra precision of log2 for the middle part of log2(x).
-// dx: the reduced argument for log2(x)
-// y6: 2^6 * y.
-// lo6_hi: the high part of 2^6 * (y - (hi + mid))
-// exp2_hi_mid: high part of 2^(hi + mid)
-LIBC_INLINE double powf_double_double(int idx_x, double dx, double y6,
- double lo6_hi,
- const DoubleDouble &exp2_hi_mid) {
- using DoubleBits = typename fputil::FPBits<double>;
-
- // Perform a second range reduction step:
- // idx2 = round(2^14 * (dx + 2^-8)) = round ( dx * 2^14 + 2^6)
- // dx2 = (1 + dx) * r2 - 1
- // Output range:
- // -0x1.3ffcp-15 <= dx2 <= 0x1.3e3dp-15
- int idx2 = static_cast<int>(
- fputil::nearest_integer(fputil::multiply_add(dx, 0x1.0p14, 0x1.0p6)));
- double dx2 = fputil::multiply_add(
- 1.0 + dx, common_constants_internal::R2[idx2], -1.0); // Exact
-
- // Degree-5 polynomial approximation of log2(1 + x)/x in double-double
- // precision. Generate by Solya with:
- // > P = fpminimax(log2(1 + x)/x, 5, [|DD...|],
- // [-0x1.3ffcp-15, 0x1.3e3dp-15]);
- // > dirtyinfnorm(log2(1 + x)/x - P, [-0x1.3ffcp-15, 0x1.3e3dp-15]);
- // 0x1.8be5...p-96.
- constexpr DoubleDouble COEFFS[] = {
- {0x1.777d0ffda25ep-56, 0x1.71547652b82fep0},
- {-0x1.777d101cf0a84p-57, -0x1.71547652b82fep-1},
- {0x1.ce04b5140d867p-56, 0x1.ec709dc3a03fdp-2},
- {0x1.137b47e635be5p-56, -0x1.71547652b82fbp-2},
- {-0x1.b5a30b3bdb318p-58, 0x1.2776c516a92a2p-2},
- {0x1.2d2fbd081e657p-57, -0x1.ec70af1929ca6p-3},
- };
-
- DoubleDouble dx_dd({0.0, dx2});
- DoubleDouble p = fputil::polyeval(dx_dd, COEFFS[0], COEFFS[1], COEFFS[2],
- COEFFS[3], COEFFS[4], COEFFS[5]);
- // log2(1 + dx2) ~ dx2 * P(dx2)
- DoubleDouble log2_x_lo = fputil::quick_mult(dx2, p);
- // Lower parts of (e_x - log2(r1)) of the first range reduction constant
- DoubleDouble log2_x_mid({LOG2_R_TD[idx_x].lo, LOG2_R_TD[idx_x].mid});
- // -log2(r2) + lower part of (e_x - log2(r1))
- DoubleDouble log2_x_m = fputil::add(LOG2_R2_DD[idx2], log2_x_mid);
- // log2(1 + dx2) - log2(r2) + lower part of (e_x - log2(r1))
- // Since we don't know which one has larger exponent to apply Fast2Sum
- // algorithm, we need to check them before calling double-double addition.
- DoubleDouble log2_x = larger_exponent(log2_x_m.hi, log2_x_lo.hi)
- ? fputil::add(log2_x_m, log2_x_lo)
- : fputil::add(log2_x_lo, log2_x_m);
- DoubleDouble lo6_hi_dd({0.0, lo6_hi});
- // 2^6 * y * (log2(1 + dx2) - log2(r2) + lower part of (e_x - log2(r1)))
- DoubleDouble prod = fputil::quick_mult(y6, log2_x);
- // 2^6 * (y * log2(x) - (hi + mid)) = 2^6 * lo
- DoubleDouble lo6 = larger_exponent(prod.hi, lo6_hi)
- ? fputil::add(prod, lo6_hi_dd)
- : fputil::add(lo6_hi_dd, prod);
-
- constexpr DoubleDouble EXP2_COEFFS[] = {
- {0, 0x1p0},
- {0x1.abc9e3b398024p-62, 0x1.62e42fefa39efp-7},
- {-0x1.5e43a5429bddbp-69, 0x1.ebfbdff82c58fp-15},
- {-0x1.d33162491268fp-77, 0x1.c6b08d704a0cp-23},
- {0x1.4fb32d240a14ep-86, 0x1.3b2ab6fba4e77p-31},
- {0x1.e84e916be83ep-97, 0x1.5d87fe78a6731p-40},
- {-0x1.9a447bfddc5e6p-103, 0x1.430912f86bfb8p-49},
- {-0x1.31a55719de47fp-113, 0x1.ffcbfc588ded9p-59},
- {-0x1.0ba57164eb36bp-122, 0x1.62c034beb8339p-68},
- {-0x1.8483eabd9642dp-132, 0x1.b5251ff97bee1p-78},
- };
-
- DoubleDouble pp = fputil::polyeval(
- lo6, EXP2_COEFFS[0], EXP2_COEFFS[1], EXP2_COEFFS[2], EXP2_COEFFS[3],
- EXP2_COEFFS[4], EXP2_COEFFS[5], EXP2_COEFFS[6], EXP2_COEFFS[7],
- EXP2_COEFFS[8], EXP2_COEFFS[9]);
- DoubleDouble rr = fputil::quick_mult(exp2_hi_mid, pp);
-
- // Make sure the sum is normalized:
- DoubleDouble r = fputil::exact_add(rr.hi, rr.lo);
- // Round to odd.
- uint64_t r_bits = cpp::bit_cast<uint64_t>(r.hi);
- if (LIBC_UNLIKELY(((r_bits & 0xfff'ffff) == 0) && (r.lo != 0.0))) {
- Sign hi_sign = DoubleBits(r.hi).sign();
- Sign lo_sign = DoubleBits(r.lo).sign();
- if (hi_sign == lo_sign) {
- ++r_bits;
- } else if ((r_bits & DoubleBits::FRACTION_MASK) > 0) {
- --r_bits;
- }
- }
-
- return cpp::bit_cast<double>(r_bits);
-}
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-} // namespace powf_internal
-
-LIBC_INLINE float powf(float x, float y) {
- using namespace powf_internal;
- using FloatBits = typename fputil::FPBits<float>;
- using DoubleBits [[maybe_unused]] = typename fputil::FPBits<double>;
-
- FloatBits xbits(x), ybits(y);
-
- uint32_t x_u = xbits.uintval();
- uint32_t x_abs = xbits.abs().uintval();
- uint32_t y_u = ybits.uintval();
- uint32_t y_abs = ybits.abs().uintval();
-
- ///////// BEGIN - Check exceptional cases ////////////////////////////////////
-
- // The single precision number that is closest to 1 is (1 - 2^-24), which has
- // log2(1 - 2^-24) ~ -1.715...p-24.
- // So if |y| > 151 * 2^24, and x is finite:
- // |y * log2(x)| = 0 or > 151.
- // Hence x^y will either overflow or underflow if x is not zero.
- if (LIBC_UNLIKELY((y_abs & 0x0007'ffff) == 0) || (y_abs > 0x4f170000)) {
- // y is signaling NaN
- if (xbits.is_signaling_nan() || ybits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FloatBits::quiet_nan().get_val();
- }
-
- // Exceptional exponents.
- if (y == 0.0f)
- return 1.0f;
-
- switch (y_abs) {
- case 0x7f80'0000: { // y = +-Inf
- if (x_abs > 0x7f80'0000) {
- // pow(NaN, +-Inf) = NaN
- return x;
- }
- if (x_abs == 0x3f80'0000) {
- // pow(+-1, +-Inf) = 1.0f
- return 1.0f;
- }
- if (x == 0.0f && y_u == 0xff80'0000) {
- // pow(+-0, -Inf) = +inf and raise FE_DIVBYZERO
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_DIVBYZERO);
- return FloatBits::inf().get_val();
- }
- // pow (|x| < 1, -inf) = +inf
- // pow (|x| < 1, +inf) = 0.0f
- // pow (|x| > 1, -inf) = 0.0f
- // pow (|x| > 1, +inf) = +inf
- return ((x_abs < 0x3f80'0000) == (y_u == 0xff80'0000))
- ? FloatBits::inf().get_val()
- : 0.0f;
- }
- default:
- // Speed up for common exponents
- float r = fputil::sqrt<float>(x);
- switch (y_u) {
- case 0x3f00'0000: // y = 0.5f
- // pow(x, 1/2) = sqrt(x)
- if (LIBC_UNLIKELY(x == 0.0f || x_u == 0xff80'0000)) {
- // pow(-0, 1/2) = +0
- // pow(-inf, 1/2) = +inf
- // Make sure it is correct for FTZ/DAZ.
- return x * x;
- }
- return (FloatBits(r).uintval() != 0x8000'0000) ? r : 0.0f;
- case 0x3f80'0000: // y = 1.0f
- return x;
- case 0x4000'0000: // y = 2.0f
- // pow(x, 2) = x^2
- return x * x;
- // TODO: Enable special case speed-up for x^(-1/2) when rsqrt is ready.
- // case 0xbf00'0000: // pow(x, -1/2) = rsqrt(x)
- // return rsqrt(x);
- }
- if (is_integer(y) && (y_u > 0x4000'0000) && (y_u <= 0x41c0'0000)) {
- // Check for exact cases when 2 < y < 25 and y is an integer.
- int msb =
- (x_abs == 0) ? (FloatBits::TOTAL_LEN - 2) : cpp::countl_zero(x_abs);
- msb = (msb > FloatBits::EXP_LEN) ? msb : FloatBits::EXP_LEN;
- int lsb = (x_abs == 0) ? 0 : cpp::countr_zero(x_abs);
- lsb = (lsb > FloatBits::FRACTION_LEN) ? FloatBits::FRACTION_LEN : lsb;
- int extra_bits = FloatBits::TOTAL_LEN - 2 - lsb - msb;
- int iter = static_cast<int>(y);
-
- if (extra_bits * iter <= FloatBits::FRACTION_LEN + 2) {
- // The result is either exact or exactly half-way.
- // But it is exactly representable in double precision.
- double x_d = static_cast<double>(x);
- double result = x_d;
- for (int i = 1; i < iter; ++i)
- result *= x_d;
- return static_cast<float>(result);
- }
- }
- if (y_abs > 0x4f17'0000) {
- // if y is NaN
- if (y_abs > 0x7f80'0000) {
- if (x_u == 0x3f80'0000) { // x = 1.0f
- // pow(1, NaN) = 1
- return 1.0f;
- }
- // pow(x, NaN) = NaN
- return y;
- }
- // x^y will be overflow / underflow in single precision. Set y to a
- // large enough exponent but not too large, so that the computations
- // won't be overflow in double precision.
- y = cpp::bit_cast<float>((y_u & FloatBits::SIGN_MASK) + 0x4f800000U);
- }
- }
- }
-
- int ex = -FloatBits::EXP_BIAS;
- uint64_t sign = 0;
-
- // y is finite and non-zero.
- if (LIBC_UNLIKELY(((x_u & 0x801f'ffffU) == 0) || x_u >= 0x7f80'0000U ||
- x_u < 0x0080'0000U)) {
- // if x is signaling NaN
- if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FloatBits::quiet_nan().get_val();
- }
-
- switch (x_u) {
- case 0x3f80'0000: // x = 1.0f
- return 1.0f;
- // TODO: Put these 2 entrypoint dependency under control flag.
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- case 0x4000'0000: // x = 2.0f
- // pow(2, y) = exp2(y)
- return math::exp2f(y);
- case 0x4120'0000: // x = 10.0f
- // pow(10, y) = exp10(y)
- return math::exp10f(y);
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- }
-
- const bool x_is_neg = x_u >= FloatBits::SIGN_MASK;
-
- if (x == 0.0f) {
- const bool out_is_neg =
- x_is_neg && is_odd_integer(FloatBits(y_u).get_val());
- if (y_u > 0x8000'0000U) {
- // pow(0, negative number) = inf
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_DIVBYZERO);
- return FloatBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
- }
- // pow(0, positive number) = 0
- return out_is_neg ? -0.0f : 0.0f;
- }
-
- if (x_abs == 0x7f80'0000) {
- // x = +-Inf
- const bool out_is_neg =
- x_is_neg && is_odd_integer(FloatBits(y_u).get_val());
- if (y_u >= FloatBits::SIGN_MASK) {
- return out_is_neg ? -0.0f : 0.0f;
- }
- return FloatBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
- }
-
- if (x_abs > 0x7f80'0000) {
- // x is NaN.
- // pow (aNaN, 0) is already taken care above.
- return x;
- }
-
- // Normalize denormal inputs.
- if (x_abs < 0x0080'0000U) {
- ex -= 64;
- x *= 0x1.0p64f;
- }
-
- // x is finite and negative, and y is a finite integer.
- if (x_is_neg) {
- if (is_integer(y)) {
- x = -x;
- if (is_odd_integer(y)) {
- // sign = -1.0;
- sign = 0x8000'0000'0000'0000ULL;
- }
- } else {
- // pow( negative, non-integer ) = NaN
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- return FloatBits::quiet_nan().get_val();
- }
- }
- }
-
- ///////// END - Check exceptional cases //////////////////////////////////////
-
-#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \
- defined(LIBC_MATH_HAS_SMALL_TABLES)
- return powf_small_tables(x, ex, sign, y);
-#else
- using namespace common_constants_internal;
-
- // x^y = 2^( y * log2(x) )
- // = 2^( y * ( e_x + log2(m_x) ) )
- // First we compute log2(x) = e_x + log2(m_x)
- x_u = FloatBits(x).uintval();
-
- // Extract exponent field of x.
- ex += (x_u >> FloatBits::FRACTION_LEN);
- double e_x = static_cast<double>(ex);
- // Use the highest 7 fractional bits of m_x as the index for look up tables.
- uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
- int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 7));
- // Add the hidden bit to the mantissa.
- // 1 <= m_x < 2
- float m_x = cpp::bit_cast<float>(x_mant | 0x3f800000);
-
- // 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).
-#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
- double dx =
- static_cast<double>(fputil::multiply_add(m_x, R[idx_x], -1.0f)); // Exact
-#else
- double dx =
- fputil::multiply_add(static_cast<double>(m_x), RD[idx_x], -1.0); // Exact
-#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
-
- // Degree-5 polynomial approximation:
- // dx * P(dx) ~ log2(1 + dx)
- // Generated by Sollya with:
- // > P = fpminimax(log2(1 + x)/x, 5, [|D...|], [-2^-8, 2^-7]);
- // > dirtyinfnorm(log2(1 + x)/x - P, [-2^-8, 2^-7]);
- // 0x1.653...p-52
- constexpr double COEFFS[] = {0x1.71547652b82fep0, -0x1.71547652b7a07p-1,
- 0x1.ec709dc458db1p-2, -0x1.715479c2266c9p-2,
- 0x1.2776ae1ddf8fp-2, -0x1.e7b2178870157p-3};
-
- double dx2 = dx * dx; // Exact
- 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 p = fputil::polyeval(dx2, c0, c1, c2);
-
- //////////////////////////////////////////////////////////////////////////////
- // NOTE: For some reason, this is significantly less efficient than above!
- //
- // > P = fpminimax(log2(1 + x)/x, 4, [|D...|], [-2^-8, 2^-7]);
- // > dirtyinfnorm(log2(1 + x)/x - P, [-2^-8, 2^-7]);
- // 0x1.d04...p-44
- // constexpr double COEFFS[] = {0x1.71547652b8133p0, -0x1.71547652d1e33p-1,
- // 0x1.ec70a098473dep-2, -0x1.7154c5ccdf121p-2,
- // 0x1.2514fd90a130ap-2};
- //
- // 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 p = fputil::polyeval(dx2, c0, c1, COEFFS[4]);
- //////////////////////////////////////////////////////////////////////////////
-
- // s = e_x - log2(r) + dx * P(dx)
- // Approximation errors:
- // |log2(x) - s| < ulp(e_x) + (bounds on dx) * (error bounds of P(dx))
- // = ulp(e_x) + 2^-7 * 2^-51
- // < 2^8 * 2^-52 + 2^-7 * 2^-43
- // ~ 2^-44 + 2^-50
- double s = fputil::multiply_add(dx, p, LOG2_R[idx_x] + e_x);
-
- // To compute 2^(y * log2(x)), we break the exponent into 3 parts:
- // y * log(2) = 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| < 150, (hi + mid) uses (8 + 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
-
- // 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 = static_cast<double>(y * 0x1.0p6f); // Exact.
- double hm = fputil::nearest_integer(s * y6);
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // lo6 = 2^6 * lo.
- double lo6_hi =
- fputil::multiply_add(y6, e_x + LOG2_R_DD[idx_x].hi, -hm); // Exact
- // Error bounds:
- // | (y*log2(x) - hm * 2^-6 - lo) / y| < err(dx * p) + err(LOG2_R_DD.lo)
- // < 2^-51 + 2^-75
- double lo6 = fputil::multiply_add(
- y6, fputil::multiply_add(dx, p, LOG2_R_DD[idx_x].lo), lo6_hi);
-#else
- // lo6 = 2^6 * lo.
- double lo6_hi =
- fputil::multiply_add(y6, e_x + LOG2_R_TD[idx_x].hi, -hm); // Exact
- // Error bounds:
- // | (y*log2(x) - hm * 2^-6 - lo) / y| < err(dx * p) + err(LOG2_R_DD.lo)
- // < 2^-51 + 2^-75
- double lo6 = fputil::multiply_add(
- y6, fputil::multiply_add(dx, p, LOG2_R_TD[idx_x].mid), lo6_hi);
-#endif
-
- // |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.
- int64_t hm_i =
- cpp::clamp<int64_t>(static_cast<int64_t>(hm), -(1 << 15), 1 << 15);
-
- int idx_y = hm_i & 0x3f;
-
- // 2^hi
- int64_t exp_hi_i = (hm_i >> 6) << DoubleBits::FRACTION_LEN;
- // 2^mid
- int64_t exp_mid_i = cpp::bit_cast<uint64_t>(EXP2_MID1[idx_y].hi);
- // (-1)^sign * 2^hi * 2^mid
- // Error <= 2^hi * 2^-53
- uint64_t exp2_hi_mid_i = static_cast<uint64_t>(exp_hi_i + exp_mid_i) + sign;
- double exp2_hi_mid = cpp::bit_cast<double>(exp2_hi_mid_i);
-
- // 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[1], EXP2_COEFFS[0]);
- double d1 = fputil::multiply_add(lo6, EXP2_COEFFS[3], EXP2_COEFFS[2]);
- double d2 = fputil::multiply_add(lo6, EXP2_COEFFS[5], EXP2_COEFFS[4]);
- double pp = fputil::polyeval(lo6_sqr, d0, d1, d2);
-
- double r = pp * exp2_hi_mid;
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- return static_cast<float>(r);
-#else
- // Ziv accuracy test.
- uint64_t r_u = cpp::bit_cast<uint64_t>(r);
- float r_upper = static_cast<float>(cpp::bit_cast<double>(r_u + ERR));
- float r_lower = static_cast<float>(cpp::bit_cast<double>(r_u - ERR));
-
- if (LIBC_LIKELY(r_upper == r_lower)) {
- // Check for overflow or underflow.
- if (LIBC_UNLIKELY(FloatBits(r_upper).get_mantissa() == 0)) {
- if (FloatBits(r_upper).is_inf()) {
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_OVERFLOW);
- } else if (r_upper == 0.0f) {
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_UNDERFLOW);
- }
- }
- return r_upper;
- }
-
- // Scale lower part of 2^(hi + mid)
- DoubleDouble exp2_hi_mid_dd;
- exp2_hi_mid_dd.lo =
- (idx_y != 0)
- ? cpp::bit_cast<double>(exp_hi_i +
- cpp::bit_cast<int64_t>(EXP2_MID1[idx_y].mid))
- : 0.0;
- exp2_hi_mid_dd.hi = exp2_hi_mid;
-
- double r_dd = powf_double_double(idx_x, dx, y6, lo6_hi, exp2_hi_mid_dd);
-
- return static_cast<float>(r_dd);
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS && LIBC_MATH_HAS_SMALL_TABLES
-}
+using LIBC_MATH_POWF_IMPL::powf;
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
+#undef LIBC_MATH_POWF_IMPL
+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_H
diff --git a/libc/src/__support/math/powf_double_eval.h b/libc/src/__support/math/powf_double_eval.h
new file mode 100644
index 00000000000000..0ed2e8b194df7b
--- /dev/null
+++ b/libc/src/__support/math/powf_double_eval.h
@@ -0,0 +1,548 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Double-precision evaluation implementation for powf(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.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/FPUtil/rounding_mode.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"
+#include "src/__support/math/powf_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+namespace pow_internal = LIBC_NAMESPACE::math::pow_internal;
+namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Accurate evaluation in DoubleDouble precision when Ziv's test fails.
+LIBC_INLINE double powf_accurate(float x, float y, double e_x, uint64_t sign) {
+ using FloatBits = fputil::FPBits<float>;
+ using DoubleBits = fputil::FPBits<double>;
+ using fputil::DoubleDouble;
+
+ uint32_t x_u = FloatBits(x).uintval();
+ uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
+ int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 7));
+ // Append hidden bit.
+ float m_x = cpp::bit_cast<float>(x_mant | 0x3f80'0000);
+ double dx = fputil::multiply_add(
+ static_cast<double>(m_x),
+ static_cast<double>(common_constants_internal::R[idx_x]), -1.0);
+
+ // Degree-11 minimax polynomial for log2(1 + x)/x on [-2^-8, 2^-7] in
+ // double-double, generated by Sollya with:
+ // > display = hexadecimal;
+ // > prec = 256;
+ // > P = fpminimax(log2(1 + x)/x, 11, [|DD...|], [-2^(-8), 2^(-7)]);
+ // > for i from 0 to 11 do {
+ // c = coeff(P, i);
+ // hi = round(c, D, RN);
+ // lo = round(c - hi, D, RN);
+ // print("{", lo, ", ", hi, "},");
+ // };
+ // > dirtyinfnorm((log2(1 + x) - x * P) / x, [-2^(-8), 2^(-7)]);
+ // 0x1.c598...p-104 < 2^-103.
+ constexpr DoubleDouble LOG2_COEFFS[] = {
+ {0x1.777d0ffda0d34p-56, 0x1.71547652b82fep0},
+ {-0x1.777d0ffd88ca8p-57, -0x1.71547652b82fep-1},
+ {0x1.d27f052d9d41dp-56, 0x1.ec709dc3a03fdp-2},
+ {-0x1.777f259c6c0ap-58, -0x1.71547652b82fep-2},
+ {0x1.e5c055599afdbp-56, 0x1.2776c50ef9bfep-2},
+ {0x1.e923056eef156p-58, -0x1.ec709dc3a03fdp-3},
+ {0x1.1f24e26b5ad0dp-57, 0x1.a61762a7add76p-3},
+ {0x1.16d9b5013b55ep-58, -0x1.71547652bf304p-3},
+ {-0x1.1551148c3d1d8p-57, 0x1.484b13f00425ep-3},
+ {-0x1.a1cf327c6703p-60, -0x1.2776ca937f335p-3},
+ {-0x1.d24f6254c2e6ap-57, 0x1.0c9211efe6d7ep-3},
+ {0x1.a81cff2946d7cp-59, -0x1.e1f45baaf7f1dp-4},
+ };
+
+ // Evaluate P(dx) ~ log2(1 + dx) / dx using Horner's scheme.
+ DoubleDouble p = LOG2_COEFFS[11];
+ for (int i = 10; i >= 0; --i)
+ p = fputil::add(LOG2_COEFFS[i], fputil::quick_mult(dx, p));
+
+ // log2(1 + dx) = dx * P(dx) in DoubleDouble.
+ DoubleDouble log2_1p = fputil::quick_mult(dx, p);
+
+ // -log2(r) is represented in TripleDouble as {lo, mid, hi}.
+ // We decompose log2(x) = e_x - log2(r) + log2(1 + dx) into:
+ // log2(x) = (e_x + hi) + (log2_1p + mid + lo)
+ // = (e_x + hi) + log2_tail.
+ // The lower part log2_tail accumulates log2_1p, mid, and lo in DoubleDouble:
+ DoubleDouble log2_tail =
+ fputil::add<false>(log2_1p, powf_internal::LOG2_R_TD[idx_x].mid);
+ log2_tail = fputil::add<false>(log2_tail, powf_internal::LOG2_R_TD[idx_x].lo);
+
+ // Range reduction for 2^(y * log2(x)):
+ // Scale by 2^6 = 64 to use the 64-entry EXP2_MID1 table:
+ // 64 * y * log2(x) = y6 * (e_x + hi) + y6 * log2_tail.
+ //
+ // Since e_x + hi has at most 8 + 20 = 28 bits of precision and y has 24 bits,
+ // the product y6 * (e_x + hi) has at most 28 + 24 = 52 bits, making y6_hi
+ // exact in double precision.
+ constexpr double SCALE = 0x1.0p6;
+ double y6 = static_cast<double>(y) * SCALE;
+ DoubleDouble prod = fputil::quick_mult(y6, log2_tail);
+
+ double e_x_hi = e_x + powf_internal::LOG2_R_TD[idx_x].hi;
+ DoubleDouble y6_hi = fputil::exact_mult(y6, e_x_hi);
+
+ // hm = round(64 * y * log2(x)) is the nearest integer index.
+ double hm_d = y6_hi.hi + prod.hi;
+ double hm = fputil::nearest_integer(hm_d);
+
+ // lo6 = 64 * y * log2(x) - hm = (y6_hi - hm) + prod:
+ DoubleDouble lo6_hi = fputil::exact_add<false>(y6_hi.hi, -hm);
+ lo6_hi.lo += y6_hi.lo;
+
+ DoubleDouble lo6 = fputil::add<false>(prod, lo6_hi);
+
+ // lo = lo6 / 64, with |lo| <= 2^-7.
+ DoubleDouble lo = fputil::quick_mult(0x1.0p-6, lo6);
+
+ int64_t hm_i = static_cast<int64_t>(hm);
+ int idx_y = static_cast<int>(hm_i & 0x3f);
+ int hi = static_cast<int>(hm_i >> 6);
+ int64_t exp_hi_i = static_cast<int64_t>(static_cast<uint64_t>(hi)
+ << DoubleBits::FRACTION_LEN);
+
+ // Degree-8 minimax polynomial for (2^x - 1)/x on [-2^-7, 2^-7] in
+ // double-double, generated by Sollya with:
+ // > display = hexadecimal;
+ // > prec = 256;
+ // > P = fpminimax((2^x - 1)/x, 8, [|DD...|], [-2^(-7), 2^(-7)]);
+ // > for i from 0 to 8 do {
+ // c = coeff(P, i);
+ // hi = round(c, D, RN);
+ // lo = round(c - hi, D, RN);
+ // print("{", lo, ", ", hi, "},");
+ // };
+ // > dirtyinfnorm((2^x - (1 + x * P)) / (2^x), [-2^(-7), 2^(-7)]);
+ // 0x1.e62d...p-106 < 2^-105.
+ constexpr DoubleDouble EXP2_COEFFS_DD[] = {
+ {0x1.abc9e3b39803ep-56, 0x1.62e42fefa39efp-1},
+ {-0x1.5e43a54066432p-57, 0x1.ebfbdff82c58fp-3},
+ {-0x1.d331626e9fe04p-59, 0x1.c6b08d704a0cp-5},
+ {0x1.4f492022eeb74p-62, 0x1.3b2ab6fba4e77p-7},
+ {0x1.035733b00abdp-66, 0x1.5d87fe78a6731p-10},
+ {0x1.9e4f81ec951dep-68, 0x1.430912f86c123p-13},
+ {0x1.7be627c299948p-74, 0x1.ffcbfc588b997p-17},
+ {-0x1.c8a006d5b39acp-74, 0x1.62c03345a64f6p-20},
+ {-0x1.fd14495e4b96p-81, 0x1.b5253eef8c0a9p-24},
+ };
+
+ DoubleDouble exp2_p = EXP2_COEFFS_DD[8];
+ for (int i = 7; i >= 0; --i)
+ exp2_p = fputil::add(EXP2_COEFFS_DD[i], fputil::quick_mult(lo, exp2_p));
+
+ DoubleDouble p_lo = fputil::quick_mult(lo, exp2_p);
+ DoubleDouble exp2_lo = fputil::exact_add(1.0, p_lo.hi);
+ exp2_lo.lo += p_lo.lo;
+
+ DoubleDouble exp2_hi_mid_dd;
+ exp2_hi_mid_dd.hi = cpp::bit_cast<double>(
+ exp_hi_i + cpp::bit_cast<int64_t>(EXP2_MID1[idx_y].hi) + sign);
+ exp2_hi_mid_dd.lo =
+ (idx_y != 0)
+ ? cpp::bit_cast<double>(
+ exp_hi_i + cpp::bit_cast<int64_t>(EXP2_MID1[idx_y].mid) + sign)
+ : 0.0;
+
+ DoubleDouble rr = fputil::quick_mult(exp2_hi_mid_dd, exp2_lo);
+
+ DoubleDouble r = fputil::exact_add(rr.hi, rr.lo);
+
+ // Round to odd to avoid double rounding when converting to float.
+ uint64_t r_bits = cpp::bit_cast<uint64_t>(r.hi);
+ if (LIBC_UNLIKELY(((r_bits & 0x0fff'ffffULL) == 0) && (r.lo != 0.0))) {
+ if (DoubleBits(r.hi).sign() == DoubleBits(r.lo).sign()) {
+ ++r_bits;
+ } else if ((r_bits & DoubleBits::FRACTION_MASK) > 0) {
+ --r_bits;
+ }
+ }
+
+ return cpp::bit_cast<double>(r_bits);
+}
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Check if x^y is an exact rounding boundary.
+// When x^y = exact_m * 2^exact_exp with exact_m <= 2^25:
+// - exact_m has at most 26 bits, so it fits in double precision without
+// rounding.
+// - For single precision, exact_exp is in [-200, 127], well within the normal
+// range [-1022, 1023] of double precision.
+// We compute exact_d = exact_m * 2^exact_exp in double precision and cast to
+// float for the final rounded result.
+LIBC_INLINE cpp::optional<float>
+check_exact_boundary(float x, float y, uint64_t sign,
+ Sign out_sign = Sign::POS) {
+ uint32_t exact_m = 0;
+ int exact_exp = 0;
+ if (LIBC_UNLIKELY(powf_internal::is_exact_rounding_boundary(x, y, exact_m,
+ exact_exp))) {
+ // Number of bits in exact_m: 2^(l - 1) <= exact_m < 2^l.
+ int l = 32 - cpp::countl_zero(exact_m);
+ // Unbiased exponent of exact_m * 2^exact_exp.
+ int unbiased_exp = exact_exp + l - 1;
+
+ if (LIBC_UNLIKELY(unbiased_exp > 127))
+ return powf_internal::set_overflow(out_sign);
+
+ // Below the minimum subnormal float (2^-149).
+ if (LIBC_UNLIKELY(unbiased_exp < -149)) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<float>();
+ }
+
+ // Scale factor = sign * 2^exact_exp.
+ using DoubleBits = fputil::FPBits<double>;
+ double scale = cpp::bit_cast<double>(
+ (static_cast<uint64_t>(exact_exp + DoubleBits::EXP_BIAS)
+ << DoubleBits::FRACTION_LEN) |
+ sign);
+ double exact_d = static_cast<double>(exact_m) * scale;
+ return static_cast<float>(exact_d);
+ }
+ return cpp::nullopt;
+}
+
+// Overview of powf(x, y) = x^y computation in double precision:
+//
+// Let x = 2^(e_x) * m_x > 0. Then:
+// x^y = 2^( y * log2(x) )
+// = 2^( y * ( e_x + log2(m_x) ) )
+// = 2^( k + f )
+// = 2^k * 2^f,
+// where:
+// k = round(y * log2(x)),
+// f = y * log2(x) - k.
+//
+// In particular, k is an integer, and |f| <= 0.5.
+// For the final result to fit in single precision:
+// -150 <= y * log2(x) <= 128,
+// and anything outside that range overflows or underflows.
+//
+// Fast pass:
+// 1. Range reduction for log2(m_x) using a 32-entry lookup table:
+// dx = r * m_x - 1 in [-2^-6, 2^-5],
+// log2(m_x) = log2(1 + dx) - log2(r),
+// where -log2(r) is retrieved from LOG2_R_32[idx_x], and log2(1 + dx)/dx is
+// approximated by a degree-6 polynomial in double precision with Estrin's
+// scheme.
+// 2. Exponent reduction:
+// z = y * log2(x),
+// k = round(z),
+// f = z - k in [-0.5, 0.5].
+// 3. 2^f is evaluated using a degree-7 polynomial.
+// Absorbing +-ERR into the constant term 1.0 +- ERR of the final FMA allows
+// computing the upper and lower bounds for Ziv's rounding test in parallel.
+// 4. If upper == lower in single precision, return upper.
+// 5. If Ziv's test fails, fall back to check_exact_boundary and powf_accurate
+// (DoubleDouble / 128-entry table).
+//
+// Accurate pass:
+// 1. Range reduction using a 128-entry table: -log2(r) = {lo, mid, hi}.
+// log2(1 + dx) is evaluated using a degree-11 polynomial in DoubleDouble.
+// 2. Decompose: log2(x) = (e_x + hi) + log2_tail, where:
+// log2_tail = log2(1 + dx) + mid + lo.
+// Because e_x + hi has <= 28 bits of precision and y has 24 bits,
+// y * (e_x + hi) is exact in double precision.
+// 3. Exponent reduction with 64-entry EXP2_MID1 table and degree-8 polynomial
+// in DoubleDouble.
+LIBC_INLINE float powf(float x, float y) {
+ using FloatBits = fputil::FPBits<float>;
+ using DoubleBits = fputil::FPBits<double>;
+ using fputil::DoubleDouble;
+ using namespace common_constants_internal;
+
+ FloatBits xbits(x), ybits(y);
+ uint32_t x_u = xbits.uintval();
+ uint32_t y_u = ybits.uintval();
+ uint32_t y_a = ybits.abs().uintval();
+
+ // Quick filter for special inputs:
+ // - x: 0, +-1, 2^k, +-Inf, NaN.
+ // - y: 0, +-1, +-2, +-0.5, +-Inf, NaN.
+ if (LIBC_UNLIKELY((x_u & 0x001F'FFFF) == 0 || (y_u & 0x007F'FFFF) == 0)) {
+ if (auto r = powf_internal::check_special_inputs(x, y);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ }
+
+ float orig_x_abs = xbits.abs().get_val();
+ float orig_y = y;
+
+ int ex = -FloatBits::EXP_BIAS;
+ uint64_t sign = 0;
+
+ // Check for exceptional cases:
+ // - |y| <= 2^-40 or |y| >= 2^40.
+ // - x < 0, x is subnormal, zero, Inf, or NaN.
+ if (LIBC_UNLIKELY(y_a <= powf_internal::Y_LOWER_BOUND ||
+ y_a >= powf_internal::Y_UPPER_BOUND ||
+ x_u >= FloatBits::inf().uintval() ||
+ x_u < FloatBits::min_normal().uintval())) {
+ if (auto r = powf_internal::check_exceptional_cases(x, y, ex, sign);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ }
+
+ Sign out_sign = (sign == 0) ? Sign::POS : Sign::NEG;
+
+ // x^y = 2^( y * log2(x) ) = 2^( y * ( e_x + log2(m_x) ) )
+ // Compute log2(x) = e_x + log2(m_x)
+ x_u = FloatBits(x).uintval();
+
+ // Extract exponent field of x.
+ ex += (x_u >> FloatBits::FRACTION_LEN);
+ double e_x = static_cast<double>(ex);
+ uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
+ // Top 5 bits of mantissa for 32-entry table lookup.
+ int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 5));
+ // Embed mantissa into double format: m_x in [1.0, 2.0).
+ uint64_t m_x_u =
+ (static_cast<uint64_t>(x_mant) << 29) | 0x3ff0'0000'0000'0000ULL;
+ double m_x = cpp::bit_cast<double>(m_x_u);
+
+ // Range reduction: dx = m_x * r - 1.0 in [-2^-6, 2^-5].
+ // Since r has <= 7 bits and m_x has 24 bits, m_x * r - 1.0 is exact in
+ // double.
+ double dx =
+ fputil::multiply_add(m_x, powf_internal::R_32_D[idx_x], -1.0); // Exact
+
+ // Degree-6 minimax polynomial approximation for log2(1 + dx)/dx on
+ // [-2^-6, 2^-5]. Generated by Sollya with:
+ // > display = hexadecimal;
+ // > p6 = fpminimax(log2(1 + x)/x, 6, [|D...|], [-2^-6, 2^-5]);
+ // > for i from 0 to 6 do print(round(coeff(p6, i), D, RN), ",");
+ // > dirtyinfnorm((log2(1 + x) - x * p6) / log2(1 + x), [-2^-6, 2^-5]);
+ // 0x1.050a...p-47
+ // > dirtyinfnorm(log2(1 + x) - x * p6, [-2^-6, 2^-5]);
+ // 0x1.712b...p-52
+ constexpr double COEFFS[] = {
+ 0x1.71547652b831fp0, -0x1.71547652b2ff5p-1, 0x1.ec709dbd02463p-2,
+ -0x1.7154787f73adep-2, 0x1.2777b7f57b29fp-2, -0x1.ec55ea50eaba7p-3,
+ 0x1.92c78e7d5a931p-3,
+ };
+
+ // Evaluate P(dx) ~ log2(1 + dx)/dx with Estrin's scheme:
+ double dx2 = dx * dx;
+ double lp0 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
+ double lp1 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
+ double lp2 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);
+
+ double dx4 = dx2 * dx2;
+ double lq0 = fputil::multiply_add(dx2, lp1, lp0);
+ double lq1 = fputil::multiply_add(dx2, COEFFS[6], lp2);
+
+ double p = fputil::multiply_add(dx4, lq1, lq0);
+
+ // log2(x) = e_x - log2(r) + dx * P(dx).
+ double log2_r_ex = powf_internal::LOG2_R_32[idx_x] + e_x;
+ double s = fputil::multiply_add(dx, p, log2_r_ex);
+
+ double y_d = static_cast<double>(y);
+ double z = y_d * s;
+
+ // y * log2(x) = k + f, with k integer and |f| <= 0.5.
+ double kd = fputil::nearest_integer(z);
+ int64_t k = static_cast<int64_t>(kd);
+ double f = fputil::multiply_add(y_d, s, -kd);
+
+ // Degree-7 polynomial approximation P(f) ~ (2^f - 1)/f on [-0.5, 0.5]
+ // Generated by Sollya with:
+ // > display = hexadecimal;
+ // > p = fpminimax((2^x - 1)/x, 7, [|D...|], [-0.5, 0.5]);
+ // > for i from 0 to 7 do print(round(coeff(p, i), D, RN), ",");
+ // > dirtyinfnorm(2^x - (1 + x * p), [-0.5, 0.5]);
+ // 0x1.0496...p-39
+ // > dirtyinfnorm((2^x - (1 + x * p)) / 2^x, [-0.5, 0.5]);
+ // 0x1.0496...p-39
+ constexpr double EXP2_COEFFS[] = {
+ 0x1.62e42fef9cdf7p-1, 0x1.ebfbdff85f49p-3, 0x1.c6b08da69c963p-5,
+ 0x1.3b2ab6a0e1172p-7, 0x1.5d87762dfc134p-10, 0x1.43099b4e829a9p-13,
+ 0x1.00c080f7699bep-16, 0x1.62c0108a065d8p-20,
+ };
+
+ // Evaluate P(f) ~ (2^f - 1)/f on [-0.5, 0.5] with Estrin's scheme:
+ double f2 = f * f;
+ double f4 = f2 * f2;
+
+ double p0 = fputil::multiply_add(f, EXP2_COEFFS[1], EXP2_COEFFS[0]);
+ double p1 = fputil::multiply_add(f, EXP2_COEFFS[3], EXP2_COEFFS[2]);
+ double p2 = fputil::multiply_add(f, EXP2_COEFFS[5], EXP2_COEFFS[4]);
+ double p3 = fputil::multiply_add(f, EXP2_COEFFS[7], EXP2_COEFFS[6]);
+
+ double q0 = fputil::multiply_add(f2, p1, p0);
+ double q1 = fputil::multiply_add(f2, p3, p2);
+
+ double poly = fputil::multiply_add(f4, q1, q0);
+
+ // Normal range: -125 <= k <= 128.
+ // (k + 125) as unsigned <= 253 covers all normal outputs.
+ uint64_t k_u = static_cast<uint64_t>(k + 125);
+
+ if (LIBC_LIKELY(k_u <= 253)) {
+ // Scale by 2^k.
+ uint64_t exp2_k_i =
+ (static_cast<uint64_t>(static_cast<int>(k) + DoubleBits::EXP_BIAS)
+ << DoubleBits::FRACTION_LEN) |
+ sign;
+ double exp2_k = cpp::bit_cast<double>(exp2_k_i);
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ // Absorbing +-ERR into 1.0 computes both bounds using two parallel FMAs.
+ constexpr double ERR = 0x1.5p-39;
+ double pp_hi = fputil::multiply_add(f, poly, 1.0 + ERR);
+ double pp_lo = fputil::multiply_add(f, poly, 1.0 - ERR);
+ double r_d = pp_hi * exp2_k;
+#else // !LIBC_TARGET_CPU_HAS_FMA
+ // Without FMA, compute unbiased pp = f * poly + 1.0 directly.
+ double pp = fputil::multiply_add(f, poly, 1.0);
+ double r_d = pp * exp2_k;
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+ float res = static_cast<float>(r_d);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return res;
+#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // Since 2^k is an exact power of 2, testing if the bounds round to the same
+ // float is equivalent to testing the fully scaled bounds.
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ float upper = static_cast<float>(pp_hi);
+ float lower = static_cast<float>(pp_lo);
+#else // !LIBC_TARGET_CPU_HAS_FMA
+ constexpr double ERR = 0x1.5p-39;
+ float upper = static_cast<float>(pp + ERR);
+ float lower = static_cast<float>(pp - ERR);
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+ if (LIBC_LIKELY(upper == lower))
+ return res;
+
+ // Accurate fallback when Ziv's test fails:
+ if (auto r = check_exact_boundary(orig_x_abs, orig_y, sign, out_sign);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
+ double r_dd = powf_accurate(x, y, e_x, sign);
+ res = static_cast<float>(r_dd);
+ if (LIBC_UNLIKELY(FloatBits(res).is_inf()))
+ return powf_internal::set_overflow(out_sign);
+ return res;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // Exceptional path: k > 128 (overflow), k < -155 (underflow), or denormal.
+ if (k > 128)
+ return powf_internal::set_overflow(out_sign);
+
+ if (k < -155)
+ return powf_internal::set_underflow(out_sign);
+
+ // Denormal and underflow path for -155 <= k <= -126:
+ if (auto r = check_exact_boundary(orig_x_abs, orig_y, sign, out_sign);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
+ int64_t exp2_k_i = (static_cast<uint64_t>(k + DoubleBits::EXP_BIAS)
+ << DoubleBits::FRACTION_LEN) |
+ sign;
+ double exp2_k = cpp::bit_cast<double>(exp2_k_i);
+
+ // Add +-2^(-126 + (53 - 24)) = +-2^-97 to mimic single precision denormal
+ // rounding in double precision.
+ // For -155 <= k <= -126, the least significant bit of the double precision
+ // sum (r_d + denorm_bias) is aligned at 2^(-97 - 52) = 2^-149, matching the
+ // least significant bit of single precision denormals 2^(-126 - 23) = 2^-149.
+ // This keeps the values in the double precision normal range while performing
+ // rounding according to the current rounding mode.
+ double denorm_bias = cpp::bit_cast<double>(0x39e0'0000'0000'0000ULL | sign);
+
+ constexpr double ERR = 0x1.5p-39;
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ // Absorbing +-ERR into 1.0 computes both bounds using two parallel FMAs.
+ double pp_hi = fputil::multiply_add(f, poly, 1.0 + ERR);
+ double pp_lo = fputil::multiply_add(f, poly, 1.0 - ERR);
+ double u_hi = fputil::multiply_add(pp_hi, exp2_k, denorm_bias);
+ double u_lo = fputil::multiply_add(pp_lo, exp2_k, denorm_bias);
+#else // !LIBC_TARGET_CPU_HAS_FMA
+ // Without FMA, compute unbiased pp = f * poly + 1.0 directly.
+ double pp = fputil::multiply_add(f, poly, 1.0);
+ double r_d = pp * exp2_k;
+ double err = ERR * r_d;
+ double u_hi = (r_d + err) + denorm_bias;
+ double u_lo = (r_d - err) + denorm_bias;
+#endif // LIBC_TARGET_CPU_HAS_FMA
+
+ if (LIBC_LIKELY(u_hi == u_lo)) {
+ if (LIBC_UNLIKELY(u_hi == denorm_bias))
+ return powf_internal::set_underflow(out_sign);
+
+ float res = static_cast<float>(u_hi - denorm_bias);
+ if (LIBC_UNLIKELY(FloatBits(res).is_normal())) {
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+ double pp = fputil::multiply_add(f, poly, 1.0);
+#endif // LIBC_TARGET_CPU_HAS_FMA
+ if (static_cast<float>(pp) < 1.0f)
+ fputil::raise_underflow_except_if_required<float>();
+ return res;
+ }
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<float>();
+ return res;
+ }
+
+ // Ziv's test failed for denormal input, fall back to accurate pass.
+ double r_dd = powf_accurate(x, y, e_x, sign);
+ float res = static_cast<float>(r_dd);
+ if (LIBC_UNLIKELY(FloatBits(res).is_normal())) {
+ if (static_cast<float>(r_dd * 0x1.0p126) < 1.0f)
+ fputil::raise_underflow_except_if_required<float>();
+ return res;
+ }
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<float>();
+ return res;
+#else // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ if (k > 128)
+ return powf_internal::set_overflow(out_sign);
+ return powf_internal::set_underflow(out_sign);
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
diff --git a/libc/src/__support/math/powf_float_eval.h b/libc/src/__support/math/powf_float_eval.h
new file mode 100644
index 00000000000000..2eda051c1420ea
--- /dev/null
+++ b/libc/src/__support/math/powf_float_eval.h
@@ -0,0 +1,263 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Float-only evaluation implementation for powf(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_FLOAT_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_FLOAT_EVAL_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.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/FPUtil/rounding_mode.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/exp2f_float_utils.h"
+#include "src/__support/math/powf_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace float_eval {
+
+namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;
+
+// Overview of float-precision pow(x, y) = x^y computation:
+//
+// Let x = 2^(e_x) * m_x > 0. Then:
+// x^y = 2^( y * log2(x) )
+// = 2^( y * ( e_x + log2(m_x) ) )
+// = 2^( k + rem )
+// = 2^k * 2^rem,
+// where:
+// k = round(y * log2(x)),
+// rem = y * log2(x) - k.
+//
+// In particular, k is an integer, and |rem| <= 0.5.
+//
+// For the final result to fit in single precision, the exponent field is
+// bounded by:
+// -150 <= y * log2(x) <= 128,
+// and anything outside that range overflows or underflows.
+//
+// Since |rem| <= 0.5:
+// 2^-0.5 <= 2^rem <= 2^0.5,
+// the relative error of x^y = 2^k * 2^rem is approximately:
+// ~ absolute error of rem * ln(2)
+// ~ absolute error of y * log2(x) * ln(2) + rel_err(exp2f_eval)
+// ~ relative error(log2(x)) * |y * log2(x)| * ln(2) + 2^-25.
+//
+// With exp2f_eval providing error ~ 2^-25, the budget for
+// relative error(log2(x)) * |y * log2(x)| * ln(2) is <= 2^-25.
+// For single-precision results, the 8-bit exponent field bounds:
+// |y * log2(x)| < 2^8 = 256.
+// Thus:
+// relative error(log2(x)) <= 2^-25 / (2^8 * ln(2)) < 2^-33.
+//
+// To compute log2(x), we perform range reduction for log2(m_x) using a compact
+// 5-bit lookup table (32 entries):
+// dx = r * m_x - 1 in [-2^-6, 2^-5].
+// Then m_x = (1 + dx) / r, and:
+// log2(m_x) = log2(1 + dx) - log2(r),
+// where -log2(r) is obtained from the look up table LOG2_R_FF_32 with
+// FloatFloat precision (~48 bits). The computation of dx = r * m_x - 1 is
+// exact.
+//
+// The expansion of log2(1 + dx) is:
+// log2(1 + dx) = dx * log2(e) + dx^2 * Q(dx).
+//
+// We compute only the leading term dx * log2(e) in FloatFloat using
+// exact_mult, and evaluate the tail Q(dx) using a degree-4 minimax polynomial
+// purely in single precision:
+// Q(dx) ~ (log2(1 + dx) - dx * log2(e)) / dx^2.
+// Generated by Sollya with:
+// > log2_e = 0x1.715476p0 + 0x1.4ae0cp-26;
+// > f = (log2(1 + x) - x * log2_e) / x^2;
+// > Q = fpminimax(f, 4, [|SG...|], [-2^-6, 2^-5]);
+// > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-6, 2^-5]);
+// 0x1.0af1...p-34
+//
+// Finally, we compute the product y * log2(x) in FloatFloat:
+// y * log2(x) = y * ((e_x - log2(r)) + log2(1 + dx)).
+// When x in [1 - 2^-5, 1) (where e_x = -1 and R[31] = 0.5 so
+// -log2(R[31]) = 1.0f), e_x - log2(r) cancels to 0, avoiding cancellation.
+// We decompose into integer k and fractional remainder rem:
+// k = round(total_hi),
+// rem = (total_hi - k) + total_lo,
+// and evaluate 2^k * 2^rem via exp2f_eval(rem, k) from exp2f_float_utils.h,
+// which provides a compact, code-size-friendly implementation without large
+// exponent lookup tables.
+
+LIBC_INLINE float powf(float x, float y) {
+ using FloatBits = fputil::FPBits<float>;
+ using namespace powf_internal;
+
+ FloatBits xbits(x), ybits(y);
+ uint32_t x_u = xbits.uintval();
+ uint32_t y_u = ybits.uintval();
+ uint32_t y_a = ybits.abs().uintval();
+
+ if (LIBC_UNLIKELY((x_u & 0x007F'FFFF) == 0 || (y_u & 0x007F'FFFF) == 0)) {
+ if (auto r = powf_internal::check_special_inputs(x, y);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ }
+
+ int e_x = 0;
+ uint32_t sign = 0;
+ if (LIBC_UNLIKELY(y_a <= powf_internal::Y_LOWER_BOUND ||
+ y_a >= powf_internal::Y_UPPER_BOUND ||
+ x_u >= FloatBits::inf().uintval() ||
+ x_u < FloatBits::min_normal().uintval())) {
+ if (auto r = powf_internal::check_exceptional_cases(x, y, e_x, sign);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ }
+
+ Sign out_sign = (sign == 0) ? Sign::POS : Sign::NEG;
+
+ // Extract exponent field and mantissa of x.
+ xbits = FloatBits(x);
+ e_x += xbits.get_biased_exponent() - FloatBits::EXP_BIAS;
+ uint32_t x_mant = xbits.get_mantissa();
+ unsigned idx_x =
+ static_cast<unsigned>(x_mant >> (FloatBits::FRACTION_LEN - 5));
+ FloatBits m_x = FloatBits(x_mant | 0x3f80'0000U);
+
+ // log2(e) in FloatFloat:
+ constexpr fputil::FloatFloat LOG2_E = {0x1.4ae0cp-26f, 0x1.715476p0f};
+
+ // Degree-4 Sollya minimax polynomial for:
+ // Q(dx) ~ (log2(1 + dx) - dx * log2(e)) / dx^2
+ // Generated by Sollya with:
+ // > log2_e = 0x1.715476p0 + 0x1.4ae0cp-26;
+ // > f = (log2(1 + x) - x * log2_e) / x^2;
+ // > Q = fpminimax(f, 4, [|SG...|], [-2^-6, 2^-5]);
+ // > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-6, 2^-5]);
+ // 0x1.0af1...p-34
+ constexpr float COEFFS[] = {-0x1.715476p-1f, 0x1.ec709ep-2f, -0x1.71618ap-2f,
+ 0x1.2810e4p-2f, -0x1.ae42b2p-3f};
+
+ // Range reduction for log2(m_x):
+ // dx = r * m_x - 1 in [-2^-6, 2^-5).
+ // The computation is exact: m_x = (1 + dx) / r.
+ float r = powf_internal::R_32[idx_x];
+#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
+ float dx = fputil::multiply_add(r, m_x.get_val(), -1.0f); // Exact
+#else // !LIBC_TARGET_CPU_HAS_FMA_FLOAT
+ // Without FMA, split m_x = c + cd where c retains the high 14 bits
+ // to perform exact multiplication: r * c - 1.0f and r * cd.
+ float c = FloatBits(m_x.uintval() & 0x3fff'e000U).get_val();
+ float cd = m_x.get_val() - c;
+ float dx_hi = fputil::multiply_add(r, c, -1.0f);
+ float dx = fputil::multiply_add(r, cd, dx_hi); // Exact
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
+
+ fputil::FloatFloat log2_r = powf_internal::LOG2_R_FF_32[idx_x];
+ float e_xf = static_cast<float>(e_x);
+ // (e_xf + log2_r.hi) is exact since log2_r.hi has 16 bits of precision and
+ // |e_x| <= 150 < 2^8.
+ float e_x_hi = e_xf + log2_r.hi;
+
+ // Early overflow and underflow check:
+ // z = y * log2(x) ~ y * (e_x_hi + dx * log2(e))
+ // With |dx| < 2^-5, |log2(1 + dx) - dx * log2(e)| <= (dx^2 / 2) * log2(e) <
+ // 2^-10. For |y * log2(x)| <= 128, the error |z - y * log2(x)| < 4. Hence:
+ // z > 132 => overflow
+ // z < -155 => underflow
+ float z = y * fputil::multiply_add(dx, LOG2_E.hi, e_x_hi);
+
+ if (LIBC_UNLIKELY(z > 132.0f))
+ return set_overflow(out_sign);
+
+ if (LIBC_UNLIKELY(z < -155.0f))
+ return set_underflow(out_sign);
+
+ // Approximate the tail Q(dx) using Estrin's scheme:
+ float dx2 = dx * dx;
+ float dx4 = dx2 * dx2;
+
+ float p0 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
+ float p1 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
+
+ float q0 = fputil::multiply_add(dx2, p1, p0);
+ float p = fputil::multiply_add(dx4, COEFFS[4], q0);
+ float p_tail = dx2 * p;
+
+ // Compute dx * log2(e) in FloatFloat.
+ fputil::FloatFloat dx_log2_e = fputil::exact_mult<float>(dx, LOG2_E.hi);
+ float dx_log2_e_lo = fputil::multiply_add(dx, LOG2_E.lo, dx_log2_e.lo);
+
+ // Combine with log2(1 + dx) in FloatFloat:
+ fputil::FloatFloat log2_x_hi =
+ fputil::exact_add<false, float>(e_x_hi, dx_log2_e.hi);
+ float log2_x_lo = (dx_log2_e_lo + p_tail) + log2_r.lo + log2_x_hi.lo;
+ fputil::FloatFloat log2_x =
+ fputil::exact_add<false, float>(log2_x_hi.hi, log2_x_lo);
+
+ fputil::FloatFloat s = fputil::exact_mult<float>(y, log2_x.hi);
+ float total_hi = s.hi;
+ float total_lo = fputil::multiply_add(y, log2_x.lo, s.lo);
+
+ // Exponent reduction for 2^(total_hi + total_lo):
+ // k = round(total_hi)
+ // rem = (total_hi - k) + total_lo, with |rem| <= 0.5.
+ float kf = fputil::nearest_integer(total_hi);
+ int64_t k = static_cast<int64_t>(kf);
+
+ float rem_hi = total_hi - kf;
+ float rem = rem_hi + total_lo;
+
+ // Normal range: -125 <= k <= 127.
+ uint64_t k_u = static_cast<uint64_t>(k + 125);
+ if (LIBC_LIKELY(k_u <= 252)) {
+ float result = exp2f_eval(rem, static_cast<int>(k));
+ return cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
+ }
+
+ // Exceptional path: k > 128 (overflow), k < -150 (underflow), or subnormal.
+ if (k > 128)
+ return set_overflow(out_sign);
+
+ if (k < -150)
+ return set_underflow(out_sign);
+
+ // Subnormal range: -150 <= k <= -126.
+ if (k <= -126) {
+ float result = exp2f_eval(rem, static_cast<int>(k));
+ float res = cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
+
+ if (LIBC_UNLIKELY(FloatBits(res).is_normal()))
+ return res;
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<float>();
+ return res;
+ }
+
+ // Boundary case: k = 128.
+ float result = exp2f_eval(rem, static_cast<int>(k));
+ if (LIBC_UNLIKELY(FloatBits(result).is_inf()))
+ return set_overflow(out_sign);
+
+ return cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
+}
+
+} // namespace float_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_FLOAT_EVAL_H
diff --git a/libc/src/__support/math/powf_small_tables.h b/libc/src/__support/math/powf_small_tables.h
deleted file mode 100644
index 3c5abbe9190ce6..00000000000000
--- a/libc/src/__support/math/powf_small_tables.h
+++ /dev/null
@@ -1,120 +0,0 @@
-//===-- Implementation header for powf using less memory --------*- C++ -*-===//
-//
-// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
-// See https://llvm.org/LICENSE.txt for license information.
-// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
-//
-//===----------------------------------------------------------------------===//
-
-#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_SMALL_TABLES_H
-#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_SMALL_TABLES_H
-
-#include "src/__support/CPP/bit.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/PolyEval.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"
-
-namespace LIBC_NAMESPACE_DECL {
-
-namespace math {
-
-namespace powf_internal {
-
-LIBC_INLINE LIBC_CONSTEXPR float powf_small_tables(float x, int ex,
- uint64_t sign, float y) {
- using FloatBits = fputil::FPBits<float>;
- using DoubleBits = fputil::FPBits<double>;
-
- constexpr double ONE_OVER_SQRT2 = 0x1.6a09e667f3bcdp-1;
-
- // x^y = 2^( y * log2(x) )
- // = 2^( y * ( e_x + log2(m_x) ) )
- // First we compute log2(x) = e_x + log2(m_x)
- uint32_t x_u = FloatBits(x).uintval();
-
- double yd = static_cast<double>(y);
-
- // Extract exponent field of x.
- ex += (x_u >> FloatBits::FRACTION_LEN);
- double e_x = static_cast<double>(ex);
-
- // Add the hidden bit to the mantissa.
- // 1 <= m_x < 2
- uint32_t x_mant = (x_u & FloatBits::FRACTION_MASK);
- double m_x = static_cast<double>(cpp::bit_cast<float>(x_mant | 0x3f800000));
- // Reduce to 1 <= mx <= sqrt(2).
- if (x_mant > 0x0045'04f3) {
- e_x += 0.5;
- m_x *= ONE_OVER_SQRT2;
- }
- // 0 <= dx <= sqrt(2) - 1.
- double dx = m_x - 1.0;
-
- // Degree-13 polynomial approximation:
- // dx * P(dx) ~ log2(1 + dx)
- // Generated by Sollya with:
- // > P = fpminimax(log2(1 + x)/x, 13, [|D...|], [0, sqrt(2) - 1]);
- // > dirtyinfnorm((log2(1 + x) - x*P)/log2(1 + x), [0, sqrt(2) - 1]);
- // 0x1.b2d...p-53
- constexpr double LOG2_COEFFS[] = {
- 0x1.71547652b82fdp0, -0x1.71547652b7a2ap-1, 0x1.ec709dc2edfa6p-2,
- -0x1.71547626a9d98p-2, 0x1.2776bf5f6f40ep-2, -0x1.ec6fbbf289ce3p-3,
- 0x1.a60bf904470a7p-3, -0x1.70ef61b01fc1ep-3, 0x1.45d3270454507p-3,
- -0x1.1c5fc05b06e8fp-3, 0x1.d0f57944a937fp-4, -0x1.413e22be24d32p-4,
- 0x1.3c84b66491ccp-5, -0x1.3df9cfe5e602ep-7};
-
- double dx2 = dx * dx;
- double c0 = fputil::multiply_add(dx, LOG2_COEFFS[1], LOG2_COEFFS[0]);
- double c1 = fputil::multiply_add(dx, LOG2_COEFFS[3], LOG2_COEFFS[2]);
- double c2 = fputil::multiply_add(dx, LOG2_COEFFS[5], LOG2_COEFFS[4]);
- double c3 = fputil::multiply_add(dx, LOG2_COEFFS[7], LOG2_COEFFS[6]);
- double c4 = fputil::multiply_add(dx, LOG2_COEFFS[9], LOG2_COEFFS[8]);
- double c5 = fputil::multiply_add(dx, LOG2_COEFFS[11], LOG2_COEFFS[10]);
- double c6 = fputil::multiply_add(dx, LOG2_COEFFS[13], LOG2_COEFFS[12]);
-
- double dx4 = dx2 * dx2;
- double d0 = fputil::multiply_add(dx2, c1, c0);
- double d1 = fputil::multiply_add(dx2, c3, c2);
- double d2 = fputil::multiply_add(dx2, c5, c4);
-
- double p = fputil::polyeval(dx4, d0, d1, d2, c6);
- // u ~ y * log2(x).
- double u = yd * fputil::multiply_add(dx, p, e_x);
-
- double hi = fputil::nearest_integer(u);
- double lo = u - hi;
- int e_hi = static_cast<int>(hi) + DoubleBits::EXP_BIAS;
- double exp_hi = cpp::bit_cast<double>(
- (static_cast<uint64_t>(e_hi) << DoubleBits::FRACTION_LEN) | sign);
- // Degree-6 polynomial approximation P(lo6) ~ 2^(lo6 / 2^6) = 2^(lo).
- // Generated by Sollya with:
- // > P = fpminimax(2^x, 6, [|1, D...|], [-0.5, 0.5]);
- // > dirtyinfnorm(2^x - P, [-0.5, 0.5]);
- // 0x1.5f7...p-29
- constexpr double EXP2_COEFFS[] = {
- 0x1.62e430c7b13a8p-1, 0x1.ebfbdd2f82f6fp-3, 0x1.c6aed4f186f34p-5,
- 0x1.3b2c96c9aa336p-7, 0x1.5f4553ff53f9p-10, 0x1.4278e5fa9de78p-13};
-
- double lo2 = lo * lo;
- double f0 = fputil::multiply_add(lo, EXP2_COEFFS[1], EXP2_COEFFS[0]);
- double f1 = fputil::multiply_add(lo, EXP2_COEFFS[3], EXP2_COEFFS[2]);
- double f2 = fputil::multiply_add(lo, EXP2_COEFFS[5], EXP2_COEFFS[4]);
-
- double pp = fputil::polyeval(lo2, f0, f1, f2);
-
- double r = fputil::multiply_add(lo, pp, 1.0);
-
- double result = r * exp_hi;
-
- return static_cast<float>(result);
-}
-
-} // namespace powf_internal
-} // namespace math
-} // namespace LIBC_NAMESPACE_DECL
-
-#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_SMALL_TABLES_H
diff --git a/libc/src/__support/math/powf_utils.h b/libc/src/__support/math/powf_utils.h
new file mode 100644
index 00000000000000..46f113b6ff0373
--- /dev/null
+++ b/libc/src/__support/math/powf_utils.h
@@ -0,0 +1,718 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 single-precision powf(x, y).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_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/triple_double.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"
+
+#if defined(LIBC_TARGET_CPU_HAS_FPU_FLOAT) || \
+ !defined(LIBC_MATH_HAS_SMALL_TABLES)
+#include "src/__support/FPUtil/sqrt.h"
+#endif // LIBC_TARGET_CPU_HAS_FPU_FLOAT || !LIBC_MATH_HAS_SMALL_TABLES
+
+#ifndef LIBC_MATH_HAS_SMALL_TABLES
+#include "src/__support/math/exp10f.h"
+#include "src/__support/math/exp2f.h"
+#endif // !LIBC_MATH_HAS_SMALL_TABLES
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace powf_internal {
+
+using fputil::DoubleDouble;
+using fputil::FloatFloat;
+using fputil::TripleDouble;
+
+// Table of -log2(R[i]) represented in triple-double precision {lo, mid, hi},
+// where R[i] = 2^-8 * ceil(2^8 * (1 - 2^-8) / (1 + i * 2^-7)) for i = 0..127.
+//
+// The least significant bit of the high part `hi` is chosen to be >= 2^-20,
+// so that the total precision of (e_x + hi) is at most 8 + 20 = 28 bits
+// (since |e_x| <= 150 < 2^8 fits in 8 bits).
+// Because the single-precision exponent y has 24 bits of mantissa, the product
+// y * (e_x + hi) spans at most 24 + 28 = 52 bits <= 53 bits, making it exact
+// in double precision.
+// The mid part provides 53 bits for fast evaluation, and the lo part provides
+// another 53 bits for the accurate fallback, yielding >= 127 bits of precision.
+//
+// Generated by Sollya with:
+// > display = hexadecimal;
+// > prec = 200;
+// > for i from 0 to 127 do {
+// if (i == 0) then {
+// print("{0.0, 0.0, 0.0},");
+// } else if (i == 127) then {
+// print("{0.0, 0.0, 1.0},");
+// } else {
+// r = 2^(-8) * ceil( 2^8 * (1 - 2^(-8)) / (1 + i*2^(-7)) );
+// l = -log2(r);
+// h = nearestint(l * 2^20) * 2^(-20);
+// m = round(l - h, D, RN);
+// lo = round(l - h - m, D, RN);
+// print("{", lo, ", ", m, ", ", h, "},");
+// };
+// };
+LIBC_INLINE_VAR constexpr TripleDouble LOG2_R_TD[128] = {
+ {0.0, 0.0, 0.0},
+ {0x1.84a2c615b70adp-79, -0x1.177c23362928cp-25, 0x1.72c8p-7},
+ {-0x1.f27b820fd03eap-76, -0x1.179e0caa9c9abp-22, 0x1.744p-6},
+ {-0x1.f27ef487c8f34p-77, -0x1.c6cea541f5b7p-23, 0x1.184cp-5},
+ {-0x1.e3f80fbc71454p-76, -0x1.66c4d4e554434p-22, 0x1.773ap-5},
+ {-0x1.9f8ef14d5f6eep-79, -0x1.70700a00fdd55p-24, 0x1.d6ecp-5},
+ {0x1.452bbce7398c1p-77, 0x1.53002a4e86631p-23, 0x1.1bb3p-4},
+ {-0x1.990555535afdp-81, 0x1.fcd15f101c142p-25, 0x1.4c56p-4},
+ {0x1.447e30ad393eep-78, 0x1.25b3eed319cedp-22, 0x1.7d6p-4},
+ {0x1.b7759da88a2dap-76, -0x1.4195120d8486fp-22, 0x1.960dp-4},
+ {0x1.cee7766ece702p-78, 0x1.45b878e27d0d9p-23, 0x1.c7b5p-4},
+ {-0x1.a55c745ecdc2fp-77, 0x1.770744593a4cbp-22, 0x1.f9c9p-4},
+ {0x1.f7ec992caa67fp-77, 0x1.c673032495d24p-22, 0x1.097ep-3},
+ {-0x1.433638c6ece3ep-77, -0x1.1eaa65b49696ep-22, 0x1.22dbp-3},
+ {0x1.58f27b6518824p-76, 0x1.b2866f2850b22p-22, 0x1.3c6f8p-3},
+ {-0x1.86bdcfdfd4a4cp-79, 0x1.8ee37cd2ea9d3p-25, 0x1.494f8p-3},
+ {-0x1.ff7044a68a7fap-80, 0x1.7e86f9c2154fbp-24, 0x1.633a8p-3},
+ {-0x1.aa21694561327p-81, 0x1.8e3cfc25f0ce6p-26, 0x1.7046p-3},
+ {-0x1.d209f2d4239c6p-87, 0x1.57f7a64ccd537p-28, 0x1.8a898p-3},
+ {-0x1.a55e97e60e632p-76, -0x1.a761c09fbd2aep-22, 0x1.97c2p-3},
+ {0x1.261179225541ep-76, 0x1.24bea9a2c66f3p-22, 0x1.b26p-3},
+ {-0x1.08fa30510fca9p-82, -0x1.60002ccfe43f5p-25, 0x1.bfc68p-3},
+ {-0x1.63ec8d56242f9p-76, 0x1.69f220e97f22cp-22, 0x1.dac2p-3},
+ {0x1.8bcdaf0534365p-76, -0x1.6164f64c210ep-22, 0x1.e858p-3},
+ {0x1.1003282896056p-78, -0x1.0c1678ae89767p-24, 0x1.01d9cp-2},
+ {0x1.01bcc7025fa92p-78, -0x1.f26a05c813d57p-22, 0x1.08bdp-2},
+ {-0x1.fe8a8648e9ebcp-80, 0x1.4d8fc561c8d44p-24, 0x1.169cp-2},
+ {0x1.08dfb23650c75p-79, -0x1.362ad8f7ca2dp-22, 0x1.1d984p-2},
+ {-0x1.f8d5a89861a5ep-79, 0x1.2b13cd6c4d042p-22, 0x1.249ccp-2},
+ {-0x1.a1c872983511ep-76, -0x1.1c8f11979a5dbp-22, 0x1.32cp-2},
+ {0x1.e8e21bff3336bp-77, 0x1.c2ab3edefe569p-23, 0x1.39de8p-2},
+ {0x1.fd1994fb2c4a1p-80, 0x1.7c3eca28e69cap-26, 0x1.4106p-2},
+ {0x1.6b94b51cf76b1p-80, -0x1.34c4e99e1c6c6p-24, 0x1.4f6fcp-2},
+ {-0x1.31d55da1d0f66p-76, -0x1.194a871b63619p-22, 0x1.56b24p-2},
+ {-0x1.378b22691e28bp-77, 0x1.e3dd5c1c885aep-23, 0x1.5dfdcp-2},
+ {0x1.99e302970e411p-83, -0x1.6ccf3b1129b7cp-23, 0x1.6552cp-2},
+ {0x1.20164a049664dp-82, -0x1.2f346e2bf924bp-23, 0x1.6cb1p-2},
+ {-0x1.d14aac4d864c3p-77, -0x1.fa61aaa59c1d8p-23, 0x1.7b8ap-2},
+ {0x1.496ab4e4b293fp-79, 0x1.90c11fd32a3abp-22, 0x1.8304cp-2},
+ {-0x1.d209f2d4239c6p-86, 0x1.57f7a64ccd537p-27, 0x1.8a898p-2},
+ {0x1.eae3326327babp-81, 0x1.249ba76fee235p-27, 0x1.9218p-2},
+ {0x1.fa05bddfded8cp-77, -0x1.aad2729b21ae5p-23, 0x1.99b08p-2},
+ {-0x1.624140d175ba2p-77, 0x1.71810a5e1818p-22, 0x1.a8ff8p-2},
+ {0x1.f1c5160c515c1p-81, -0x1.6172fe015e13cp-27, 0x1.b0b68p-2},
+ {-0x1.86a6204eec8cp-79, 0x1.5ec6c1bfbf89ap-24, 0x1.b877cp-2},
+ {0x1.718f761dd3915p-78, 0x1.678bf6cdedf51p-24, 0x1.c0438p-2},
+ {-0x1.d4ee66c3700e4p-76, 0x1.c2d45fe43895ep-22, 0x1.c819cp-2},
+ {-0x1.7d14533586306p-77, -0x1.9ee52ed49d71dp-22, 0x1.cffbp-2},
+ {0x1.5ce9fb5a7bb5bp-81, 0x1.5786af187a96bp-27, 0x1.d7e6cp-2},
+ {-0x1.ae6face57ad3bp-77, 0x1.3ab0dc56138c9p-23, 0x1.dfdd8p-2},
+ {0x1.5ac93b443d55fp-78, 0x1.fe538ab34efb5p-22, 0x1.e7df4p-2},
+ {0x1.f1753e0ae1e8fp-76, -0x1.e4fee07aa4b68p-22, 0x1.efec8p-2},
+ {0x1.cdfd4c297069bp-76, -0x1.172f32fe67287p-22, 0x1.f804cp-2},
+ {0x1.97a0e8f3ba742p-79, -0x1.9a83ff9ab9cc8p-22, 0x1.00144p-1},
+ {-0x1.800450f5b2357p-78, -0x1.68cb06cece193p-22, 0x1.042bep-1},
+ {-0x1.a839041241fe7p-78, 0x1.8cd71ddf82e2p-22, 0x1.08494p-1},
+ {0x1.ed0b8eeccca86p-78, 0x1.5e18ab2df3ae6p-22, 0x1.0c6cap-1},
+ {0x1.3dd41df9689b3p-79, 0x1.5dee4d9d8a273p-25, 0x1.1096p-1},
+ {-0x1.990555535afdp-82, 0x1.fcd15f101c142p-26, 0x1.14c56p-1},
+ {-0x1.1773d02c9055cp-77, -0x1.2474b0f992ba1p-23, 0x1.18faep-1},
+ {-0x1.4aeef330c53c1p-78, 0x1.4b5a92a606047p-24, 0x1.1d368p-1},
+ {0x1.8e6ff749ebacbp-77, 0x1.16186fcf54bbdp-22, 0x1.21786p-1},
+ {0x1.c09d761c548ebp-84, 0x1.18efabeb7d722p-27, 0x1.25c0ap-1},
+ {0x1.aaa73a428e1e4p-78, -0x1.e5fc7d238691dp-24, 0x1.2a0f4p-1},
+ {-0x1.af2f3d8b63fbap-79, 0x1.f5809faf6283cp-22, 0x1.2e644p-1},
+ {-0x1.af2f3d8b63fbap-79, 0x1.f5809faf6283cp-22, 0x1.2e644p-1},
+ {0x1.78de359f2bb88p-77, 0x1.c6e1dcd0cb449p-22, 0x1.32bfep-1},
+ {-0x1.415ae1a715618p-76, 0x1.76e0e8f74b4d5p-22, 0x1.37222p-1},
+ {-0x1.4991b5375621fp-79, -0x1.cb82c89692d99p-24, 0x1.3b8b2p-1},
+ {-0x1.827d37deb2236p-76, -0x1.63161c5432aebp-22, 0x1.3ffaep-1},
+ {0x1.9576edac01c78p-77, 0x1.458104c41b901p-22, 0x1.44716p-1},
+ {0x1.9576edac01c78p-77, 0x1.458104c41b901p-22, 0x1.44716p-1},
+ {-0x1.05a27b81e2219p-77, -0x1.cd9d0cde578d5p-22, 0x1.48efp-1},
+ {0x1.237616778b4bap-82, 0x1.b9884591add87p-26, 0x1.4d738p-1},
+ {0x1.3b7d7e5d148bbp-76, 0x1.c6042978605ffp-22, 0x1.51ff2p-1},
+ {-0x1.cc3f936a5977cp-79, -0x1.fc4c96b37dcf6p-22, 0x1.56922p-1},
+ {0x1.20164a049664dp-83, -0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
+ {0x1.20164a049664dp-83, -0x1.2f346e2bf924bp-24, 0x1.5b2c4p-1},
+ {-0x1.a212919a92f7ap-77, 0x1.c4e4fbb68a4d1p-22, 0x1.5fcdcp-1},
+ {-0x1.b64b03f7230ddp-77, -0x1.9d499bd9b3226p-23, 0x1.6476ep-1},
+ {-0x1.1ec6379e6e3b9p-77, -0x1.f89b355ede26fp-23, 0x1.69278p-1},
+ {-0x1.1ec6379e6e3b9p-77, -0x1.f89b355ede26fp-23, 0x1.69278p-1},
+ {-0x1.4ba44c03bfbbdp-78, 0x1.53c7e319f6e92p-24, 0x1.6ddfcp-1},
+ {-0x1.c36fc650d030fp-77, -0x1.b291f070528c7p-22, 0x1.729fep-1},
+ {-0x1.69e5693a7f067p-80, 0x1.2967a451a7b48p-25, 0x1.7767cp-1},
+ {-0x1.69e5693a7f067p-80, 0x1.2967a451a7b48p-25, 0x1.7767cp-1},
+ {0x1.6598aae91499ap-76, 0x1.244fcff690fcep-22, 0x1.7c37ap-1},
+ {0x1.99d61ec432837p-77, 0x1.46fd97f5dc572p-23, 0x1.810fap-1},
+ {0x1.99d61ec432837p-77, 0x1.46fd97f5dc572p-23, 0x1.810fap-1},
+ {0x1.855c42078f81bp-76, -0x1.f3a7352663e5p-22, 0x1.85efep-1},
+ {-0x1.59408e815107p-77, 0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
+ {-0x1.59408e815107p-77, 0x1.b3cda690370b5p-23, 0x1.8ad84p-1},
+ {0x1.33b318085e50ap-78, 0x1.3226b211bf1d9p-23, 0x1.8fc92p-1},
+ {0x1.343fe7c9cb4aep-79, 0x1.d24b136c101eep-23, 0x1.94c28p-1},
+ {0x1.343fe7c9cb4aep-79, 0x1.d24b136c101eep-23, 0x1.94c28p-1},
+ {-0x1.d19522e56fe6p-76, 0x1.7c40c7907e82ap-22, 0x1.99c48p-1},
+ {-0x1.23b9d8ea55c3ep-77, -0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
+ {-0x1.23b9d8ea55c3ep-77, -0x1.e81781d97ee91p-22, 0x1.9ecf6p-1},
+ {0x1.829440c24aeb6p-78, -0x1.6a77813f94e01p-22, 0x1.a3e3p-1},
+ {-0x1.624140d175ba2p-76, -0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
+ {-0x1.624140d175ba2p-76, -0x1.1cfdeb43cfdp-22, 0x1.a8ffap-1},
+ {0x1.afa6f024fb045p-77, -0x1.f983f74d3138fp-23, 0x1.ae256p-1},
+ {-0x1.603ad3a5d326dp-78, -0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
+ {-0x1.603ad3a5d326dp-78, -0x1.e278ae1a1f51fp-23, 0x1.b3546p-1},
+ {-0x1.0c1e0e5855d6ap-77, -0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
+ {-0x1.0c1e0e5855d6ap-77, -0x1.97552b7b5ea45p-23, 0x1.b88ccp-1},
+ {0x1.c817ad56baa16p-78, -0x1.19b4f3c72c4f8p-24, 0x1.bdceap-1},
+ {0x1.44c47ac1bf62bp-77, 0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
+ {0x1.44c47ac1bf62bp-77, 0x1.f7402d26f1a12p-23, 0x1.c31a2p-1},
+ {-0x1.69b9465eae1e6p-78, -0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
+ {-0x1.69b9465eae1e6p-78, -0x1.2056d5dd31d96p-23, 0x1.c86f8p-1},
+ {-0x1.24a6d9d1d1904p-79, -0x1.6e46335aae723p-24, 0x1.cdcecp-1},
+ {-0x1.3826144575ac4p-76, -0x1.beb244c59f331p-22, 0x1.d3382p-1},
+ {-0x1.3826144575ac4p-76, -0x1.beb244c59f331p-22, 0x1.d3382p-1},
+ {0x1.dbc96b3b12b25p-81, 0x1.16c071e93fd97p-27, 0x1.d8abap-1},
+ {0x1.dbc96b3b12b25p-81, 0x1.16c071e93fd97p-27, 0x1.d8abap-1},
+ {0x1.68a8ccdbd1f33p-77, 0x1.d8175819530c2p-22, 0x1.de298p-1},
+ {0x1.68a8ccdbd1f33p-77, 0x1.d8175819530c2p-22, 0x1.de298p-1},
+ {0x1.e586711df5ea1p-79, 0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
+ {0x1.e586711df5ea1p-79, 0x1.51bd552842c1cp-23, 0x1.e3b2p-1},
+ {-0x1.bc25adf042483p-79, 0x1.914e204f19d94p-22, 0x1.e9452p-1},
+ {-0x1.bc25adf042483p-79, 0x1.914e204f19d94p-22, 0x1.e9452p-1},
+ {0x1.d7d82b65c5686p-76, 0x1.c55d997da24fdp-22, 0x1.eee32p-1},
+ {0x1.d7d82b65c5686p-76, 0x1.c55d997da24fdp-22, 0x1.eee32p-1},
+ {-0x1.3f108c0857ca3p-77, -0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
+ {-0x1.3f108c0857ca3p-77, -0x1.685c2d2298a6ep-22, 0x1.f48c4p-1},
+ {-0x1.bd800bca7a221p-78, 0x1.7a4887bd74039p-22, 0x1.fa406p-1},
+ {0.0, 0.0, 1.0},
+};
+
+// Lookup table for 5-bit (32 entries) range reduction for float_eval and
+// double_eval fast path:
+// r(0) = 1.0f
+// r(i) = 2^-6 * ceil(2^6 * (1 - 2^-6) / (1 + i * 2^-5)), i = 1..31.
+// The constants are chosen so that dx = r * m_x - 1 is exact in
+// single precision (with hardware FMA or 14-bit Sterbenz split), and
+// -2^-6 <= dx < 2^-5.
+//
+// Generated by Sollya with:
+// > display = hexadecimal;
+// > for i from 0 to 31 do {
+// if (i == 0) then {
+// print("0x1.0p+0f,");
+// } else {
+// r = 2^(-6) * ceil( 2^6 * (1 - 2^(-6)) / (1 + i*2^(-5)) );
+// print(r @ "f,");
+// };
+// };
+LIBC_INLINE_VAR constexpr float R_32[32] = {
+ 0x1.0p+0f, 0x1.fp-1f, 0x1.ep-1f, 0x1.dp-1f, 0x1.cp-1f, 0x1.b8p-1f,
+ 0x1.bp-1f, 0x1.ap-1f, 0x1.98p-1f, 0x1.9p-1f, 0x1.8p-1f, 0x1.78p-1f,
+ 0x1.7p-1f, 0x1.68p-1f, 0x1.6p-1f, 0x1.58p-1f, 0x1.5p-1f, 0x1.5p-1f,
+ 0x1.48p-1f, 0x1.4p-1f, 0x1.38p-1f, 0x1.38p-1f, 0x1.3p-1f, 0x1.28p-1f,
+ 0x1.2p-1f, 0x1.2p-1f, 0x1.18p-1f, 0x1.18p-1f, 0x1.1p-1f, 0x1.1p-1f,
+ 0x1.08p-1f, 0x1.0p-1f,
+};
+
+// R_32 represented in double precision for double_eval.
+LIBC_INLINE_VAR constexpr double R_32_D[32] = {
+ 0x1.0p+0, 0x1.fp-1, 0x1.ep-1, 0x1.dp-1, 0x1.cp-1, 0x1.b8p-1, 0x1.bp-1,
+ 0x1.ap-1, 0x1.98p-1, 0x1.9p-1, 0x1.8p-1, 0x1.78p-1, 0x1.7p-1, 0x1.68p-1,
+ 0x1.6p-1, 0x1.58p-1, 0x1.5p-1, 0x1.5p-1, 0x1.48p-1, 0x1.4p-1, 0x1.38p-1,
+ 0x1.38p-1, 0x1.3p-1, 0x1.28p-1, 0x1.2p-1, 0x1.2p-1, 0x1.18p-1, 0x1.18p-1,
+ 0x1.1p-1, 0x1.1p-1, 0x1.08p-1, 0x1.0p-1,
+};
+
+// Table of -log2(R_32[i]) represented in double precision for double_eval.
+//
+// Generated by Sollya with:
+// > display = hexadecimal;
+// > prec = 200;
+// > for i from 0 to 31 do {
+// if (i == 0) then {
+// print("0x0.0000000000000p+0,");
+// } else if (i == 31) then {
+// print("0x1.0000000000000p+0,");
+// } else {
+// r = 2^(-6) * ceil( 2^6 * (1 - 2^(-6)) / (1 + i*2^(-5)) );
+// print(round(-log2(r), D, RN), ",");
+// };
+// };
+LIBC_INLINE_VAR constexpr double LOG2_R_32[32] = {
+ 0x0.0000000000000p+0, 0x1.77394c9d958d5p-5, 0x1.7d60496cfbb4cp-4,
+ 0x1.22dadc2ab3497p-3, 0x1.8a8980abfbd32p-3, 0x1.bfc67a7fff4ccp-3,
+ 0x1.f5fd8a9063e35p-3, 0x1.32bfee370ee68p-2, 0x1.4f6fbb2cec598p-2,
+ 0x1.6cb0f6865c8eap-2, 0x1.a8ff971810a5ep-2, 0x1.c819dc2d45fe4p-2,
+ 0x1.e7df5fe538ab3p-2, 0x1.042bd4b9a7c99p-1, 0x1.14c560fe68af9p-1,
+ 0x1.25c0a0463bebp-1, 0x1.37222bb70747cp-1, 0x1.37222bb70747cp-1,
+ 0x1.48eef19317991p-1, 0x1.5b2c3da19723bp-1, 0x1.6ddfc2a78fc63p-1,
+ 0x1.6ddfc2a78fc63p-1, 0x1.810fa51bf65fdp-1, 0x1.94c287492c4dbp-1,
+ 0x1.a8ff971810a5ep-1, 0x1.a8ff971810a5ep-1, 0x1.bdce9dcc96187p-1,
+ 0x1.bdce9dcc96187p-1, 0x1.d338120a6dd9dp-1, 0x1.d338120a6dd9dp-1,
+ 0x1.e9452c8a71028p-1, 0x1.0000000000000p+0,
+};
+
+// Table of -log2(R_32[i]) represented in FloatFloat precision {lo, hi} for
+// float_eval.
+// We choose the precision of the high part to be 24 - 8 = 16 bits, so that
+// e_xf + LOG2_R_FF_32[i].hi
+// is exact in float for |e_x| <= 150.
+// For i = 31, R_32[31] = 0.5, so -log2(R_32[31]) = 1.0f is exact.
+//
+// Generated by Sollya with:
+// > display = hexadecimal;
+// > prec = 200;
+// > for i from 0 to 31 do {
+// if (i == 0) then {
+// print("{0x0.000000p+0f, 0x0.000000p+0f},");
+// } else if (i == 31) then {
+// print("{0x0.000000p+0f, 0x1.000000p+0f},");
+// } else {
+// r = 2^(-6) * ceil( 2^6 * (1 - 2^(-6)) / (1 + i*2^(-5)) );
+// l = -log2(r);
+// h = round(1 + l, 17, RN) - 1;
+// lo = round(l - h, SG, RN);
+// print("{" @ lo @ "f, " @ h @ "f},");
+// };
+// };
+LIBC_INLINE_VAR constexpr FloatFloat LOG2_R_FF_32[32] = {
+ {0x0.000000p+0f, 0x0.000000p+0f}, {-0x1.acd89ap-19f, 0x1.774p-5f},
+ {0x1.25b3eep-22f, 0x1.7d6p-4f}, {0x1.6e155ap-18f, 0x1.22d8p-3f},
+ {0x1.80abfcp-19f, 0x1.8a88p-3f}, {-0x1.858p-19f, 0x1.bfc8p-3f},
+ {-0x1.3ab7cep-18f, 0x1.f6p-3f}, {-0x1.1c8f12p-22f, 0x1.32cp-2f},
+ {-0x1.134c4ep-20f, 0x1.4f7p-2f}, {0x1.ed0cbap-19f, 0x1.6cbp-2f},
+ {-0x1.a39fbep-20f, 0x1.a9p-2f}, {0x1.dc2d46p-18f, 0x1.c818p-2f},
+ {-0x1.40358ep-19f, 0x1.e7ep-2f}, {-0x1.5a32c2p-20f, 0x1.042cp-1f},
+ {-0x1.3e032ep-18f, 0x1.14c6p-1f}, {0x1.408c78p-18f, 0x1.25c0p-1f},
+ {0x1.5db83ap-20f, 0x1.3722p-1f}, {0x1.5db83ap-20f, 0x1.3722p-1f},
+ {0x1.e3263p-18f, 0x1.48eep-1f}, {0x1.ed0cbap-20f, 0x1.5b2cp-1f},
+ {-0x1.eac382p-20f, 0x1.6de0p-1f}, {-0x1.eac382p-20f, 0x1.6de0p-1f},
+ {-0x1.6b9026p-19f, 0x1.811p-1f}, {0x1.0e9258p-18f, 0x1.94c2p-1f},
+ {-0x1.a39fbep-19f, 0x1.a9p-1f}, {-0x1.a39fbep-19f, 0x1.a9p-1f},
+ {0x1.3b992cp-18f, 0x1.bdcep-1f}, {0x1.3b992cp-18f, 0x1.bdcep-1f},
+ {0x1.20a6dep-21f, 0x1.d338p-1f}, {0x1.20a6dep-21f, 0x1.d338p-1f},
+ {-0x1.a6eb1ep-18f, 0x1.e946p-1f}, {0x0.000000p+0f, 0x1.000000p+0f},
+};
+
+// Upper bound for y = 150 / |log2(1 - 2^-24)|, generated by Sollya:
+// > y = round(-150 / log2(1 - 2^(-24)), SG, RU);
+// > y;
+// 0x1.9fe368p30
+// > printsingle(y);
+// 0x4ecff1b4
+constexpr uint32_t Y_UPPER_BOUND = 0x4ecf'f1b4;
+// Lower bound for y = 2^-25 / 150, generated by Sollya:
+// > y = round(2^(-25) / 150, SG, RD);
+// > y;
+// 0x1.b4e81ap-33
+// > printsingle(y);
+// 0x2f5a740d
+constexpr uint32_t Y_LOWER_BOUND = 0x2f5a'740d;
+
+// 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(float x) {
+ using FPBits = fputil::FPBits<float>;
+ FPBits xbits(x);
+ uint32_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(float x) {
+ if (x == 0.0f)
+ return true;
+ using FPBits = fputil::FPBits<float>;
+ FPBits xbits(x);
+ uint32_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 float set_overflow(Sign sign = Sign::POS) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_overflow_except_if_required<float>();
+ using FPBits = fputil::FPBits<float>;
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_TOWARDZERO)
+ return FPBits::max_normal(sign).get_val();
+ if (rounding == FE_DOWNWARD)
+ return sign.is_neg() ? FPBits::inf(Sign::NEG).get_val()
+ : FPBits::max_normal(Sign::POS).get_val();
+ if (rounding == FE_UPWARD)
+ return sign.is_neg() ? FPBits::max_normal(Sign::NEG).get_val()
+ : FPBits::inf(Sign::POS).get_val();
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return FPBits::inf(sign).get_val();
+}
+
+LIBC_INLINE float set_underflow(Sign sign = Sign::POS) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<float>();
+ using FPBits = fputil::FPBits<float>;
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_UPWARD && sign.is_pos())
+ return FPBits::min_subnormal(Sign::POS).get_val();
+ if (rounding == FE_DOWNWARD && sign.is_neg())
+ return FPBits::min_subnormal(Sign::NEG).get_val();
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return FPBits::zero(sign).get_val();
+}
+
+// Fast checks for special inputs:
+// x = 0, +-1, 2^k, 10, +- Inf
+// y = 0, +-1, 2, 0.5
+LIBC_ALWAYS_INLINE cpp::optional<float> check_special_inputs(float x, float y) {
+ using FPBits = fputil::FPBits<float>;
+ 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();
+
+ // 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.0f || y == 0.0f)
+ return 1.0f;
+
+ if (x == 0.0f) {
+ if (y_abs.is_nan())
+ return y;
+ if (y_abs.is_inf())
+ return y_sign ? FPBits::inf().get_val() : 0.0f;
+ 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.0f : 0.0f;
+ }
+
+ if (y == 1.0f)
+ return x;
+
+ if (y == 2.0f)
+ return x * x;
+
+#ifdef LIBC_TARGET_CPU_HAS_FPU_FLOAT
+ if (y == 0.5f && !x_sign)
+ return fputil::sqrt<float>(x);
+#endif // LIBC_TARGET_CPU_HAS_FPU_FLOAT
+
+ // Exact power of 2:
+ // pow(2^k, y) = 2^(k * y).
+ // This is exact if k * y is an integer.
+ uint32_t y_a = y_abs.uintval();
+ if (x_abs.is_normal() && xbits.get_mantissa() == 0 && y_a > Y_LOWER_BOUND &&
+ y_a < Y_UPPER_BOUND) {
+ if (!x_sign || is_integer(y)) {
+ Sign out_sign = (x_sign && is_odd_integer(y)) ? Sign::NEG : Sign::POS;
+ int e_x = xbits.get_exponent();
+ double ey = static_cast<double>(e_x) * static_cast<double>(y);
+ double ey_int = fputil::nearest_integer(ey);
+ if (ey_int == ey) {
+ if (ey > 127.0)
+ return set_overflow(out_sign);
+ if (ey < -149.0)
+ return set_underflow(out_sign);
+ int k = static_cast<int>(ey);
+ if (k >= -126)
+ return FPBits::create_value(out_sign, static_cast<uint32_t>(k + 127),
+ 0)
+ .get_val();
+ return FPBits::create_value(out_sign, 0, 1U << (k + 149)).get_val();
+ }
+ }
+ }
+
+#ifndef LIBC_MATH_HAS_SMALL_TABLES
+ if (x == 2.0f)
+ return math::exp2f(y);
+
+ if (x == 10.0f)
+ return math::exp10f(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)) {
+ float hi = static_cast<float>(e_x) * y;
+ if (hi >= 128.0f)
+ return set_overflow();
+ if (hi < -150.0f)
+ return set_underflow();
+ if (hi >= -126.0f)
+ return math::exp2f(hi);
+ }
+ }
+#endif // !LIBC_MATH_HAS_SMALL_TABLES
+
+ 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, ex, sign may be updated).
+template <typename SignType = uint64_t>
+LIBC_ALWAYS_INLINE cpp::optional<float>
+check_exceptional_cases(float &x, float &y, int &ex, SignType &sign) {
+ using FloatBits = fputil::FPBits<float>;
+ FloatBits xbits(x), ybits(y);
+ bool x_sign = xbits.sign() == Sign::NEG;
+ bool y_sign = ybits.sign() == Sign::NEG;
+
+ FloatBits x_abs = xbits.abs();
+ FloatBits y_abs = ybits.abs();
+
+ uint32_t x_a = x_abs.uintval();
+ uint32_t y_a = y_abs.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 FloatBits::quiet_nan().get_val();
+ }
+
+ // Extreme |y| >= Y_UPPER_BOUND: |y * log2(x)| >= 150 for all x != 1.
+ if (y_a >= Y_UPPER_BOUND) {
+ if (x_abs.is_nan())
+ return x;
+ if (ybits.is_nan())
+ return y;
+
+ if (x_a == FloatBits::one().uintval())
+ return 1.0f;
+
+ bool is_overflow = (x_a < FloatBits::one().uintval()) == y_sign;
+ if (y_abs.is_inf() || x_a == FloatBits::inf().uintval())
+ return is_overflow ? FloatBits::inf().get_val() : 0.0f;
+
+ return is_overflow ? set_overflow() : set_underflow();
+ }
+
+ // y is finite and non-zero.
+
+ if (x_a == FloatBits::inf().uintval()) {
+ // pow(+-Inf, y):
+ // - y < 0: returns +-0.0 (negative if x = -Inf and y is odd integer).
+ // - y > 0: returns +-Inf (negative if x = -Inf and y is odd integer).
+ bool out_is_neg = x_sign && is_odd_integer(y);
+ Sign out_sign = out_is_neg ? Sign::NEG : Sign::POS;
+ return y_sign ? FloatBits::zero(out_sign).get_val()
+ : FloatBits::inf(out_sign).get_val();
+ }
+
+ if (x_a > FloatBits::inf().uintval()) {
+ // x is NaN.
+ return x;
+ }
+
+ // Normalize denormal inputs by scaling with 2^64.
+ if (x_a < FloatBits::min_normal().uintval()) {
+ ex -= 64;
+ x *= 0x1.0p64f;
+ }
+
+ // Handle negative base x < 0:
+ // - If y is an integer: (-|x|)^y = (-1)^y * |x|^y.
+ // - If y is not an integer: (-|x|)^y is undefined in real numbers;
+ // raises FE_INVALID, sets errno to EDOM, and returns quiet NaN.
+ if (x_sign) {
+ if (is_integer(y)) {
+ x = -x;
+ if (is_odd_integer(y)) {
+ if constexpr (sizeof(SignType) == 8)
+ sign = 0x8000'0000'0000'0000ULL;
+ else
+ sign = 0x8000'0000U;
+ }
+ } else {
+ // pow( negative, non-integer ) = NaN
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FloatBits::quiet_nan().get_val();
+ }
+ }
+
+ if (y_a <= Y_LOWER_BOUND) {
+ volatile float one = 1.0f;
+ volatile float eps = ((x_a < FloatBits::one().uintval()) == !y_sign)
+ ? -0x1.0p-50f
+ : 0x1.0p-50f;
+ return one + eps;
+ }
+
+ return cpp::nullopt;
+}
+
+#ifndef LIBC_MATH_HAS_SMALL_TABLES
+// Check if x^y is an exact rounding boundary case (a 25-bit dyadic float,
+// which is either an exact 24-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.
+//
+// 1. x = 2^e: x^y = 2^(e * y) is in F_25 iff e * y is an integer.
+// 2. x = m * 2^e, y = n * 2^f with m, n odd integers:
+// x^y in F_25 implies:
+// - 0 <= y <= 15, n <= 15.
+// - f >= -3.
+// - 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^25.
+//
+// Returns true and sets exact_m and exact_exp if x^y is an exact boundary.
+LIBC_INLINE bool is_exact_rounding_boundary(float x, float y, uint32_t &exact_m,
+ int &exact_exp) {
+ using FPBits = fputil::FPBits<float>;
+
+ FPBits xbits(x);
+ int x_e = 0;
+ uint32_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) - 8;
+ x_mant = (x_mant << shift) & FPBits::FRACTION_MASK;
+ x_e = -126 - 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 * static_cast<double>(y);
+ if (fputil::nearest_integer(ey) == ey) {
+ exact_m = 1;
+ if (ey > 150.0)
+ exact_exp = 200;
+ else if (ey < -200.0)
+ exact_exp = -200;
+ else
+ exact_exp = static_cast<int>(ey);
+ return true;
+ }
+ return false;
+ }
+
+ // Case 2: x is not a power of 2.
+ if (y < 0.0f || y > 15.0f)
+ return false;
+
+ // Decompose y = n * 2^f with n odd.
+ FPBits ybits(y);
+ uint32_t y_mant = ybits.get_mantissa() | (1U << FPBits::FRACTION_LEN);
+ int y_exp = ybits.get_exponent() - static_cast<int>(FPBits::FRACTION_LEN);
+ int tz_y = cpp::countr_zero(y_mant);
+ uint32_t n = y_mant >> tz_y;
+ int f = y_exp + tz_y;
+
+ if (n > 15 || f < -3)
+ return false;
+
+ // Decompose x = m * 2^e with m odd.
+ uint32_t full_x_mant = x_mant | (1U << FPBits::FRACTION_LEN);
+ int tz_x = cpp::countr_zero(full_x_mant);
+ uint32_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, 2, 3}.
+ double ey = static_cast<double>(e) * static_cast<double>(y);
+ if (fputil::nearest_integer(ey) != ey)
+ return false;
+
+ // Check if m^(2^f) is an integer by taking repeated square roots -f times.
+ int count = -f;
+ uint32_t cur = m;
+ for (int i = 0; i < count; ++i) {
+ uint32_t s =
+ static_cast<uint32_t>(fputil::sqrt<double>(static_cast<double>(cur)));
+ if (s * s != cur)
+ return false;
+ cur = s;
+ }
+
+ // Compute res = cur^n and check if res <= 2^25.
+ uint32_t res = 1;
+ for (uint32_t i = 0; i < n; ++i) {
+ if (res > (1U << 25) / 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);
+ uint32_t res = 1;
+ for (int i = 0; i < y_int; ++i) {
+ if (res > (1U << 25) / m)
+ return false;
+ res *= m;
+ }
+ exact_m = res;
+ exact_exp = e * y_int;
+ return true;
+ }
+}
+#endif // !LIBC_MATH_HAS_SMALL_TABLES
+
+} // namespace powf_internal
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_UTILS_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 1fd0772aa50314..f1743784478bc4 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -3075,6 +3075,8 @@ add_fp_unittest(
powf_test.cpp
DEPENDS
libc.src.math.powf
+ libc.src.__support.math.powf_double_eval
+ libc.src.__support.math.powf_float_eval
libc.src.__support.FPUtil.fp_bits
)
diff --git a/libc/test/src/math/powf_test.cpp b/libc/test/src/math/powf_test.cpp
index fe8ef4fa85117c..b97a34e03292f2 100644
--- a/libc/test/src/math/powf_test.cpp
+++ b/libc/test/src/math/powf_test.cpp
@@ -1,15 +1,22 @@
-//===-- Unittests for powf ------------------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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 powf.
+///
+//===----------------------------------------------------------------------===//
#include "hdr/math_macros.h"
#include "hdr/stdint_proxy.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/macros/optimization.h"
+#include "src/__support/math/powf_double_eval.h"
+#include "src/__support/math/powf_float_eval.h"
#include "src/math/powf.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
@@ -21,116 +28,150 @@
#define TOLERANCE 0
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-using LlvmLibcPowfTest = LIBC_NAMESPACE::testing::FPTest<float>;
using LIBC_NAMESPACE::testing::tlog;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
-TEST_F(LlvmLibcPowfTest, TrickyInputs) {
- constexpr int N = 13;
- constexpr mpfr::BinaryInput<float> INPUTS[N] = {
- {0x1.290bbp-124f, 0x1.1e6d92p-25f},
- {0x1.2e9fb6p+5f, -0x1.1b82b6p-18f},
- {0x1.6877f6p+60f, -0x1.75f1c6p-4f},
- {0x1.0936acp-63f, -0x1.55200ep-15f},
- {0x1.d6d72ap+43f, -0x1.749ccap-5f},
- {0x1.4afb2ap-40f, 0x1.063198p+0f},
- {0x1.0124dep+0f, -0x1.fdb016p+9f},
- {0x1.1058p+0f, 0x1.ap+64f},
- {0x1.1058p+0f, -0x1.ap+64f},
- {0x1.1058p+0f, 0x1.ap+64f},
- {0x1.fa32d4p-1f, 0x1.67a62ep+12f},
- {-0x1.8p-49, 0x1.8p+1},
- {0x1.8p-48, 0x1.8p+1},
- };
-
- for (int i = 0; i < N; ++i) {
- float x = INPUTS[i].x;
- float y = INPUTS[i].y;
- EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, INPUTS[i],
- LIBC_NAMESPACE::powf(x, y), TOLERANCE + 0.5);
+class PowfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
+public:
+ void test_tricky_inputs(float (*func)(float, float), double ulp_tolerance,
+ bool all_rounding) {
+ constexpr int N = 13;
+ constexpr mpfr::BinaryInput<float> INPUTS[N] = {
+ {0x1.290bbp-124f, 0x1.1e6d92p-25f},
+ {0x1.2e9fb6p+5f, -0x1.1b82b6p-18f},
+ {0x1.6877f6p+60f, -0x1.75f1c6p-4f},
+ {0x1.0936acp-63f, -0x1.55200ep-15f},
+ {0x1.d6d72ap+43f, -0x1.749ccap-5f},
+ {0x1.4afb2ap-40f, 0x1.063198p+0f},
+ {0x1.0124dep+0f, -0x1.fdb016p+9f},
+ {0x1.1058p+0f, 0x1.ap+64f},
+ {0x1.1058p+0f, -0x1.ap+64f},
+ {0x1.1058p+0f, 0x1.ap+64f},
+ {0x1.fa32d4p-1f, 0x1.67a62ep+12f},
+ {-0x1.8p-49, 0x1.8p+1},
+ {0x1.8p-48, 0x1.8p+1},
+ };
+
+ for (int i = 0; i < N; ++i) {
+ if (all_rounding) {
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, INPUTS[i],
+ func(INPUTS[i].x, INPUTS[i].y),
+ ulp_tolerance);
+ } else {
+ EXPECT_MPFR_MATCH(mpfr::Operation::Pow, INPUTS[i],
+ func(INPUTS[i].x, INPUTS[i].y), ulp_tolerance);
+ }
+ }
}
-}
-
-TEST_F(LlvmLibcPowfTest, InFloatRange) {
- constexpr uint32_t X_COUNT = 1'23;
- constexpr uint32_t X_START = FPBits(0.25f).uintval();
- constexpr uint32_t X_STOP = FPBits(4.0f).uintval();
- constexpr uint32_t X_STEP = (X_STOP - X_START) / X_COUNT;
-
- constexpr uint32_t Y_COUNT = 1'37;
- constexpr uint32_t Y_START = FPBits(0.25f).uintval();
- constexpr uint32_t Y_STOP = FPBits(4.0f).uintval();
- constexpr uint32_t Y_STEP = (Y_STOP - Y_START) / Y_COUNT;
-
- auto test = [&](mpfr::RoundingMode rounding_mode) {
- mpfr::ForceRoundingMode __r(rounding_mode);
- if (!__r.success)
- return;
-
- uint64_t fails = 0;
- uint64_t count = 0;
- uint64_t cc = 0;
- float mx, my, mr = 0.0;
- double tol = 0.5;
-
- for (uint32_t i = 0, v = X_START; i <= X_COUNT; ++i, v += X_STEP) {
- float x = FPBits(v).get_val();
- if (FPBits(v).is_nan() || FPBits(v).is_inf() || x < 0.0)
- continue;
-
- for (uint32_t j = 0, w = Y_START; j <= Y_COUNT; ++j, w += Y_STEP) {
- float y = FPBits(w).get_val();
- if (FPBits(w).is_nan() || FPBits(w).is_inf())
- continue;
- libc_errno = 0;
- float result = LIBC_NAMESPACE::powf(x, y);
- ++cc;
- if (FPBits(result).is_nan() || FPBits(result).is_inf())
+ void test_in_range(float (*func)(float, float), double ulp_tolerance,
+ bool all_rounding) {
+ constexpr uint32_t X_COUNT = 1'23;
+ constexpr uint32_t X_START = FPBits(0.25f).uintval();
+ constexpr uint32_t X_STOP = FPBits(4.0f).uintval();
+ constexpr uint32_t X_STEP = (X_STOP - X_START) / X_COUNT;
+
+ constexpr uint32_t Y_COUNT = 1'37;
+ constexpr uint32_t Y_START = FPBits(0.25f).uintval();
+ constexpr uint32_t Y_STOP = FPBits(4.0f).uintval();
+ constexpr uint32_t Y_STEP = (Y_STOP - Y_START) / Y_COUNT;
+
+ auto test = [&](mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode __r(rounding_mode);
+ if (!__r.success)
+ return;
+
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ float mx = 0.0f, my = 0.0f, mr = 0.0f;
+ double tol = ulp_tolerance;
+
+ for (uint32_t i = 0, v = X_START; i <= X_COUNT; ++i, v += X_STEP) {
+ float x = FPBits(v).get_val();
+ if (FPBits(v).is_nan() || FPBits(v).is_inf() || x < 0.0f)
continue;
- ++count;
- mpfr::BinaryInput<float> inputs{x, y};
-
- if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Pow, inputs,
- result, TOLERANCE + 0.5,
- rounding_mode)) {
- ++fails;
- while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
- mpfr::Operation::Pow, inputs, result, tol, rounding_mode)) {
- mx = x;
- my = y;
- mr = result;
-
- if (tol > 1000.0)
- break;
-
- tol *= 2.0;
+ for (uint32_t j = 0, w = Y_START; j <= Y_COUNT; ++j, w += Y_STEP) {
+ float y = FPBits(w).get_val();
+ if (FPBits(w).is_nan() || FPBits(w).is_inf())
+ continue;
+
+ libc_errno = 0;
+ float result = func(x, y);
+ ++cc;
+ if (FPBits(result).is_nan() || FPBits(result).is_inf())
+ continue;
+
+ ++count;
+ mpfr::BinaryInput<float> inputs{x, y};
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Pow, inputs,
+ result, ulp_tolerance,
+ rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Pow, inputs, result, tol, rounding_mode)) {
+ mx = x;
+ my = y;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
}
}
}
- }
- if (fails || (count < cc)) {
- tlog << " Powf failed: " << fails << "/" << count << "/" << cc
- << " tests.\n"
- << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
- }
- if (fails) {
- mpfr::BinaryInput<float> inputs{mx, my};
- EXPECT_MPFR_MATCH(mpfr::Operation::Pow, inputs, mr, 0.5, rounding_mode);
- }
- };
+ if (fails || (count < cc)) {
+ tlog << " Powf failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n"
+ << " Max ULPs is at most: " << static_cast<uint64_t>(tol)
+ << ".\n";
+ }
+ if (fails) {
+ mpfr::BinaryInput<float> inputs{mx, my};
+ EXPECT_MPFR_MATCH(mpfr::Operation::Pow, inputs, mr, ulp_tolerance,
+ rounding_mode);
+ }
+ };
- tlog << " Test Rounding To Nearest...\n";
- test(mpfr::RoundingMode::Nearest);
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
- tlog << " Test Rounding Downward...\n";
- test(mpfr::RoundingMode::Downward);
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (all_rounding) {
+ tlog << " Test Rounding Downward...\n";
+ test(mpfr::RoundingMode::Downward);
- tlog << " Test Rounding Upward...\n";
- test(mpfr::RoundingMode::Upward);
+ tlog << " Test Rounding Upward...\n";
+ test(mpfr::RoundingMode::Upward);
- tlog << " Test Rounding Toward Zero...\n";
- test(mpfr::RoundingMode::TowardZero);
-}
+ tlog << " Test Rounding Toward Zero...\n";
+ test(mpfr::RoundingMode::TowardZero);
+ }
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ }
+};
+
+#define LIST_POWF_TESTS(suffix, func, tricky_ulp, range_ulp, all_rounding) \
+ using LlvmLibcPowfTest##suffix = PowfTest; \
+ TEST_F(LlvmLibcPowfTest##suffix, TrickyInputs) { \
+ test_tricky_inputs(&func, tricky_ulp, all_rounding); \
+ } \
+ TEST_F(LlvmLibcPowfTest##suffix, InFloatRange) { \
+ test_in_range(&func, range_ulp, all_rounding); \
+ } \
+ static_assert(true, "Require semicolon.")
+
+LIST_POWF_TESTS(Default, LIBC_NAMESPACE::powf,
+ /*tricky_ulp=*/TOLERANCE + 0.5, /*range_ulp=*/TOLERANCE + 0.5,
+ /*all_rounding=*/true);
+LIST_POWF_TESTS(DoubleEval, LIBC_NAMESPACE::math::double_eval::powf,
+ /*tricky_ulp=*/TOLERANCE + 0.5, /*range_ulp=*/TOLERANCE + 0.5,
+ /*all_rounding=*/true);
+LIST_POWF_TESTS(FloatEval, LIBC_NAMESPACE::math::float_eval::powf,
+ /*tricky_ulp=*/1.0, /*range_ulp=*/1.5,
+ /*all_rounding=*/false);
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 0bfe73a2e643a7..6037e373fad14f 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -5363,6 +5363,8 @@ add_fp_unittest(
powf_test.cpp
DEPENDS
libc.src.math.powf
+ libc.src.__support.math.powf_double_eval
+ libc.src.__support.math.powf_float_eval
libc.src.__support.FPUtil.fp_bits
)
diff --git a/libc/test/src/math/smoke/powf_test.cpp b/libc/test/src/math/smoke/powf_test.cpp
index 65c4da35b6db07..c20c02ddd89578 100644
--- a/libc/test/src/math/smoke/powf_test.cpp
+++ b/libc/test/src/math/smoke/powf_test.cpp
@@ -1,271 +1,301 @@
-//===-- Unittests for powf ------------------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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 powf.
+///
+//===----------------------------------------------------------------------===//
+#include "hdr/fenv_macros.h"
#include "hdr/math_macros.h"
#include "hdr/stdint_proxy.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/math/powf_double_eval.h"
+#include "src/__support/math/powf_float_eval.h"
#include "src/math/powf.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
-using LlvmLibcPowfTest = LIBC_NAMESPACE::testing::FPTest<float>;
using LIBC_NAMESPACE::fputil::testing::ForceRoundingMode;
using LIBC_NAMESPACE::fputil::testing::RoundingMode;
-TEST_F(LlvmLibcPowfTest, SpecialNumbers) {
- constexpr float neg_odd_integer = -3.0f;
- constexpr float neg_even_integer = -6.0f;
- constexpr float neg_non_integer = -1.1f;
- constexpr float pos_odd_integer = 5.0f;
- constexpr float pos_even_integer = 8.0f;
- constexpr float pos_non_integer = 1.3f;
- constexpr float one_half = 0.5f;
+class PowfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
+public:
+ void test_special_numbers(float (*func)(float, float), int tolerance = 0) {
+ constexpr float neg_odd_integer = -3.0f;
+ constexpr float neg_even_integer = -6.0f;
+ constexpr float neg_non_integer = -1.1f;
+ constexpr float pos_odd_integer = 5.0f;
+ constexpr float pos_even_integer = 8.0f;
+ constexpr float pos_non_integer = 1.3f;
+ constexpr float one_half = 0.5f;
- for (int i = 0; i < N_ROUNDING_MODES; ++i) {
- ForceRoundingMode __r(ROUNDING_MODES[i]);
- if (!__r.success)
- continue;
+ for (int i = 0; i < N_ROUNDING_MODES; ++i) {
+ ForceRoundingMode __r(ROUNDING_MODES[i]);
+ if (!__r.success)
+ continue;
- // pow( sNaN, exponent)
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- aNaN, LIBC_NAMESPACE::powf(sNaN, neg_odd_integer), FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- aNaN, LIBC_NAMESPACE::powf(sNaN, neg_even_integer), FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- aNaN, LIBC_NAMESPACE::powf(sNaN, pos_odd_integer), FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- aNaN, LIBC_NAMESPACE::powf(sNaN, pos_even_integer), FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, one_half),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, zero),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, neg_zero),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, inf),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, neg_inf),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(sNaN, aNaN),
- FE_INVALID);
+ // pow( sNaN, exponent)
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, sNaN), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, neg_odd_integer),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, neg_even_integer),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, pos_odd_integer),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, pos_even_integer),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, one_half), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, zero), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, neg_zero), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, inf), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, neg_inf), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN, aNaN), FE_INVALID);
- // pow( 0.0f, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(zero, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- inf, LIBC_NAMESPACE::powf(zero, neg_odd_integer), FE_DIVBYZERO);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- inf, LIBC_NAMESPACE::powf(zero, neg_even_integer), FE_DIVBYZERO);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- inf, LIBC_NAMESPACE::powf(zero, neg_non_integer), FE_DIVBYZERO);
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(zero, pos_odd_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(zero, pos_even_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(zero, pos_non_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(zero, one_half));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(zero, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(zero, neg_zero));
- EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::powf(zero, inf));
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::powf(zero, neg_inf),
- FE_DIVBYZERO);
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(zero, aNaN));
+ // pow( 0.0f, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(zero, sNaN), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(zero, neg_odd_integer),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(zero, neg_even_integer),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(zero, neg_non_integer),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ(zero, func(zero, pos_odd_integer));
+ EXPECT_FP_EQ(zero, func(zero, pos_even_integer));
+ EXPECT_FP_EQ(zero, func(zero, pos_non_integer));
+ EXPECT_FP_EQ(zero, func(zero, one_half));
+ EXPECT_FP_EQ(1.0f, func(zero, zero));
+ EXPECT_FP_EQ(1.0f, func(zero, neg_zero));
+ EXPECT_FP_EQ(0.0f, func(zero, inf));
+ EXPECT_FP_EQ(inf, func(zero, neg_inf));
+ EXPECT_FP_IS_NAN(func(zero, aNaN));
- // pow( -0.0f, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(neg_zero, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- neg_inf, LIBC_NAMESPACE::powf(neg_zero, neg_odd_integer), FE_DIVBYZERO);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- inf, LIBC_NAMESPACE::powf(neg_zero, neg_even_integer), FE_DIVBYZERO);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- inf, LIBC_NAMESPACE::powf(neg_zero, neg_non_integer), FE_DIVBYZERO);
- EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::powf(neg_zero, pos_odd_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(neg_zero, pos_even_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(neg_zero, pos_non_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(neg_zero, one_half));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(neg_zero, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(neg_zero, neg_zero));
- EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::powf(neg_zero, inf));
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::powf(neg_zero, neg_inf),
- FE_DIVBYZERO);
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(neg_zero, aNaN));
+ // pow( -0.0f, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(neg_zero, sNaN), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, func(neg_zero, neg_odd_integer),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(neg_zero, neg_even_integer),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(neg_zero, neg_non_integer),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ(neg_zero, func(neg_zero, pos_odd_integer));
+ EXPECT_FP_EQ(zero, func(neg_zero, pos_even_integer));
+ EXPECT_FP_EQ(zero, func(neg_zero, pos_non_integer));
+ EXPECT_FP_EQ(zero, func(neg_zero, one_half));
+ EXPECT_FP_EQ(1.0f, func(neg_zero, zero));
+ EXPECT_FP_EQ(1.0f, func(neg_zero, neg_zero));
+ EXPECT_FP_EQ(0.0f, func(neg_zero, inf));
+ EXPECT_FP_EQ(inf, func(neg_zero, neg_inf));
+ EXPECT_FP_IS_NAN(func(neg_zero, aNaN));
- // pow( 1.0f, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(1.0f, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, neg_zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, 1.0f));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, -1.0f));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, neg_odd_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, neg_even_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, neg_non_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, pos_odd_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, pos_even_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, pos_non_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, inf));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, neg_inf));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(1.0f, aNaN));
+ // pow( 1.0f, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(1.0f, sNaN), FE_INVALID);
+ EXPECT_FP_EQ(1.0f, func(1.0f, neg_odd_integer));
+ EXPECT_FP_EQ(1.0f, func(1.0f, neg_even_integer));
+ EXPECT_FP_EQ(1.0f, func(1.0f, neg_non_integer));
+ EXPECT_FP_EQ(1.0f, func(1.0f, pos_odd_integer));
+ EXPECT_FP_EQ(1.0f, func(1.0f, pos_even_integer));
+ EXPECT_FP_EQ(1.0f, func(1.0f, pos_non_integer));
+ EXPECT_FP_EQ(1.0f, func(1.0f, one_half));
+ EXPECT_FP_EQ(1.0f, func(1.0f, zero));
+ EXPECT_FP_EQ(1.0f, func(1.0f, neg_zero));
+ EXPECT_FP_EQ(1.0f, func(1.0f, inf));
+ EXPECT_FP_EQ(1.0f, func(1.0f, neg_inf));
+ EXPECT_FP_EQ(1.0f, func(1.0f, aNaN));
- // pow( -1.0f, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(-1.0f, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(-1.0f, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(-1.0f, neg_zero));
- EXPECT_FP_EQ(-1.0f, LIBC_NAMESPACE::powf(-1.0f, 1.0f));
- EXPECT_FP_EQ(-1.0f, LIBC_NAMESPACE::powf(-1.0f, -1.0f));
- EXPECT_FP_EQ(-1.0f, LIBC_NAMESPACE::powf(-1.0f, neg_odd_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(-1.0f, neg_even_integer));
- EXPECT_FP_IS_NAN_WITH_EXCEPTION(
- LIBC_NAMESPACE::powf(-1.0f, neg_non_integer), FE_INVALID);
- EXPECT_FP_EQ(-1.0f, LIBC_NAMESPACE::powf(-1.0f, pos_odd_integer));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(-1.0f, pos_even_integer));
- EXPECT_FP_IS_NAN_WITH_EXCEPTION(
- LIBC_NAMESPACE::powf(-1.0f, pos_non_integer), FE_INVALID);
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(-1.0f, inf));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(-1.0f, neg_inf));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(-1.0f, aNaN));
+ // pow( -1.0f, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(-1.0f, sNaN), FE_INVALID);
+ EXPECT_FP_EQ(-1.0f, func(-1.0f, neg_odd_integer));
+ EXPECT_FP_EQ(1.0f, func(-1.0f, neg_even_integer));
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(-1.0f, neg_non_integer),
+ FE_INVALID);
+ EXPECT_FP_EQ(-1.0f, func(-1.0f, pos_odd_integer));
+ EXPECT_FP_EQ(1.0f, func(-1.0f, pos_even_integer));
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(-1.0f, pos_non_integer),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(-1.0f, one_half), FE_INVALID);
+ EXPECT_FP_EQ(1.0f, func(-1.0f, zero));
+ EXPECT_FP_EQ(1.0f, func(-1.0f, neg_zero));
+ EXPECT_FP_EQ(1.0f, func(-1.0f, inf));
+ EXPECT_FP_EQ(1.0f, func(-1.0f, neg_inf));
+ EXPECT_FP_IS_NAN(func(-1.0f, aNaN));
- // pow( inf, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(inf, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(inf, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(inf, neg_zero));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(inf, 1.0f));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(inf, -1.0f));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(inf, neg_odd_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(inf, neg_even_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(inf, neg_non_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(inf, pos_odd_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(inf, pos_even_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(inf, pos_non_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(inf, one_half));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(inf, inf));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(inf, neg_inf));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(inf, aNaN));
+ // pow( inf, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(inf, sNaN), FE_INVALID);
+ EXPECT_FP_EQ(0.0f, func(inf, neg_odd_integer));
+ EXPECT_FP_EQ(0.0f, func(inf, neg_even_integer));
+ EXPECT_FP_EQ(0.0f, func(inf, neg_non_integer));
+ EXPECT_FP_EQ(inf, func(inf, pos_odd_integer));
+ EXPECT_FP_EQ(inf, func(inf, pos_even_integer));
+ EXPECT_FP_EQ(inf, func(inf, pos_non_integer));
+ EXPECT_FP_EQ(inf, func(inf, one_half));
+ EXPECT_FP_EQ(1.0f, func(inf, zero));
+ EXPECT_FP_EQ(1.0f, func(inf, neg_zero));
+ EXPECT_FP_EQ(inf, func(inf, inf));
+ EXPECT_FP_EQ(0.0f, func(inf, neg_inf));
+ EXPECT_FP_IS_NAN(func(inf, aNaN));
- // pow( -inf, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(neg_inf, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(neg_inf, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(neg_inf, neg_zero));
- EXPECT_FP_EQ(neg_inf, LIBC_NAMESPACE::powf(neg_inf, 1.0f));
- EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::powf(neg_inf, -1.0f));
- EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::powf(neg_inf, neg_odd_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(neg_inf, neg_even_integer));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(neg_inf, neg_non_integer));
- EXPECT_FP_EQ(neg_inf, LIBC_NAMESPACE::powf(neg_inf, pos_odd_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(neg_inf, pos_even_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(neg_inf, pos_non_integer));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(neg_inf, one_half));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(neg_inf, inf));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(neg_inf, neg_inf));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(neg_inf, aNaN));
+ // pow( -inf, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(neg_inf, sNaN), FE_INVALID);
+ EXPECT_FP_EQ(-0.0f, func(neg_inf, neg_odd_integer));
+ EXPECT_FP_EQ(0.0f, func(neg_inf, neg_even_integer));
+ EXPECT_FP_EQ(0.0f, func(neg_inf, neg_non_integer));
+ EXPECT_FP_EQ(neg_inf, func(neg_inf, pos_odd_integer));
+ EXPECT_FP_EQ(inf, func(neg_inf, pos_even_integer));
+ EXPECT_FP_EQ(inf, func(neg_inf, pos_non_integer));
+ EXPECT_FP_EQ(inf, func(neg_inf, one_half));
+ EXPECT_FP_EQ(1.0f, func(neg_inf, zero));
+ EXPECT_FP_EQ(1.0f, func(neg_inf, neg_zero));
+ EXPECT_FP_EQ(inf, func(neg_inf, inf));
+ EXPECT_FP_EQ(0.0f, func(neg_inf, neg_inf));
+ EXPECT_FP_IS_NAN(func(neg_inf, aNaN));
- // pow ( aNaN, exponent )
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(aNaN, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(aNaN, zero));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(aNaN, neg_zero));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, 1.0f));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, -1.0f));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, neg_odd_integer));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, neg_even_integer));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, neg_non_integer));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, pos_odd_integer));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, pos_even_integer));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, pos_non_integer));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, inf));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, neg_inf));
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(aNaN, aNaN));
+ // pow( aNaN, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(aNaN, sNaN), FE_INVALID);
+ EXPECT_FP_IS_NAN(func(aNaN, neg_odd_integer));
+ EXPECT_FP_IS_NAN(func(aNaN, neg_even_integer));
+ EXPECT_FP_IS_NAN(func(aNaN, neg_non_integer));
+ EXPECT_FP_IS_NAN(func(aNaN, pos_odd_integer));
+ EXPECT_FP_IS_NAN(func(aNaN, pos_even_integer));
+ EXPECT_FP_IS_NAN(func(aNaN, pos_non_integer));
+ EXPECT_FP_IS_NAN(func(aNaN, one_half));
+ EXPECT_FP_EQ(1.0f, func(aNaN, zero));
+ EXPECT_FP_EQ(1.0f, func(aNaN, neg_zero));
+ EXPECT_FP_IS_NAN(func(aNaN, inf));
+ EXPECT_FP_IS_NAN(func(aNaN, neg_inf));
+ EXPECT_FP_IS_NAN(func(aNaN, aNaN));
- // pow ( base, inf )
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(0.1f, inf));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(-0.1f, inf));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(1.1f, inf));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(-1.1f, inf));
+ // Exact powers of 2:
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(2.0f, sNaN), FE_INVALID);
+ EXPECT_FP_EQ(0x1.0p15f, func(2.0f, 15.0f));
+ EXPECT_FP_EQ(0x1.0p126f, func(2.0f, 126.0f));
+ EXPECT_FP_EQ(0x1.0p-45f, func(2.0f, -45.0f));
+ EXPECT_FP_EQ(0x1.0p-126f, func(2.0f, -126.0f));
+ EXPECT_FP_EQ(0x1.0p-149f, func(2.0f, -149.0f));
- // pow ( base, -inf )
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(0.1f, neg_inf));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf(-0.1f, neg_inf));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(1.1f, neg_inf));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf(-1.1f, neg_inf));
+ // Powers of 10:
+ EXPECT_FP_EQ(1.0f, func(10.0f, 0.0f));
+ EXPECT_FP_EQ(10.0f, func(10.0f, 1.0f));
+ EXPECT_FP_EQ(100.0f, func(10.0f, 2.0f));
+ if (tolerance == 0) {
+ EXPECT_FP_EQ(1000.0f, func(10.0f, 3.0f));
+ EXPECT_FP_EQ(10000.0f, func(10.0f, 4.0f));
+ EXPECT_FP_EQ(100000.0f, func(10.0f, 5.0f));
+ EXPECT_FP_EQ(1000000.0f, func(10.0f, 6.0f));
+ EXPECT_FP_EQ(10000000.0f, func(10.0f, 7.0f));
+ EXPECT_FP_EQ(100000000.0f, func(10.0f, 8.0f));
+ EXPECT_FP_EQ(1000000000.0f, func(10.0f, 9.0f));
+ EXPECT_FP_EQ(10000000000.0f, func(10.0f, 10.0f));
+ } else {
+ auto check_close = [tolerance](float act, float exp) {
+ uint32_t a = FPBits(act).uintval();
+ uint32_t e = FPBits(exp).uintval();
+ EXPECT_LE(a >= e ? a - e : e - a, static_cast<uint32_t>(tolerance));
+ };
+ check_close(func(10.0f, 3.0f), 1000.0f);
+ check_close(func(10.0f, 4.0f), 10000.0f);
+ check_close(func(10.0f, 5.0f), 100000.0f);
+ check_close(func(10.0f, 6.0f), 1000000.0f);
+ check_close(func(10.0f, 7.0f), 10000000.0f);
+ check_close(func(10.0f, 8.0f), 100000000.0f);
+ check_close(func(10.0f, 9.0f), 1000000000.0f);
+ check_close(func(10.0f, 10.0f), 10000000000.0f);
+ }
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(10.0f, sNaN), FE_INVALID);
- // Exact powers of 2:
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(2.0f, sNaN),
- FE_INVALID);
- EXPECT_FP_EQ(0x1.0p15f, LIBC_NAMESPACE::powf(2.0f, 15.0f));
- EXPECT_FP_EQ(0x1.0p126f, LIBC_NAMESPACE::powf(2.0f, 126.0f));
- EXPECT_FP_EQ(0x1.0p-45f, LIBC_NAMESPACE::powf(2.0f, -45.0f));
- EXPECT_FP_EQ(0x1.0p-126f, LIBC_NAMESPACE::powf(2.0f, -126.0f));
- EXPECT_FP_EQ(0x1.0p-149f, LIBC_NAMESPACE::powf(2.0f, -149.0f));
+ // Overflow / Underflow:
+ if (ROUNDING_MODES[i] != RoundingMode::Downward &&
+ ROUNDING_MODES[i] != RoundingMode::TowardZero) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(3.1f, 201.0f), FE_OVERFLOW);
+ }
+ if (ROUNDING_MODES[i] != RoundingMode::Upward) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(0.0f, func(3.1f, -201.0f), FE_UNDERFLOW);
+ }
+ }
- // Exact powers of 10:
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(10.0f, 0.0f));
- EXPECT_FP_EQ(10.0f, LIBC_NAMESPACE::powf(10.0f, 1.0f));
- EXPECT_FP_EQ(100.0f, LIBC_NAMESPACE::powf(10.0f, 2.0f));
- EXPECT_FP_EQ(1000.0f, LIBC_NAMESPACE::powf(10.0f, 3.0f));
- EXPECT_FP_EQ(10000.0f, LIBC_NAMESPACE::powf(10.0f, 4.0f));
- EXPECT_FP_EQ(100000.0f, LIBC_NAMESPACE::powf(10.0f, 5.0f));
- EXPECT_FP_EQ(1000000.0f, LIBC_NAMESPACE::powf(10.0f, 6.0f));
- EXPECT_FP_EQ(10000000.0f, LIBC_NAMESPACE::powf(10.0f, 7.0f));
- EXPECT_FP_EQ(100000000.0f, LIBC_NAMESPACE::powf(10.0f, 8.0f));
- EXPECT_FP_EQ(1000000000.0f, LIBC_NAMESPACE::powf(10.0f, 9.0f));
- EXPECT_FP_EQ(10000000000.0f, LIBC_NAMESPACE::powf(10.0f, 10.0f));
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf(10.0f, sNaN),
- FE_INVALID);
+ EXPECT_FP_EQ(-0.0f, func(-0.015625f, 25.0f));
+ EXPECT_FP_EQ(0.0f, func(-0.015625f, 26.0f));
+ }
- // Overflow / Underflow:
- if (ROUNDING_MODES[i] != RoundingMode::Downward &&
- ROUNDING_MODES[i] != RoundingMode::TowardZero) {
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::powf(3.1f, 201.0f),
- FE_OVERFLOW);
- }
- if (ROUNDING_MODES[i] != RoundingMode::Upward) {
- EXPECT_FP_EQ_WITH_EXCEPTION(0.0f, LIBC_NAMESPACE::powf(3.1f, -201.0f),
- FE_UNDERFLOW);
+ void test_subnormal_base(float (*func)(float, float), int tolerance = 0) {
+ EXPECT_FP_EQ(0x1.0p-32f, func(0x1.0p-128f, 0.25f));
+ EXPECT_FP_EQ(0x1.0p96f, func(0x1.0p-128f, -0.75f));
+ if (tolerance == 0) {
+ EXPECT_FP_EQ(0x1.90a962p-33f, func(0x1.8p-130f, 0.25f));
+ EXPECT_FP_EQ(0x1.47238cp+32f, func(0x1.8p-130f, -0.25f));
+ } else {
+ uint32_t act1 = FPBits(func(0x1.8p-130f, 0.25f)).uintval();
+ uint32_t exp1 = FPBits(0x1.90a962p-33f).uintval();
+ EXPECT_LE(act1 >= exp1 ? act1 - exp1 : exp1 - act1,
+ static_cast<uint32_t>(tolerance));
+ uint32_t act2 = FPBits(func(0x1.8p-130f, -0.25f)).uintval();
+ uint32_t exp2 = FPBits(0x1.47238cp+32f).uintval();
+ EXPECT_LE(act2 >= exp2 ? act2 - exp2 : exp2 - act2,
+ static_cast<uint32_t>(tolerance));
}
}
- EXPECT_FP_EQ(-0.0f, LIBC_NAMESPACE::powf(-0.015625f, 25.0f));
- EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::powf(-0.015625f, 26.0f));
-}
-
-TEST_F(LlvmLibcPowfTest, SubnormalBase) {
- EXPECT_FP_EQ(0x1.0p-32f, LIBC_NAMESPACE::powf(0x1.0p-128f, 0.25f));
- EXPECT_FP_EQ(0x1.0p96f, LIBC_NAMESPACE::powf(0x1.0p-128f, -0.75f));
- EXPECT_FP_EQ(0x1.90a962p-33f, LIBC_NAMESPACE::powf(0x1.8p-130f, 0.25f));
- EXPECT_FP_EQ(0x1.47238cp+32f, LIBC_NAMESPACE::powf(0x1.8p-130f, -0.25f));
-}
-
#ifdef LIBC_TEST_FTZ_DAZ
+ void test_ftz(float (*func)(float, float)) {
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ);
+ volatile float x = -min_denormal;
+ volatile float y = 0.5f;
+ EXPECT_FP_IS_NAN(func(x, y));
+ volatile float two = 2.0f;
+ volatile float d = min_denormal;
+ EXPECT_FP_EQ(1.0f, func(two, d));
+ }
-using namespace LIBC_NAMESPACE::testing;
-
-TEST_F(LlvmLibcPowfTest, FTZMode) {
- ModifyMXCSR mxcsr(FTZ);
-
- EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf(-min_denormal, 0.5f));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(2.0f, min_denormal));
-}
-
-TEST_F(LlvmLibcPowfTest, DAZMode) {
- ModifyMXCSR mxcsr(DAZ);
+ void test_daz(float (*func)(float, float)) {
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::DAZ);
+ volatile float x = -min_denormal;
+ volatile float y = 0.5f;
+ EXPECT_FP_EQ(0.0f, func(x, y));
+ volatile float two = 2.0f;
+ volatile float d = min_denormal;
+ EXPECT_FP_EQ(1.0f, func(two, d));
+ }
- EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::powf(-min_denormal, 0.5f));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(2.0f, min_denormal));
-}
+ void test_ftzdaz(float (*func)(float, float)) {
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ |
+ LIBC_NAMESPACE::testing::DAZ);
+ volatile float x = -min_denormal;
+ volatile float y = 0.5f;
+ EXPECT_FP_EQ(0.0f, func(x, y));
+ volatile float two = 2.0f;
+ volatile float d = min_denormal;
+ EXPECT_FP_EQ(1.0f, func(two, d));
+ }
+#endif // LIBC_TEST_FTZ_DAZ
+};
-TEST_F(LlvmLibcPowfTest, FTZDAZMode) {
- ModifyMXCSR mxcsr(FTZ | DAZ);
+#ifdef LIBC_TEST_FTZ_DAZ
+#define LIST_POWF_FTZ_DAZ_TESTS(suffix, func) \
+ TEST_F(LlvmLibcPowfTest##suffix, FTZMode) { test_ftz(&func); } \
+ TEST_F(LlvmLibcPowfTest##suffix, DAZMode) { test_daz(&func); } \
+ TEST_F(LlvmLibcPowfTest##suffix, FTZDAZMode) { test_ftzdaz(&func); }
+#else
+#define LIST_POWF_FTZ_DAZ_TESTS(suffix, func)
+#endif // LIBC_TEST_FTZ_DAZ
- EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::powf(-min_denormal, 0.5f));
- EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::powf(2.0f, min_denormal));
-}
+#define LIST_POWF_TESTS(suffix, func, tolerance) \
+ using LlvmLibcPowfTest##suffix = PowfTest; \
+ TEST_F(LlvmLibcPowfTest##suffix, SpecialNumbers) { \
+ test_special_numbers(&func, tolerance); \
+ } \
+ TEST_F(LlvmLibcPowfTest##suffix, SubnormalBase) { \
+ test_subnormal_base(&func, tolerance); \
+ } \
+ LIST_POWF_FTZ_DAZ_TESTS(suffix, func) \
+ static_assert(true, "Require semicolon.")
-#endif
+LIST_POWF_TESTS(Default, LIBC_NAMESPACE::powf, /*tolerance=*/0);
+LIST_POWF_TESTS(DoubleEval, LIBC_NAMESPACE::math::double_eval::powf,
+ /*tolerance=*/0);
+LIST_POWF_TESTS(FloatEval, LIBC_NAMESPACE::math::float_eval::powf,
+ /*tolerance=*/1);
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 33fc35455fef30..7df5b115f2a5ab 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -10185,29 +10185,84 @@ libc_support_library(
)
libc_support_library(
- name = "__support_math_powf",
- hdrs = [
- "src/__support/math/powf.h",
- "src/__support/math/powf_small_tables.h",
- ],
+ name = "__support_math_powf_utils",
+ hdrs = ["src/__support/math/powf_utils.h"],
deps = [
":__support_common",
- ":__support_cpp_algorithm",
":__support_cpp_bit",
+ ":__support_cpp_optional",
":__support_fputil_double_double",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
- ":__support_fputil_multiply_add",
":__support_fputil_nearest_integer",
- ":__support_fputil_polyeval",
+ ":__support_fputil_rounding_mode",
":__support_fputil_sqrt",
":__support_fputil_triple_double",
+ ":__support_macros_attributes",
":__support_macros_config",
":__support_macros_optimization",
":__support_math_common_constants",
":__support_math_exp10f",
":__support_math_exp2f",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_powf_double_eval",
+ hdrs = ["src/__support/math/powf_double_eval.h"],
+ deps = [
+ ":__support_common",
+ ":__support_cpp_algorithm",
+ ":__support_cpp_bit",
+ ":__support_fputil_double_double",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_manipulation_functions",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer",
+ ":__support_fputil_polyeval",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_common_constants",
":__support_math_exp_constants",
+ ":__support_math_powf_utils",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_powf_float_eval",
+ hdrs = ["src/__support/math/powf_float_eval.h"],
+ deps = [
+ ":__support_common",
+ ":__support_cpp_bit",
+ ":__support_fputil_double_double",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer",
+ ":__support_fputil_rounding_mode",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_common_constants",
+ ":__support_math_exp2f_float_utils",
+ ":__support_math_powf_utils",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_powf",
+ hdrs = ["src/__support/math/powf.h"],
+ deps = [
+ ":__support_common",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_powf_double_eval",
+ ":__support_math_powf_float_eval",
],
)
More information about the libc-commits
mailing list