[libc][math] Implement C23 half precision pow function #159906
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
Thank you for submitting a Pull Request (PR) to the LLVM Project! This PR will be automatically labeled and the relevant teams will be notified. If you wish to, you can add reviewers by using the "Reviewers" section on this page. If this is not working for you, it is probably because you do not have write permissions for the repository. In which case you can instead tag reviewers by name in a comment by using If you have received no comments on your PR for a week, you can request a review by "ping"ing the PR by adding a comment “Ping”. The common courtesy "ping" rate is once a week. Please remember that you are asking for valuable time from other developers. If you have further questions, they may be answered by the LLVM GitHub User Guide. You can also ask questions in a comment on this PR, on the LLVM Discord or on the forums. |
|
@llvm/pr-subscribers-backend-risc-v @llvm/pr-subscribers-libc Author: Aِl-Mi'raj (AnonMiraj) ChangesThis patch adds half precision
Special cases (NaN, inf, zero, negative base with non-integer) are handled. CC: @lntue , @overmighty Closes #132198 Patch is 33.42 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/159906.diff 15 Files Affected:
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index 4b6f3337036aa..7dd5118665456 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -451,6 +451,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.nextupl
libc.src.math.pow
libc.src.math.powf
+ libc.src.math.powf16
libc.src.math.remainder
libc.src.math.remainderf
libc.src.math.remainderl
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index d24cc740d4234..9aae829cd4004 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -451,6 +451,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.nextupl
libc.src.math.pow
libc.src.math.powf
+ libc.src.math.powf16
libc.src.math.remainder
libc.src.math.remainderf
libc.src.math.remainderl
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 1bc5df9d45a93..584323d086079 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -589,6 +589,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.nextupl
libc.src.math.pow
libc.src.math.powf
+ libc.src.math.powf16
libc.src.math.remainder
libc.src.math.remainderf
libc.src.math.remainderl
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index ec01030c77d4f..a5759d1901cb8 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -405,6 +405,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.nextupl
libc.src.math.pow
libc.src.math.powf
+ libc.src.math.powf16
libc.src.math.remainder
libc.src.math.remainderf
libc.src.math.remainderl
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 54ea983d64839..b19d1033f2b75 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -599,6 +599,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.nextupl
libc.src.math.pow
libc.src.math.powf
+ libc.src.math.powf16
libc.src.math.remainder
libc.src.math.remainderf
libc.src.math.remainderl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 1fc9a2b901c14..03785ce2f8779 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -629,6 +629,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.nextupl
libc.src.math.pow
libc.src.math.powf
+ libc.src.math.powf16
libc.src.math.remainder
libc.src.math.remainderf
libc.src.math.remainderl
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index 4e398676bf91e..ff1608e3ec3dd 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -2083,6 +2083,12 @@ functions:
arguments:
- type: float
- type: float
+ - name: powf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
- name: powi
standards: llvm_libc_ext
return_type: double
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 187bc92e5c2c6..3a1ab17341265 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -451,6 +451,7 @@ add_math_entrypoint_object(nextupf128)
add_math_entrypoint_object(pow)
add_math_entrypoint_object(powf)
+add_math_entrypoint_object(powf16)
add_math_entrypoint_object(powi)
add_math_entrypoint_object(powif)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 419f562713b5f..fa5e0ed5b8451 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1581,6 +1581,30 @@ add_entrypoint_object(
libc.src.__support.math.expxf16_utils
)
+add_entrypoint_object(
+ powf16
+ SRCS
+ powf16.cpp
+ HDRS
+ ../pow.h
+ DEPENDS
+ .common_constants
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.optimization
+ libc.src.__support.math.expxf16_utils
+ COMPILE_OPTIONS
+ -O3
+)
+
add_entrypoint_object(
powf
SRCS
diff --git a/libc/src/math/generic/powf16.cpp b/libc/src/math/generic/powf16.cpp
new file mode 100644
index 0000000000000..8e2e6b552632f
--- /dev/null
+++ b/libc/src/math/generic/powf16.cpp
@@ -0,0 +1,296 @@
+//===-- Half-precision x^y function ---------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/powf16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/types.h"
+#include "src/__support/math/expxf16_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace {
+
+bool is_odd_integer(float16 x) {
+ using FPBits = fputil::FPBits<float16>;
+ FPBits xbits(x);
+ uint16_t x_u = xbits.uintval();
+ unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
+ unsigned lsb = static_cast<unsigned>(
+ cpp::countr_zero(static_cast<uint32_t>(x_u | FPBits::EXP_MASK)));
+ constexpr unsigned UNIT_EXPONENT =
+ static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
+ return (x_e + lsb == UNIT_EXPONENT);
+}
+
+bool is_integer(float16 x) {
+ using FPBits = fputil::FPBits<float16>;
+ FPBits xbits(x);
+ uint16_t x_u = xbits.uintval();
+ unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
+ unsigned lsb = static_cast<unsigned>(
+ cpp::countr_zero(static_cast<uint32_t>(x_u | FPBits::EXP_MASK)));
+ constexpr unsigned UNIT_EXPONENT =
+ static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
+ return (x_e + lsb >= UNIT_EXPONENT);
+}
+
+} // namespace
+
+LLVM_LIBC_FUNCTION(float16, powf16, (float16 x, float16 y)) {
+ using FPBits = fputil::FPBits<float16>;
+
+ FPBits xbits(x), ybits(y);
+ bool x_sign = xbits.is_neg();
+ bool y_sign = ybits.is_neg();
+
+ FPBits x_abs = xbits.abs();
+ FPBits y_abs = ybits.abs();
+
+ uint16_t x_u = xbits.uintval();
+ uint16_t x_a = x_abs.uintval();
+ uint16_t y_a = y_abs.uintval();
+ bool result_sign = false;
+
+ ///////// BEGIN - Check exceptional cases ////////////////////////////////////
+ // If x or y is signaling NaN
+ if (xbits.is_signaling_nan() || ybits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ //
+ if (LIBC_UNLIKELY(ybits.is_zero() || x_u == FPBits::one().uintval() ||
+ x_u >= FPBits::inf().uintval() ||
+ x_u < FPBits::min_normal().uintval())) {
+ // pow(x, 0) = 1
+ if (ybits.is_zero()) {
+ return fputil::cast<float16>(1.0f);
+ }
+
+ // pow(1, Y) = 1
+ if (x_u == FPBits::one().uintval()) {
+ return fputil::cast<float16>(1.0f);
+ }
+
+ switch (y_a) {
+
+ case 0x3800U: { // y = +-0.5
+ if (LIBC_UNLIKELY(
+ (x == 0.0 || x_u == FPBits::inf(Sign::NEG).uintval()))) {
+ // pow(-0, 1/2) = +0
+ // pow(-inf, 1/2) = +inf
+ // Make sure it works correctly for FTZ/DAZ.
+ return fputil::cast<float16>(y_sign ? (1.0 / (x * x)) : (x * x));
+ }
+ return fputil::cast<float16>(y_sign ? (1.0 / fputil::sqrt<float16>(x))
+ : fputil::sqrt<float16>(x));
+ }
+ case 0x3c00U: // y = +-1.0
+ return fputil::cast<float16>(y_sign ? (1.0 / x) : x);
+
+ case 0x4000U: // y = +-2.0
+ return fputil::cast<float16>(y_sign ? (1.0 / (x * x)) : (x * x));
+ }
+ // TODO: Speed things up with pow(2, y) = exp2(y) and pow(10, y) = exp10(y).
+
+ // Propagate remaining quiet NaNs.
+ if (xbits.is_quiet_nan()) {
+ return x;
+ }
+ if (ybits.is_quiet_nan()) {
+ return y;
+ }
+
+ // x = -1: special case for integer exponents
+ if (x_u == FPBits::one(Sign::NEG).uintval()) {
+ if (is_integer(y)) {
+ if (is_odd_integer(y)) {
+ return fputil::cast<float16>(-1.0f);
+ } else {
+ return fputil::cast<float16>(1.0f);
+ }
+ }
+ // pow(-1, non-integer) = NaN
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // x = 0 cases
+ if (xbits.is_zero()) {
+ if (y_sign) {
+ // pow(+-0, negative) = +-inf and raise FE_DIVBYZERO
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ bool result_neg = x_sign && ybits.is_finite() && is_odd_integer(y);
+ return FPBits::inf(result_neg ? Sign::NEG : Sign::POS).get_val();
+ } else {
+ // pow(+-0, positive) = +-0
+ bool out_is_neg = x_sign && is_odd_integer(y);
+ return out_is_neg ? FPBits::zero(Sign::NEG).get_val()
+ : FPBits::zero(Sign::POS).get_val();
+ }
+ }
+
+ if (xbits.is_inf()) {
+ bool out_is_neg = x_sign && ybits.is_finite() && is_odd_integer(y);
+ if (y_sign) // pow(+-inf, negative) = +-0
+ return out_is_neg ? FPBits::zero(Sign::NEG).get_val()
+ : FPBits::zero(Sign::POS).get_val();
+ // pow(+-inf, positive) = +-inf
+ return FPBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
+ }
+
+ // y = +-inf cases
+ if (ybits.is_inf()) {
+ // pow(1, inf) handled above.
+ bool x_abs_less_than_one = x_a < FPBits::one().uintval();
+ if ((x_abs_less_than_one && !y_sign) ||
+ (!x_abs_less_than_one && y_sign)) {
+ return fputil::cast<float16>(0.0f);
+ } else {
+ return FPBits::inf().get_val();
+ }
+ }
+
+ // pow( negative, non-integer ) = NaN
+ if (x_sign && !is_integer(y)) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // For negative x with integer y, compute pow(|x|, y) and adjust sign
+ if (x_sign) {
+ x = -x;
+ if (is_odd_integer(y)) {
+ result_sign = true;
+ }
+ }
+ }
+ ///////// END - Check exceptional cases //////////////////////////////////////
+
+ // x^y = 2^( y * log2(x) )
+ // = 2^( y * ( e_x + log2(m_x) ) )
+ // First we compute log2(x) = e_x + log2(m_x)
+
+ using namespace math::expxf16_internal;
+ FPBits x_bits(x);
+
+ uint16_t x_u_log = x_bits.uintval();
+
+ // Extract exponent field of x.
+ int m = -FPBits::EXP_BIAS;
+
+ // When x is subnormal, normalize it by multiplying by 2^FRACTION_LEN.
+ if ((x_u_log & FPBits::EXP_MASK) == 0U) {
+ constexpr float NORMALIZE_EXP =
+ static_cast<float>(1U << FPBits::FRACTION_LEN);
+ x_bits = FPBits(x_bits.get_val() * fputil::cast<float16>(NORMALIZE_EXP));
+ x_u_log = x_bits.uintval();
+ m -= FPBits::FRACTION_LEN;
+ }
+
+ // Extract the mantissa and index into small lookup tables.
+ uint16_t mant = x_bits.get_mantissa();
+ // Use the highest 5 fractional bits of the mantissa as the index f.
+ int f = mant >> 5;
+
+ m += (x_u_log >> FPBits::FRACTION_LEN);
+
+ // Add the hidden bit to the mantissa.
+ // 1 <= m_x < 2
+ x_bits.set_biased_exponent(FPBits::EXP_BIAS);
+ float mant_f = x_bits.get_val();
+
+ // Range reduction for log2(m_x):
+ // v = r * m_x - 1, where r is a power of 2 from a lookup table.
+ // The computation is exact for half-precision, and -2^-5 <= v < 2^-4.
+ // Then m_x = (1 + v) / r, and log2(m_x) = log2(1 + v) - log2(r).
+
+ float v = fputil::multiply_add(mant_f, ONE_OVER_F_F[f], -1.0f);
+
+ // For half-precision accuracy, we use a degree-2 polynomial approximation:
+ // P(v) ~ log2(1 + v) / v
+ // Generated by Sollya with:
+ // > P = fpminimax(log2(1+x)/x, 2, [|D...|], [-2^-5, 2^-4]);
+ // The coefficients are rounded from the Sollya output.
+ float log2p1_d_over_f =
+ v * fputil::polyeval(v, 0x1.715476p+0f, -0x1.71771ap-1f, 0x1.ecb38ep-2f);
+
+ // log2(1.mant) = log2(f) + log2(1 + v)
+ float log2_1_mant = LOG2F_F[f] + log2p1_d_over_f;
+
+ // Complete log2(x) = e_x + log2(m_x)
+ float log2_x = static_cast<float>(m) + log2_1_mant;
+
+ // z = y * log2(x)
+ // Now compute 2^z = 2^(n + r), with n integer and r in [-0.5, 0.5].
+ double z = fputil::cast<double>(y) * log2_x;
+
+ // Check for overflow/underflow for half-precision.
+ // Half-precision range is approximately 2^-24 to 2^15.
+ if (z > 15.0) {
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+ }
+ if (z < -24.0) {
+ fputil::raise_except_if_required(FE_UNDERFLOW);
+ return fputil::cast<float16>(0.0f);
+ }
+
+ double n = fputil::nearest_integer(z);
+ double r = z - n;
+
+ // Compute 2^r using a degree-7 polynomial for r in [-0.5, 0.5].
+ // Generated by Sollya with:
+ // > P = fpminimax(2^x, 7, [|D...|], [-0.5, 0.5]);
+ // The polynomial coefficients are rounded from the Sollya output.
+ constexpr double EXP2_COEFFS[] = {
+ 0x1p+0, // 1.0
+ 0x1.62e42fefa39efp-1, // ln(2)
+ 0x1.ebfbdff82c58fp-3, // ln(2)^2 / 2
+ 0x1.c6b08d704a0c0p-5, // ln(2)^3 / 6
+ 0x1.3b2ab6fba4e77p-7, // ln(2)^4 / 24
+ 0x1.5d87fe78a6737p-10, // ln(2)^5 / 120
+ 0x1.430912f86a805p-13, // ln(2)^6 / 720
+ 0x1.10e4104ac8015p-17 // ln(2)^7 / 5040
+ };
+
+ double exp2_r = fputil::polyeval(
+ r, EXP2_COEFFS[0], EXP2_COEFFS[1], EXP2_COEFFS[2], EXP2_COEFFS[3],
+ EXP2_COEFFS[4], EXP2_COEFFS[5], EXP2_COEFFS[6], EXP2_COEFFS[7]);
+
+ // Compute 2^n by direct bit manipulation.
+ int n_int = static_cast<int>(n);
+ uint64_t exp_bits = static_cast<uint64_t>(n_int + 1023) << 52;
+ double pow2_n = cpp::bit_cast<double>(exp_bits);
+
+ float16 result = fputil::cast<float16>(pow2_n * exp2_r);
+
+ if (result_sign) {
+ FPBits result_bits(result);
+ result_bits.set_sign(Sign::NEG);
+ result = result_bits.get_val();
+ }
+
+ return result;
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/powf16.h b/libc/src/math/powf16.h
new file mode 100644
index 0000000000000..bd50f16edeede
--- /dev/null
+++ b/libc/src/math/powf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for powf16 ------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_POWF16_H
+#define LLVM_LIBC_SRC_MATH_POWF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 powf16(float16 x, float16 y);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_POWF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index e15df147c3c35..eec048e094bd8 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2470,6 +2470,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ powf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ powf16_test.cpp
+ DEPENDS
+ libc.src.math.powf16
+ libc.src.__support.FPUtil.fp_bits
+)
add_fp_unittest(
atan2f_test
NEED_MPFR
diff --git a/libc/test/src/math/powf16_test.cpp b/libc/test/src/math/powf16_test.cpp
new file mode 100644
index 0000000000000..030f0e9cf8be3
--- /dev/null
+++ b/libc/test/src/math/powf16_test.cpp
@@ -0,0 +1,118 @@
+//===-- Unittests for powf16 ----------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/powf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcPowF16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+using LIBC_NAMESPACE::testing::tlog;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+TEST_F(LlvmLibcPowF16Test, TrickyInputs) {
+ // These values are in half precision.
+ constexpr mpfr::BinaryInput<float16> INPUTS[] = {
+ {static_cast<float16>(0x1.08p-2f), static_cast<float16>(0x1.0cp-1f)},
+ {static_cast<float16>(0x1.66p-1f), static_cast<float16>(0x1.f1p+1f)},
+ {static_cast<float16>(0x1.c04p-1f), static_cast<float16>(0x1.2p+12f)},
+ {static_cast<float16>(0x1.aep-1f), static_cast<float16>(0x1.f9p-1f)},
+ {static_cast<float16>(0x1.ffcp-1f), static_cast<float16>(0x1.fffp-2f)},
+ {static_cast<float16>(0x1.f55p-1f), static_cast<float16>(0x1.88p+12f)},
+ {static_cast<float16>(0x1.e84p-1f), static_cast<float16>(0x1.2cp+13f)},
+ };
+
+ for (auto input : INPUTS) {
+ float16 x = input.x;
+ float16 y = input.y;
+ EXPECT_MPFR_MATCH(mpfr::Operation::Pow, input, LIBC_NAMESPACE::powf16(x, y),
+ 1.0); // 1 ULP tolerance is enough for f16
+ }
+}
+
+TEST_F(LlvmLibcPowF16Test, InFloat16Range) {
+ constexpr uint16_t X_COUNT = 63;
+ constexpr uint16_t X_START = FPBits(static_cast<float16>(0.25)).uintval();
+ constexpr uint16_t X_STOP = FPBits(static_cast<float16>(4.0)).uintval();
+ constexpr uint16_t X_STEP = (X_STOP - X_START) / X_COUNT;
+
+ constexpr uint16_t Y_COUNT = 59;
+ constexpr uint16_t Y_START = FPBits(static_cast<float16>(0.25)).uintval();
+ constexpr uint16_t Y_STOP = FPBits(static_cast<float16>(4.0)).uintval();
+ constexpr uint16_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;
+ float16 mx = 0.0, my = 0.0, mr = 0.0;
+ double tol = 1.0; // start with 1 ULP for half precision
+
+ for (uint16_t i = 0, v = X_START; i <= X_COUNT; ++i, v += X_STEP) {
+ float16 x = FPBits(v).get_val();
+ if (FPBits(x).is_inf_or_nan() || x < static_cast<float16>(0.0))
+ continue;
+
+ for (uint16_t j = 0, w = Y_START; j <= Y_COUNT; ++j, w += Y_STEP) {
+ float16 y = FPBits(w).get_val();
+ if (FPBits(y).is_inf_or_nan())
+ continue;
+
+ float16 result = LIBC_NAMESPACE::powf16(x, y);
+ ++cc;
+ if (FPBits(result).is_inf_or_nan())
+ continue;
+
+ ++count;
+ mpfr::BinaryInput<float16> inputs{x, y};
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Pow, inputs,
+ result, 1.0, 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 > 128.0) // half precision is only ~11 bits
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ }
+ if (fails || (count < cc)) {
+ tlog << " powf16 failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n"
+ << " Max ULPs is at most: " << static_cast<uint64_t>(tol) ...
[truncated]
|
lntue
reviewed
Sep 21, 2025
lntue
reviewed
Sep 21, 2025
lntue
reviewed
Sep 21, 2025
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This patch adds half precision
powf16(x, y).We use the formula
log2(x)is computed with exponent + small table + short polynomial. computed in float2^zis done by splitting into integer part (direct exponent bits) and small remainder with degree-7 polynomial. computed in doubleSpecial cases (NaN, inf, zero, negative base with non-integer) are handled.
CC: @lntue , @overmighty
Closes #132198