[libc-commits] [libc] [llvm] [libc][math][c23] Add sinhbf16 function (PR #224645)

[via libc-commits]

================
@@ -0,0 +1,171 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// General sinh/cosh functions for bfloat16.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINHBF16COSHBF16_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINHBF16COSHBF16_UTILS_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/macros/attributes.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace sinhbf16coshbf16_internal {
+
+LIBC_INLINE_VAR constexpr float LOG2F_E = 0x1.715476p+0f;
+
+LIBC_INLINE_VAR constexpr float LOGF_2 = 0x1.62e43p-1f;
+
+LIBC_INLINE_VAR constexpr cpp::array<uint32_t, 32> EXP2_MID_5_BITS = {
+    0x3f80'0000U, 0x3f82'cd87U, 0x3f85'aac3U, 0x3f88'980fU, 0x3f8b'95c2U,
+    0x3f8e'a43aU, 0x3f91'c3d3U, 0x3f94'f4f0U, 0x3f98'37f0U, 0x3f9b'8d3aU,
+    0x3f9e'f532U, 0x3fa2'7043U, 0x3fa5'fed7U, 0x3fa9'a15bU, 0x3fad'583fU,
+    0x3fb1'23f6U, 0x3fb5'04f3U, 0x3fb8'fbafU, 0x3fbd'08a4U, 0x3fc1'2c4dU,
+    0x3fc5'672aU, 0x3fc9'b9beU, 0x3fce'248cU, 0x3fd2'a81eU, 0x3fd7'44fdU,
+    0x3fdb'fbb8U, 0x3fe0'ccdfU, 0x3fe5'b907U, 0x3fea'c0c7U, 0x3fef'e4baU,
+    0x3ff5'257dU, 0x3ffa'83b3U,
+};
+
+// This function correctly calculates sinh(x) and cosh(x) by calculating exp(x)
+// and exp(-x) simultaneously.
+// To compute e^x, we perform the following range reduction:
+// find hi, mid, lo such that:
+//   x = (hi + mid) * log(2) + lo, in which
+//     hi is an integer,
+//     0 <= mid * 2^5 < 32 is an integer
+//     -2^(-5) <= lo * log2(e) <= 2^-5.
+// In particular,
+//   hi + mid = round(x * log2(e) * 2^5) * 2^(-5).
+// Then,
+//   e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.
+// We store 2^mid in the lookup table EXP2_MID_5_BITS, and compute 2^hi * 2^mid
+// by adding hi to the exponent field of 2^mid.
+// e^lo is computed using a degree-3 minimax polynomial generated by Sollya:
+//   e^lo ~ P(lo)
+//        = 1 + lo + c2 * lo^2 + ... + c5 * lo^5
+//        = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)
+//        = P_even + lo * P_odd
+// To compute e^(-x), notice that:
+//   e^(-x) = 2^(-(hi + mid)) * e^(-lo)
+//          ~ 2^(-(hi + mid)) * P(-lo)
+//          = 2^(-(hi + mid)) * (P_even - lo * P_odd)
+// So:
+//   sinh(x) = (e^x - e^(-x)) / 2
+//           ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -
+//                    2^(-(hi + mid)) * (P_even - lo * P_odd))
+//           = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +
+//                    lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))
+// And similarly:
+//   cosh(x) = (e^x + e^(-x)) / 2
+//           ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +
+//                    lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))
+// The main point of these formulas is that the expensive part of calculating
+// the polynomials approximating lower parts of e^x and e^(-x) is shared and
+// only done once.
+template <bool IsSinh> LIBC_INLINE float eval_sinh_or_cosh(float xf) {
+  float kf = fputil::nearest_integer(xf * (LOG2F_E * 0x1.0p+5f));
+  int x_hi_mid_p = static_cast<int>(kf);
+  int x_hi_mid_m = -x_hi_mid_p;
+
+  unsigned x_hi_p = static_cast<unsigned>(x_hi_mid_p) >> 5;
+  unsigned x_hi_m = static_cast<unsigned>(x_hi_mid_m) >> 5;
+  unsigned x_mid_p = static_cast<unsigned>(x_hi_mid_p) & 0x1f;
+  unsigned x_mid_m = static_cast<unsigned>(x_hi_mid_m) & 0x1f;
+
+  uint32_t exp2_hi_mid_bits_p =
+      EXP2_MID_5_BITS[x_mid_p] +
+      static_cast<uint32_t>(x_hi_p << fputil::FPBits<float>::FRACTION_LEN);
+  uint32_t exp2_hi_mid_bits_m =
+      EXP2_MID_5_BITS[x_mid_m] +
+      static_cast<uint32_t>(x_hi_m << fputil::FPBits<float>::FRACTION_LEN);
+  // exp2_hi_mid_p = 2^(hi + mid)
+  float exp2_hi_mid_p = fputil::FPBits<float>(exp2_hi_mid_bits_p).get_val();
+  // exp2_hi_mid_m = 2^(-(hi + mid))
+  float exp2_hi_mid_m = fputil::FPBits<float>(exp2_hi_mid_bits_m).get_val();
+
+  // exp2_hi_mid_sum = 2^(hi + mid) + 2^(-(hi + mid))
+  float exp2_hi_mid_sum = exp2_hi_mid_p + exp2_hi_mid_m;
+  // exp2_hi_mid_diff = 2^(hi + mid) - 2^(-(hi + mid))
+  float exp2_hi_mid_diff = exp2_hi_mid_p - exp2_hi_mid_m;
+
+  // lo = x - (hi + mid) = round(x * log2(e) * 2^5) * log(2) * (-2^(-5)) + x
+  float lo = fputil::multiply_add(kf, LOGF_2 * -0x1.0p-5f, xf);
+  float lo_sq = lo * lo;
+
+  // Degree-3 minimax polynomial generated by Sollya with the following
+  // commands:
+  //   > display = hexadecimal;
+  //   > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
+  //   > 1 + x * P;
+  constexpr cpp::array<float, 4> COEFFS = {0x1p+0f, 0x1p+0f, 0x1.0004p-1f,
+                                           0x1.555778p-3f};
+  float half_p_odd =
----------------
lntue wrote:

The polynomial is small enough, maybe you can just use Estrin's scheme directly.

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