[libc-commits] [libc] [libc][math] Implement C23 half precision pow function (PR #159906)

[via libc-commits]

================
@@ -0,0 +1,397 @@
+//===-- 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/common_constants.h"
+#include "src/__support/math/exp10f_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+using namespace common_constants_internal;
+namespace {
+static 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;
+}
+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();
+  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 ||
----------------
lntue wrote:

Add comments on these `y_a` values.

https://github.com/llvm/llvm-project/pull/159906