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

[Muhammad Bassiouni via llvm-commits]

================
@@ -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/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 {
+
+LIBC_INLINE static 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) {
----------------
bassiounix wrote:

mark as `constexpr`.

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