[libc-commits] [libc] [llvm] [libc][math] Implement C23 half precision pow function (PR #159906)
via libc-commits
libc-commits at lists.llvm.org
Fri Jun 26 05:26:59 PDT 2026
https://github.com/AnonMiraj updated https://github.com/llvm/llvm-project/pull/159906
>From 5f350776ccac624cb27bf4a4cae60637d46b9e68 Mon Sep 17 00:00:00 2001
From: anonmiraj <nabilmalek48 at gmail.com>
Date: Sat, 20 Sep 2025 07:58:40 +0300
Subject: [PATCH 1/2] [libc][math] Implement C23 half precision pow function
fix formatting
add missing argument in math.yml
remove optimzation option
make tests exaustive
update log poly to double
fix powf16 tests
update powf16 approach
update tests
fix build failure
handle LIBC_TARGET_CPU_HAS_FMA_DOUBLE
ensure f16 literals are used in smoke test
Add LIBC_INLINE and clarify special exponents
Refactor powf16 implementation to header-only in src/__support/math folder.
add to shared tests and fix formating
fix typo
remove static
guard float16 entrypoints
add todo
fix nits
remove constexpr
---
libc/config/gpu/amdgpu/entrypoints.txt | 1 +
libc/config/gpu/nvptx/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/include/math.yaml | 8 +
libc/shared/math.h | 1 +
libc/shared/math/powf16.h | 29 ++
libc/src/__support/math/CMakeLists.txt | 21 +
libc/src/__support/math/powf16.h | 418 ++++++++++++++++++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/generic/CMakeLists.txt | 12 +
libc/src/math/generic/powf16.cpp | 18 +
libc/src/math/powf16.h | 21 +
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 1 +
libc/test/src/math/CMakeLists.txt | 11 +
libc/test/src/math/exhaustive/CMakeLists.txt | 16 +
libc/test/src/math/exhaustive/powf16_test.cpp | 72 +++
libc/test/src/math/powf16_test.cpp | 113 +++++
libc/test/src/math/smoke/CMakeLists.txt | 10 +
libc/test/src/math/smoke/powf16_test.cpp | 210 +++++++++
.../llvm-project-overlay/libc/BUILD.bazel | 29 ++
23 files changed, 997 insertions(+)
create mode 100644 libc/shared/math/powf16.h
create mode 100644 libc/src/__support/math/powf16.h
create mode 100644 libc/src/math/generic/powf16.cpp
create mode 100644 libc/src/math/powf16.h
create mode 100644 libc/test/src/math/exhaustive/powf16_test.cpp
create mode 100644 libc/test/src/math/powf16_test.cpp
create mode 100644 libc/test/src/math/smoke/powf16_test.cpp
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index 60a98cf5c773f..d924ef7a3852b 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -617,6 +617,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.nextdownf16
libc.src.math.nexttowardf16
libc.src.math.nextupf16
+ libc.src.math.powf16
libc.src.math.remainderf16
libc.src.math.remquof16
libc.src.math.rintf16
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index 70bdbde97e868..86be49457f967 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -619,6 +619,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.nextdownf16
libc.src.math.nexttowardf16
libc.src.math.nextupf16
+ libc.src.math.powf16
libc.src.math.remainderf16
libc.src.math.remquof16
libc.src.math.rintf16
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 570f7666ba67b..97f9c76460fe3 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -801,6 +801,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
# clang-12 and after: https://godbolt.org/z/8ceT9454c
# libc.src.math.nexttowardf16
libc.src.math.nextupf16
+ libc.src.math.powf16
libc.src.math.remainderf16
libc.src.math.remquof16
libc.src.math.rintf16
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 2508de5cfe5a4..4fe7d4e6dc08a 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -817,6 +817,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.nextdownf16
libc.src.math.nexttowardf16
libc.src.math.nextupf16
+ libc.src.math.powf16
libc.src.math.remainderf16
libc.src.math.remquof16
libc.src.math.rintf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5ace968a739c8..596f5bbf8d82d 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -887,6 +887,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.nextdownf16
libc.src.math.nexttowardf16
libc.src.math.nextupf16
+ libc.src.math.powf16
libc.src.math.remainderf16
libc.src.math.remquof16
libc.src.math.rintf16
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index 5bd0d9a4da76d..5301b7f7070f7 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -2253,6 +2253,14 @@ functions:
arguments:
- type: float
- type: float
+ - name: powf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: powi
standards:
- llvm_libc_ext
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 2e8a574e4d042..6e1b840e62d9d 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -402,6 +402,7 @@
#include "math/nextupl.h"
#include "math/pow.h"
#include "math/powf.h"
+#include "math/powf16.h"
#include "math/remainder.h"
#include "math/remainderbf16.h"
#include "math/remainderf.h"
diff --git a/libc/shared/math/powf16.h b/libc/shared/math/powf16.h
new file mode 100644
index 0000000000000..a44ff80a5ab90
--- /dev/null
+++ b/libc/shared/math/powf16.h
@@ -0,0 +1,29 @@
+//===-- Shared powf16 function ----------------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_POWF16_H
+#define LLVM_LIBC_SHARED_MATH_POWF16_H
+
+#include "include/llvm-libc-macros/float16-macros.h"
+#include "shared/libc_common.h"
+
+#ifdef LIBC_TYPES_HAS_FLOAT16
+
+#include "src/__support/math/powf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::powf16;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_TYPES_HAS_FLOAT16
+
+#endif // LLVM_LIBC_SHARED_MATH_POWF16_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index a57086a31dc03..5f2f0240c60b7 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4261,6 +4261,27 @@ add_header_library(
libc.src.errno.errno
)
+
+add_header_library(
+ powf16
+ HDRS
+ powf16.h
+ DEPENDS
+ .common_constants
+ .exp10f_utils
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.optimization
+)
+
add_header_library(
expm1f
HDRS
diff --git a/libc/src/__support/math/powf16.h b/libc/src/__support/math/powf16.h
new file mode 100644
index 0000000000000..6c59c2d061b57
--- /dev/null
+++ b/libc/src/__support/math/powf16.h
@@ -0,0 +1,418 @@
+//===-- 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___SUPPORT_MATH_POWF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF16_H
+#include "include/llvm-libc-macros/float16-macros.h"
+
+#ifdef LIBC_TYPES_HAS_FLOAT16
+
+#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/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/types.h"
+#include "src/__support/math/common_constants.h"
+#include "src/__support/math/exp10f_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace powf16_impl {
+// TODO: mark as constexpr when nearest_integer issue is resolved
+LIBC_INLINE double exp2_range_reduced(double x) {
+ // k = round(x * 32) => (hi + mid) * 2^5
+ double kf = fputil::nearest_integer(x * 32.0);
+ int k = static_cast<int>(kf);
+ // dx = lo = x - (hi + mid) = x - k * 2^(-5)
+ double dx = fputil::multiply_add(-0x1.0p-5, kf, x); // -2^-5 * k + x
+
+ // hi = k >> MID_BITS
+ // exp_hi = hi shifted into double exponent field
+ int64_t hi = static_cast<int64_t>(k >> ExpBase::MID_BITS);
+ int64_t exp_hi = static_cast<int64_t>(
+ static_cast<uint64_t>(hi) << fputil::FPBits<double>::FRACTION_LEN);
+
+ // mh_bits = bits for 2^hi * 2^mid (lookup contains base bits for 2^mid)
+ int tab_index = k & ExpBase::MID_MASK; // mid index in [0, 31]
+ int64_t mh_bits = ExpBase::EXP_2_MID[tab_index] + exp_hi;
+
+ // mh = 2^(hi + mid)
+ double mh = fputil::FPBits<double>(static_cast<uint64_t>(mh_bits)).get_val();
+
+ // Degree-5 polynomial approximating (2^x - 1)/x generating by Sollya with:
+ // > P = fpminimax((2^x - 1)/x, 5, [|D...|], [-1/32. 1/32]);
+ constexpr double COEFFS[5] = {0x1.62e42fefa39efp-1, 0x1.ebfbdff8131c4p-3,
+ 0x1.c6b08d7061695p-5, 0x1.3b2b1bee74b2ap-7,
+ 0x1.5d88091198529p-10};
+
+ double dx_sq = dx * dx;
+ double c1 = fputil::multiply_add(dx, COEFFS[0], 1.0); // 1 + ln2*dx
+ double c2 =
+ fputil::multiply_add(dx, COEFFS[2], COEFFS[1]); // COEFF1 + COEFF2*dx
+ double c3 =
+ fputil::multiply_add(dx, COEFFS[4], COEFFS[3]); // COEFF3 + COEFF4*dx
+ double p = fputil::multiply_add(dx_sq, c3, c2); // c2 + c3*dx^2
+
+ // 2^x = 2^(hi+mid) * 2^dx
+ // ≈ mh * (1 + dx * P(dx))
+ // = mh + (mh * dx) * P(dx)
+ double result = fputil::multiply_add(p, dx_sq * mh, c1 * mh);
+
+ return result;
+}
+
+LIBC_INLINE 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);
+}
+
+LIBC_INLINE constexpr 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 powf16_impl
+
+// TODO : optimize by using float whenever possible
+LIBC_INLINE float16 powf16(float16 x, float16 y) {
+ using namespace powf16_impl;
+ using namespace common_constants_internal;
+ 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();
+ uint16_t y_u = ybits.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() || xbits.is_nan() ||
+ ybits.is_nan() || x_u == FPBits::one().uintval() ||
+ x_u == FPBits::zero().uintval() || x_u >= FPBits::inf().uintval() ||
+ y_u >= FPBits::inf().uintval() ||
+ x_u < FPBits::min_normal().uintval() || y_a == 0x3400U || // 0.25
+ y_a == 0x3800U || // 0.5
+ y_a == 0x3A00U || // 0.75
+ y_a == 0x3D00U || // 1.25
+ y_a == 0x3E00U || // 1.5
+ y_a == 0x4000U || // 2.0
+ y_a == 0x4100U || // 2.5
+ y_a == 0x4300U || // 3.5
+ is_integer(y))) {
+ // pow(x, 0) = 1
+ if (ybits.is_zero()) {
+ return 1.0f16;
+ }
+
+ // pow(1, Y) = 1
+ if (x_u == FPBits::one().uintval()) {
+ return 1.0f16;
+ }
+ // 4. Handle remaining NaNs
+ // pow(NaN, y) = NaN (for y != 0)
+ if (xbits.is_nan()) {
+ return x;
+ }
+ // pow(x, NaN) = NaN (for x != 1)
+ if (ybits.is_nan()) {
+ return y;
+ }
+ switch (y_a) {
+ case 0x3400U: // y = ±0.25 (1/4)
+ case 0x3800U: // y = ±0.5 (1/2)
+ case 0x3A00U: // y = ±0.75 (3/4)
+ case 0x3D00U: // y = ±1.25 (5/4)
+ case 0x3E00U: // y = ±1.5 (3/2)
+ case 0x4100U: // y = ±2.5 (5/2)
+ case 0x4300U: // y = ±3.5 (7/2)
+ {
+ if (xbits.is_zero()) {
+ if (y_sign) {
+ // pow(±0, negative) handled below
+ break;
+ } else {
+ // pow(±0, positive_fractional) = +0
+ return FPBits::zero(Sign::POS).get_val();
+ }
+ }
+
+ if (x_sign && !xbits.is_zero()) {
+ break; // pow(negative, non-integer) = NaN
+ }
+
+ double x_d = static_cast<double>(x);
+ double sqrt_x = fputil::sqrt<double>(x_d);
+ double fourth_root = fputil::sqrt<double>(sqrt_x);
+ double result_d = 0.0;
+
+ // Compute based on exponent value
+ switch (y_a) {
+ case 0x3400U: // 0.25 = x^(1/4)
+ result_d = fourth_root;
+ break;
+ case 0x3800U: // 0.5 = x^(1/2)
+ result_d = sqrt_x;
+ break;
+ case 0x3A00U: // 0.75 = x^(1/2) * x^(1/4)
+ result_d = sqrt_x * fourth_root;
+ break;
+ case 0x3D00U: // 1.25 = x * x^(1/4)
+ result_d = x_d * fourth_root;
+ break;
+ case 0x3E00U: // 1.5 = x * x^(1/2)
+ result_d = x_d * sqrt_x;
+ break;
+ case 0x4100U: // 2.5 = x^2 * x^(1/2)
+ result_d = x_d * x_d * sqrt_x;
+ break;
+ case 0x4300U: // 3.5 = x^3 * x^(1/2)
+ result_d = x_d * x_d * x_d * sqrt_x;
+ break;
+ }
+
+ result_d = y_sign ? (1.0 / result_d) : result_d;
+ return fputil::cast<float16>(result_d);
+ }
+ case 0x3c00U: // y = +-1.0
+ return fputil::cast<float16>(y_sign ? (1.0 / x) : x);
+
+ case 0x4000U: // y = +-2.0
+ double result_d = static_cast<double>(x) * static_cast<double>(x);
+ return fputil::cast<float16>(y_sign ? (1.0 / (result_d)) : (result_d));
+ }
+ // TODO: Speed things up with pow(2, y) = exp2(y) and pow(10, y) = exp10(y).
+ //
+ // pow(-1, y) for integer y
+ if (x_u == FPBits::one(Sign::NEG).uintval()) {
+ if (is_integer(y)) {
+ if (is_odd_integer(y)) {
+ return -1.0f16;
+ } else {
+ return 1.0f16;
+ }
+ }
+ // pow(-1, non-integer) = NaN
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // pow(±0, y) 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)) {
+ // |x| < 1 and y = +inf => 0.0
+ // |x| > 1 and y = -inf => 0.0
+ return 0.0f16;
+ } else {
+ // |x| > 1 and y = +inf => +inf
+ // |x| < 1 and y = -inf => +inf
+ return FPBits::inf(Sign::POS).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();
+ }
+
+ bool result_sign = false;
+ if (x_sign && is_integer(y)) {
+ result_sign = is_odd_integer(y);
+ }
+
+ if (is_integer(y)) {
+ double base = x_abs.get_val();
+ double res = 1.0;
+ int yi = static_cast<int>(y_abs.get_val());
+
+ // Fast exponentiation by squaring
+ while (yi > 0) {
+ if (yi & 1)
+ res *= base;
+ base *= base;
+ yi = yi >> 1;
+ }
+
+ if (y_sign) {
+ res = 1.0 / res;
+ }
+
+ if (result_sign) {
+ res = -res;
+ }
+
+ if (FPBits(fputil::cast<float16>(res)).is_inf()) {
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ res = result_sign ? -0x1.0p20 : 0x1.0p20;
+ }
+
+ float16 final_res = fputil::cast<float16>(res);
+ return final_res;
+ }
+ }
+
+ ///////// END - Check exceptional cases //////////////////////////////////////
+
+ // Core computation: x^y = 2^( y * log2(x) )
+ // We compute log2(x) = log(x) / log(2) using a polynomial approximation.
+
+ // The exponent part (m) is added later to get the final log(x).
+ FPBits x_bits(x);
+ uint16_t x_u_log = x_bits.uintval();
+
+ // Extract exponent field of x.
+ int m = x_bits.get_exponent();
+
+ // When x is subnormal, normalize it by adjusting m.
+ if ((x_u_log & FPBits::EXP_MASK) == 0U) {
+ unsigned leading_zeros =
+ cpp::countl_zero(static_cast<uint32_t>(x_u_log)) - (32 - 16);
+
+ constexpr unsigned SUBNORMAL_SHIFT_CORRECTION = 5;
+ unsigned shift = leading_zeros - SUBNORMAL_SHIFT_CORRECTION;
+
+ x_bits.set_mantissa(static_cast<uint16_t>(x_u_log << shift));
+
+ m = 1 - FPBits::EXP_BIAS - static_cast<int>(shift);
+ }
+
+ // Extract the mantissa and index into small lookup tables.
+ uint16_t mant = x_bits.get_mantissa();
+ // Use the highest 7 fractional bits of the mantissa as the index f.
+ int f = mant >> (FPBits::FRACTION_LEN - 7);
+
+ // Reconstruct the mantissa value m_x so it's in the range [1.0, 2.0).
+ x_bits.set_biased_exponent(FPBits::EXP_BIAS);
+ double mant_d = x_bits.get_val();
+ // Degree-5 polynomial approximation
+ // of log2 generated by Sollya with:
+ // > P = fpminimax(log2(1 + x)/x, 4, [|1, D...|], [-2^-8, 2^-7]);
+ constexpr double COEFFS[5] = {0x1.71547652b8133p0, -0x1.71547652d1e33p-1,
+ 0x1.ec70a098473dep-2, -0x1.7154c5ccdf121p-2,
+ 0x1.2514fd90a130ap-2};
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double v = fputil::multiply_add<double>(mant_d, RD[f], -1.0);
+#else
+ double c = fputil::FPBits<double>(fputil::FPBits<double>(mant_d).uintval() &
+ 0x3fff'e000'0000'0000)
+ .get_val();
+ double v = fputil::multiply_add(RD[f], mant_d - c, CD[f]);
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double extra_factor = static_cast<double>(m) + LOG2_R[f];
+ double vsq = v * v;
+ double c0 = fputil::multiply_add(v, COEFFS[0], 0.0);
+ double c1 = fputil::multiply_add(v, COEFFS[2], COEFFS[1]);
+ double c2 = fputil::multiply_add(v, COEFFS[4], COEFFS[3]);
+
+ double log2_x = fputil::polyeval(vsq, c0, c1, c2);
+
+ double y_d = fputil::cast<double>(y);
+ double z = fputil::multiply_add(y_d, log2_x, y_d * extra_factor);
+
+ // Check for underflow
+ // Float16 min normal is 2^-14, smallest subnormal is 2^-24
+ if (LIBC_UNLIKELY(z < -25.0)) {
+ fputil::raise_except_if_required(FE_UNDERFLOW);
+ return result_sign ? FPBits::zero(Sign::NEG).get_val()
+ : FPBits::zero(Sign::POS).get_val();
+ }
+
+ // Check for overflow
+ // Float16 max is ~2^16
+ double result_d = 0.0;
+ if (LIBC_UNLIKELY(z > 16.0)) {
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ result_d = result_sign ? -0x1.0p20 : 0x1.0p20;
+ } else {
+ result_d = exp2_range_reduced(z);
+ }
+
+ if (result_sign) {
+
+ result_d = -result_d;
+ }
+
+ float16 result = fputil::cast<float16>((result_d));
+ return result;
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_TYPES_HAS_FLOAT16
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 4823ceefd7475..69565644a92d4 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -497,6 +497,7 @@ add_math_entrypoint_object(nextupbf16)
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 5a5bd84213a55..32bcb46b36c04 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1374,6 +1374,18 @@ add_entrypoint_object(
libc.src.__support.math.expm1f16
)
+
+add_entrypoint_object(
+ powf16
+ SRCS
+ powf16.cpp
+ HDRS
+ ../powf16.h
+ DEPENDS
+ libc.src.__support.math.powf16
+ libc.src.errno.errno
+)
+
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..068423bcc2f8a
--- /dev/null
+++ b/libc/src/math/generic/powf16.cpp
@@ -0,0 +1,18 @@
+//===-- 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 "src/__support/math/powf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(float16, powf16, (float16 x, float16 y)) {
+ return math::powf16(x, y);
+}
+
+} // 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/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index 9b990dd5a5c84..9b39ab44b8fa9 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -407,6 +407,7 @@ add_fp_unittest(
libc.src.__support.math.nextafterf128
libc.src.__support.math.pow
libc.src.__support.math.powf
+ libc.src.__support.math.powf16
libc.src.__support.math.remainder
libc.src.__support.math.remainderbf16
libc.src.__support.math.remainderf
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 5877729ae7931..ffc962d9c5104 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -61,6 +61,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {
EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::shared::logbf16(1.0f16));
EXPECT_EQ(0L, LIBC_NAMESPACE::shared::llogbf16(1.0f16));
+ EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::powf16(1.0f16, 1.0f16));
EXPECT_FP_EQ(0x1.921fb6p+0f16, LIBC_NAMESPACE::shared::acosf16(0.0f16));
EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::shared::sinf16(0.0f16));
EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::shared::tanf16(0.0f16));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 31f6468e01110..c244c30e511a6 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2819,6 +2819,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/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index 2d1301d3a1e66..aca3704e7487c 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -394,6 +394,22 @@ add_fp_unittest(
-lpthread
)
+add_fp_unittest(
+ powf16_test
+ NO_RUN_POSTBUILD
+ NEED_MPFR
+ SUITE
+ libc_math_exhaustive_tests
+ SRCS
+ powf16_test.cpp
+ DEPENDS
+ .exhaustive_test
+ libc.src.math.powf16
+ libc.src.__support.FPUtil.fp_bits
+ LINK_LIBRARIES
+ -lpthread
+)
+
add_fp_unittest(
hypotf_test
NO_RUN_POSTBUILD
diff --git a/libc/test/src/math/exhaustive/powf16_test.cpp b/libc/test/src/math/exhaustive/powf16_test.cpp
new file mode 100644
index 0000000000000..61c99baa77508
--- /dev/null
+++ b/libc/test/src/math/exhaustive/powf16_test.cpp
@@ -0,0 +1,72 @@
+//===-- Exhaustive test 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 "exhaustive_test.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/powf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+struct Powf16Checker : public virtual LIBC_NAMESPACE::testing::Test {
+ using FloatType = float16;
+ using FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+ using StorageType = typename FPBits::StorageType;
+
+ uint64_t check(uint16_t x_start, uint16_t x_stop, uint16_t y_start,
+ uint16_t y_stop, mpfr::RoundingMode rounding) {
+ mpfr::ForceRoundingMode r(rounding);
+ if (!r.success)
+ return true;
+ uint16_t xbits = x_start;
+ uint64_t failed = 0;
+ do {
+ float16 x = FPBits(xbits).get_val();
+ uint16_t ybits = y_start;
+ do {
+ float16 y = FPBits(ybits).get_val();
+ mpfr::BinaryInput<float16> input{x, y};
+ bool correct = TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Pow, input, LIBC_NAMESPACE::powf16(x, y), 0.5,
+ rounding);
+ failed += (!correct);
+ } while (ybits++ < y_stop);
+ } while (xbits++ < x_stop);
+ return failed;
+ }
+};
+
+using LlvmLibcPowf16ExhaustiveTest =
+ LlvmLibcExhaustiveMathTest<Powf16Checker, 1 << 8>;
+
+// Range: x in [0, inf], y in [0, inf]
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7C00U;
+
+TEST_F(LlvmLibcPowf16ExhaustiveTest, PositiveRange) {
+ test_full_range_all_roundings(POS_START, POS_STOP, POS_START, POS_STOP);
+}
+
+// Range: x in [-0, -inf], y in [0, inf]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xFC00U;
+
+TEST_F(LlvmLibcPowf16ExhaustiveTest, NegativeBasePositiveExponent) {
+ test_full_range_all_roundings(NEG_START, NEG_STOP, POS_START, POS_STOP);
+}
+
+// Range: x in [0, inf], y in [-0, -inf]
+TEST_F(LlvmLibcPowf16ExhaustiveTest, PositiveBaseNegativeExponent) {
+ test_full_range_all_roundings(POS_START, POS_STOP, NEG_START, NEG_STOP);
+}
+
+// Range: x in [-0, -inf], y in [-0, -inf]
+TEST_F(LlvmLibcPowf16ExhaustiveTest, NegativeRange) {
+ test_full_range_all_roundings(NEG_START, NEG_STOP, NEG_START, NEG_STOP);
+}
diff --git a/libc/test/src/math/powf16_test.cpp b/libc/test/src/math/powf16_test.cpp
new file mode 100644
index 0000000000000..1610e66184026
--- /dev/null
+++ b/libc/test/src/math/powf16_test.cpp
@@ -0,0 +1,113 @@
+//===-- 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/__support/CPP/bit.h"
+#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 FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+static constexpr float16 SELECTED_VALS[] = {
+ 0.83984375f16, 1.414f16, 0.0625f16, 2.5f16,
+ 3.140625f16, 15.5f16, 2.f16, 3.25f16};
+
+// Test tricky inputs for selected x values against all possible y values.
+TEST_F(LlvmLibcPowF16Test, TrickyInput_SelectedX_AllY) {
+ for (float16 x_base : SELECTED_VALS) {
+ // Only test non-negative x_base
+ if (FPBits(x_base).is_neg())
+ continue;
+
+ // Loop through normal and subnormal values only (0x0001 to 0x7BFF)
+ for (uint16_t y_u = 1; y_u <= 0x7BFFU; ++y_u) {
+ float16 y_base = FPBits(y_u).get_val();
+
+ // Case 1: (+x, +y) - Standard positive case
+ mpfr::BinaryInput<float16> input1{x_base, y_base};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input1,
+ LIBC_NAMESPACE::powf16(x_base, y_base),
+ 0.5);
+
+ // Case 2: (+x, -y) - Always valid for positive x
+ float16 y_neg = -y_base;
+ mpfr::BinaryInput<float16> input2{x_base, y_neg};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input2,
+ LIBC_NAMESPACE::powf16(x_base, y_neg),
+ 0.5);
+ }
+
+ // Case 3: (-x, +y) - Only test with positive integer y values
+ for (int y_int = 1; y_int <= 2048; ++y_int) {
+ float16 y_val = static_cast<float16>(y_int);
+ float16 x_neg = -x_base;
+ mpfr::BinaryInput<float16> input{x_neg, y_val};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input,
+ LIBC_NAMESPACE::powf16(x_neg, y_val), 0.5);
+ }
+
+ // Case 4: (-x, -y) - Only test with negative integer y values
+ for (int y_int = -2048; y_int < 0; ++y_int) {
+ float16 y_val = static_cast<float16>(y_int);
+ float16 x_neg = -x_base;
+ mpfr::BinaryInput<float16> input{x_neg, y_val};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input,
+ LIBC_NAMESPACE::powf16(x_neg, y_val), 0.5);
+ }
+ }
+}
+
+// Test tricky inputs for selected y values against all possible x values.
+TEST_F(LlvmLibcPowF16Test, TrickyInput_SelectedY_AllX) {
+ for (float16 y_base : SELECTED_VALS) {
+ // Only test non-negative y_base
+ if (FPBits(y_base).is_neg())
+ continue;
+
+ // Loop through normal and subnormal values only (0x0001 to 0x7BFF)
+ for (uint16_t x_u = 1; x_u <= 0x7BFFU; ++x_u) {
+ float16 x_base = FPBits(x_u).get_val();
+
+ // Case 1: (+x, +y) - Standard positive case
+ mpfr::BinaryInput<float16> input1{x_base, y_base};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input1,
+ LIBC_NAMESPACE::powf16(x_base, y_base),
+ 0.5);
+
+ // Case 2: (+x, -y) - Always valid for positive x
+ float16 y_neg = -y_base;
+ mpfr::BinaryInput<float16> input2{x_base, y_neg};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input2,
+ LIBC_NAMESPACE::powf16(x_base, y_neg),
+ 0.5);
+ }
+
+ // Case 3: (-x, +y) - Only test with positive integer x values
+ for (int x_int = 1; x_int <= 2048; ++x_int) {
+ float16 x_val = static_cast<float16>(x_int);
+ float16 x_neg = -x_val;
+ mpfr::BinaryInput<float16> input{x_neg, y_base};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input,
+ LIBC_NAMESPACE::powf16(x_neg, y_base),
+ 0.5);
+ }
+
+ // Case 4: (-x, -y) - Only test with negative integer x values
+ for (int x_int = 1; x_int <= 2048; ++x_int) {
+ float16 x_val = static_cast<float16>(x_int);
+ float16 x_neg = -x_val;
+ float16 y_neg = -y_base;
+ mpfr::BinaryInput<float16> input{x_neg, y_neg};
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Pow, input,
+ LIBC_NAMESPACE::powf16(x_neg, y_neg), 0.5);
+ }
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 73d2265322c4d..d3ed4a9d1798e 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -5305,6 +5305,16 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ powf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ powf16_test.cpp
+ DEPENDS
+ libc.src.math.powf16
+ libc.src.__support.FPUtil.fp_bits
+)
add_fp_unittest(
totalorder_test
SUITE
diff --git a/libc/test/src/math/smoke/powf16_test.cpp b/libc/test/src/math/smoke/powf16_test.cpp
new file mode 100644
index 0000000000000..f74bb6e57877e
--- /dev/null
+++ b/libc/test/src/math/smoke/powf16_test.cpp
@@ -0,0 +1,210 @@
+//===-- 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 "hdr/fenv_macros.h"
+#include "src/math/powf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcPowF16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+using LIBC_NAMESPACE::fputil::testing::ForceRoundingMode;
+using LIBC_NAMESPACE::fputil::testing::RoundingMode;
+
+TEST_F(LlvmLibcPowF16Test, SpecialNumbers) {
+ constexpr float16 NEG_ODD_INTEGER = -3.0f16;
+ constexpr float16 NEG_EVEN_INTEGER = -6.0f16;
+ constexpr float16 NEG_NON_INTEGER = -1.5f16;
+ constexpr float16 POS_ODD_INTEGER = 5.0f16;
+ constexpr float16 POS_EVEN_INTEGER = 8.0f16;
+ constexpr float16 POS_NON_INTEGER = 1.5f16;
+ constexpr float16 ONE_HALF = 0.5f16;
+
+ for (int i = 0; i < N_ROUNDING_MODES; ++i) {
+
+ ForceRoundingMode __r(ROUNDING_MODES[i]);
+ if (!__r.success)
+ continue;
+
+ // powf16( sNaN, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ aNaN, LIBC_NAMESPACE::powf16(sNaN, NEG_ODD_INTEGER), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ aNaN, LIBC_NAMESPACE::powf16(sNaN, NEG_EVEN_INTEGER), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ aNaN, LIBC_NAMESPACE::powf16(sNaN, POS_ODD_INTEGER), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ aNaN, LIBC_NAMESPACE::powf16(sNaN, POS_EVEN_INTEGER), FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, ONE_HALF),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, zero),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, neg_zero),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, inf),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, neg_inf),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(sNaN, aNaN),
+ FE_INVALID);
+
+ // powf16( 0.0f16, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(zero, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::powf16(zero, NEG_ODD_INTEGER), FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::powf16(zero, NEG_EVEN_INTEGER), FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::powf16(zero, NEG_NON_INTEGER), FE_DIVBYZERO);
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(zero, POS_ODD_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(zero, POS_EVEN_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(zero, POS_NON_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(zero, ONE_HALF));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(zero, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(zero, neg_zero));
+ EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::powf16(zero, inf));
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::powf16(zero, neg_inf),
+ FE_DIVBYZERO);
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(zero, aNaN));
+
+ // powf16( -0.0f16, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(neg_zero, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ neg_inf, LIBC_NAMESPACE::powf16(neg_zero, NEG_ODD_INTEGER),
+ FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::powf16(neg_zero, NEG_EVEN_INTEGER), FE_DIVBYZERO);
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::powf16(neg_zero, NEG_NON_INTEGER), FE_DIVBYZERO);
+ EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::powf16(neg_zero, POS_ODD_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(neg_zero, POS_EVEN_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(neg_zero, POS_NON_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(neg_zero, ONE_HALF));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(neg_zero, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(neg_zero, neg_zero));
+ EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::powf16(neg_zero, inf));
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::powf16(neg_zero, neg_inf),
+ FE_DIVBYZERO);
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(neg_zero, aNaN));
+
+ // powf16( 1.0f16, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(1.0f16, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, neg_zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, 1.0f16));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, -1.0f16));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, NEG_ODD_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, NEG_EVEN_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, NEG_NON_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, POS_ODD_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, POS_EVEN_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, POS_NON_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, inf));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, neg_inf));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(1.0f16, aNaN));
+
+ // powf16( -1.0f16, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(-1.0f16, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, neg_zero));
+ EXPECT_FP_EQ(-1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, 1.0f16));
+ EXPECT_FP_EQ(-1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, -1.0f16));
+ EXPECT_FP_EQ(-1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, NEG_ODD_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, NEG_EVEN_INTEGER));
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(
+ LIBC_NAMESPACE::powf16(-1.0f16, NEG_NON_INTEGER), FE_INVALID);
+ EXPECT_FP_EQ(-1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, POS_ODD_INTEGER));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, POS_EVEN_INTEGER));
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(
+ LIBC_NAMESPACE::powf16(-1.0f16, POS_NON_INTEGER), FE_INVALID);
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, inf));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(-1.0f16, neg_inf));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(-1.0f16, aNaN));
+
+ // powf16( inf, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(inf, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(inf, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(inf, neg_zero));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(inf, 1.0f16));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(inf, -1.0f16));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(inf, NEG_ODD_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(inf, NEG_EVEN_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(inf, NEG_NON_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(inf, POS_ODD_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(inf, POS_EVEN_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(inf, POS_NON_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(inf, ONE_HALF));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(inf, inf));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(inf, neg_inf));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(inf, aNaN));
+
+ // powf16( -inf, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(neg_inf, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(neg_inf, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(neg_inf, neg_zero));
+ EXPECT_FP_EQ(neg_inf, LIBC_NAMESPACE::powf16(neg_inf, 1.0f16));
+ EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::powf16(neg_inf, -1.0f16));
+ EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::powf16(neg_inf, NEG_ODD_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(neg_inf, NEG_EVEN_INTEGER));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(neg_inf, NEG_NON_INTEGER));
+ EXPECT_FP_EQ(neg_inf, LIBC_NAMESPACE::powf16(neg_inf, POS_ODD_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(neg_inf, POS_EVEN_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(neg_inf, POS_NON_INTEGER));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(neg_inf, ONE_HALF));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(neg_inf, inf));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(neg_inf, neg_inf));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(neg_inf, aNaN));
+
+ // powf16 ( aNaN, exponent )
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::powf16(aNaN, sNaN),
+ FE_INVALID);
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(aNaN, zero));
+ EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::powf16(aNaN, neg_zero));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, 1.0f16));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, -1.0f16));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, NEG_ODD_INTEGER));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, NEG_EVEN_INTEGER));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, NEG_NON_INTEGER));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, POS_ODD_INTEGER));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, POS_EVEN_INTEGER));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, POS_NON_INTEGER));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, inf));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, neg_inf));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::powf16(aNaN, aNaN));
+
+ // powf16 ( base, inf )
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(0.1f16, inf));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(-0.1f16, inf));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(1.1f16, inf));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(-1.1f16, inf));
+
+ // powf16 ( base, -inf )
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(0.1f16, neg_inf));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::powf16(-0.1f16, neg_inf));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(1.1f16, neg_inf));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::powf16(-1.1f16, neg_inf));
+
+ // Overflow / Underflow:
+ if (ROUNDING_MODES[i] != RoundingMode::Downward &&
+ ROUNDING_MODES[i] != RoundingMode::TowardZero) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::powf16(3.1f16, 21.0f16),
+ FE_OVERFLOW);
+ }
+ if (ROUNDING_MODES[i] != RoundingMode::Upward) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ 0.0f16, LIBC_NAMESPACE::powf16(3.1f16, -21.0f16), FE_UNDERFLOW);
+ }
+ }
+}
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 68bb1f1d60bf2..bfe9593be45d7 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -8695,6 +8695,25 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_powf16",
+ hdrs = ["src/__support/math/powf16.h"],
+ deps = [
+ ":__support_cpp_bit",
+ ":__support_fputil_cast",
+ ":__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_macros_optimization",
+ ":__support_math_common_constants",
+ ":__support_math_exp10f_utils",
+ ],
+)
+
libc_support_library(
name = "__support_math_expm1f16",
hdrs = ["src/__support/math/expm1f16.h"],
@@ -12414,6 +12433,14 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "powf16",
+ additional_deps = [
+ ":__support_math_powf16",
+ ":errno",
+ ],
+)
+
libc_math_function(
name = "remainder",
additional_deps = [
@@ -12421,6 +12448,8 @@ libc_math_function(
],
)
+libc_math_function(name = "remainder")
+
libc_math_function(
name = "remainderbf16",
additional_deps = [
>From e388540a0e7a33b5d52f54dfb3d99340312e9370 Mon Sep 17 00:00:00 2001
From: Anonmiraj <ezzibrahimx at gmail.com>
Date: Fri, 26 Jun 2026 15:26:28 +0300
Subject: [PATCH 2/2] remove power by squraing and some optimzation
---
libc/src/__support/math/CMakeLists.txt | 3 +
libc/src/__support/math/powf16.h | 116 +++++++-----------
.../llvm-project-overlay/libc/BUILD.bazel | 8 +-
3 files changed, 51 insertions(+), 76 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 5f2f0240c60b7..75ec60cb12f9b 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4271,6 +4271,7 @@ add_header_library(
.exp10f_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.include.llvm-libc-macros.float16_macros
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
@@ -4279,7 +4280,9 @@ add_header_library(
libc.src.__support.FPUtil.nearest_integer
libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.config
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
)
add_header_library(
diff --git a/libc/src/__support/math/powf16.h b/libc/src/__support/math/powf16.h
index 6c59c2d061b57..beb24f63a0127 100644
--- a/libc/src/__support/math/powf16.h
+++ b/libc/src/__support/math/powf16.h
@@ -33,8 +33,7 @@ namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace powf16_impl {
-// TODO: mark as constexpr when nearest_integer issue is resolved
-LIBC_INLINE double exp2_range_reduced(double x) {
+LIBC_INLINE constexpr double exp2_range_reduced(double x) {
// k = round(x * 32) => (hi + mid) * 2^5
double kf = fputil::nearest_integer(x * 32.0);
int k = static_cast<int>(kf);
@@ -102,7 +101,9 @@ LIBC_INLINE constexpr bool is_integer(float16 x) {
} // namespace powf16_impl
-// TODO : optimize by using float whenever possible
+// 0.5 ULP correctly rounding of float16 needs the ~2^-52 precision
+// to survive the final double->float16 rounding. Computing them in float was
+// measured to round x^+-2.5 1 ULP off via double rounding.
LIBC_INLINE float16 powf16(float16 x, float16 y) {
using namespace powf16_impl;
using namespace common_constants_internal;
@@ -161,20 +162,20 @@ LIBC_INLINE float16 powf16(float16 x, float16 y) {
return y;
}
switch (y_a) {
- case 0x3400U: // y = ±0.25 (1/4)
- case 0x3800U: // y = ±0.5 (1/2)
- case 0x3A00U: // y = ±0.75 (3/4)
- case 0x3D00U: // y = ±1.25 (5/4)
- case 0x3E00U: // y = ±1.5 (3/2)
- case 0x4100U: // y = ±2.5 (5/2)
- case 0x4300U: // y = ±3.5 (7/2)
+ case 0x3400U: // y = +-0.25 (1/4)
+ case 0x3800U: // y = +-0.5 (1/2)
+ case 0x3A00U: // y = +-0.75 (3/4)
+ case 0x3D00U: // y = +-1.25 (5/4)
+ case 0x3E00U: // y = +-1.5 (3/2)
+ case 0x4100U: // y = +-2.5 (5/2)
+ case 0x4300U: // y = +-3.5 (7/2)
{
if (xbits.is_zero()) {
if (y_sign) {
- // pow(±0, negative) handled below
+ // pow(+-0, negative) handled below
break;
} else {
- // pow(±0, positive_fractional) = +0
+ // pow(+-0, positive_fractional) = +0
return FPBits::zero(Sign::POS).get_val();
}
}
@@ -213,18 +214,19 @@ LIBC_INLINE float16 powf16(float16 x, float16 y) {
break;
}
- result_d = y_sign ? (1.0 / result_d) : result_d;
- return fputil::cast<float16>(result_d);
+ return fputil::cast<float16>(y_sign ? (1.0 / result_d) : result_d);
}
case 0x3c00U: // y = +-1.0
- return fputil::cast<float16>(y_sign ? (1.0 / x) : x);
+ return fputil::cast<float16>(y_sign ? (1.0 / static_cast<double>(x))
+ : static_cast<double>(x));
- case 0x4000U: // y = +-2.0
- double result_d = static_cast<double>(x) * static_cast<double>(x);
- return fputil::cast<float16>(y_sign ? (1.0 / (result_d)) : (result_d));
+ case 0x4000U: { // y = +-2.0
+ double sq = static_cast<double>(x) * static_cast<double>(x);
+ return fputil::cast<float16>(y_sign ? (1.0 / sq) : sq);
+ }
}
// TODO: Speed things up with pow(2, y) = exp2(y) and pow(10, y) = exp10(y).
- //
+
// pow(-1, y) for integer y
if (x_u == FPBits::one(Sign::NEG).uintval()) {
if (is_integer(y)) {
@@ -240,7 +242,7 @@ LIBC_INLINE float16 powf16(float16 x, float16 y) {
return FPBits::quiet_nan().get_val();
}
- // pow(±0, y) cases
+ // pow(+-0, y) cases
if (xbits.is_zero()) {
if (y_sign) {
// pow(+-0, negative) = +-inf and raise FE_DIVBYZERO
@@ -287,39 +289,25 @@ LIBC_INLINE float16 powf16(float16 x, float16 y) {
return FPBits::quiet_nan().get_val();
}
- bool result_sign = false;
- if (x_sign && is_integer(y)) {
+ if (x_sign)
result_sign = is_odd_integer(y);
- }
- if (is_integer(y)) {
- double base = x_abs.get_val();
- double res = 1.0;
- int yi = static_cast<int>(y_abs.get_val());
-
- // Fast exponentiation by squaring
- while (yi > 0) {
- if (yi & 1)
+ if (!y_sign && is_integer(y)) {
+ int n = static_cast<int>(y_abs.get_val());
+ if (n >= 2 && n <= 7) {
+ double base = x_abs.get_val();
+ double res = 1.0;
+ if (n & 4)
+ res *= base;
+ res *= res;
+ if (n & 2)
+ res *= base;
+ res *= res;
+ if (n & 1)
res *= base;
- base *= base;
- yi = yi >> 1;
- }
-
- if (y_sign) {
- res = 1.0 / res;
- }
-
- if (result_sign) {
- res = -res;
- }
- if (FPBits(fputil::cast<float16>(res)).is_inf()) {
- fputil::raise_except_if_required(FE_OVERFLOW);
- res = result_sign ? -0x1.0p20 : 0x1.0p20;
+ return fputil::cast<float16>(result_sign ? -res : res);
}
-
- float16 final_res = fputil::cast<float16>(res);
- return final_res;
}
}
@@ -329,7 +317,7 @@ LIBC_INLINE float16 powf16(float16 x, float16 y) {
// We compute log2(x) = log(x) / log(2) using a polynomial approximation.
// The exponent part (m) is added later to get the final log(x).
- FPBits x_bits(x);
+ FPBits x_bits = x_abs;
uint16_t x_u_log = x_bits.uintval();
// Extract exponent field of x.
@@ -382,31 +370,13 @@ LIBC_INLINE float16 powf16(float16 x, float16 y) {
double y_d = fputil::cast<double>(y);
double z = fputil::multiply_add(y_d, log2_x, y_d * extra_factor);
- // Check for underflow
- // Float16 min normal is 2^-14, smallest subnormal is 2^-24
- if (LIBC_UNLIKELY(z < -25.0)) {
- fputil::raise_except_if_required(FE_UNDERFLOW);
- return result_sign ? FPBits::zero(Sign::NEG).get_val()
- : FPBits::zero(Sign::POS).get_val();
- }
-
- // Check for overflow
- // Float16 max is ~2^16
- double result_d = 0.0;
- if (LIBC_UNLIKELY(z > 16.0)) {
- fputil::raise_except_if_required(FE_OVERFLOW);
- result_d = result_sign ? -0x1.0p20 : 0x1.0p20;
- } else {
- result_d = exp2_range_reduced(z);
- }
-
- if (result_sign) {
+ if (LIBC_UNLIKELY(z > 16.0)) // |x^y| > 2^16 -> overflow
+ return fputil::cast<float16>(result_sign ? -0x1.0p20 : 0x1.0p20);
+ if (LIBC_UNLIKELY(z < -25.0)) // |x^y| < 2^-25 -> flush to zero
+ return fputil::cast<float16>(result_sign ? -0x1.0p-30 : 0x1.0p-30);
- result_d = -result_d;
- }
-
- float16 result = fputil::cast<float16>((result_d));
- return result;
+ double result_d = exp2_range_reduced(z);
+ return fputil::cast<float16>(result_sign ? -result_d : result_d);
}
} // namespace math
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index bfe9593be45d7..8d82dc4c899ba 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -8706,11 +8706,15 @@ libc_support_library(
":__support_fputil_multiply_add",
":__support_fputil_nearest_integer",
":__support_fputil_polyeval",
- ":__support_fputil_rounding_mode",
":__support_fputil_sqrt",
+ ":__support_macros_config",
":__support_macros_optimization",
+ ":__support_macros_properties_types",
":__support_math_common_constants",
":__support_math_exp10f_utils",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ":llvm_libc_macros_float16_macros",
],
)
@@ -12448,8 +12452,6 @@ libc_math_function(
],
)
-libc_math_function(name = "remainder")
-
libc_math_function(
name = "remainderbf16",
additional_deps = [
More information about the libc-commits
mailing list