[libc-commits] [libc] [llvm] [libc][math][c++23] Add log10p1bf16 math function (PR #235343)

[via libc-commits]

================
@@ -0,0 +1,161 @@
+//===-- Implementation header for log10p1bf16 -------------------*- 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_LOG10P1BF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_LOG10P1BF16_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+namespace log10p1bf16_internal {
+
+// The top five mantissa bits select F = 1 + f * 2^-5, where f is in [0, 31].
+// BFloat16 has seven explicit mantissa bits, so the table uses five of those
+// bits for lookup and the remaining two bits to refine the residual.
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for f from 0 to 31 do print(round(log10(1 + f * 2^-5), SG, RN));
+LIBC_INLINE_VAR constexpr float LOG10F_F[32] = {
+    0x0p+0f,        0x1.b5e908p-7f, 0x1.af5f92p-6f, 0x1.3ed11ap-5f,
+    0x1.a30a9ep-5f, 0x1.02428cp-4f, 0x1.31b306p-4f, 0x1.5fe804p-4f,
+    0x1.8cf184p-4f, 0x1.b8de4ep-4f, 0x1.e3bc1ap-4f, 0x1.06cbd6p-3f,
+    0x1.1b3e72p-3f, 0x1.2f3b6ap-3f, 0x1.42c7e8p-3f, 0x1.55e8c6p-3f,
+    0x1.68a288p-3f, 0x1.7af974p-3f, 0x1.8cf184p-3f, 0x1.9e8e7cp-3f,
+    0x1.afd3e4p-3f, 0x1.c0c514p-3f, 0x1.d1653p-3f,  0x1.e1b734p-3f,
+    0x1.f1bdeep-3f, 0x1.00be06p-2f, 0x1.087a08p-2f, 0x1.101432p-2f,
+    0x1.178da6p-2f, 0x1.1ee778p-2f, 0x1.2622bp-2f,  0x1.2d404cp-2f,
+};
+
+// ONE_OVER_F_F[f] stores the reciprocal of the same representative F. It is
+// used to form the residual v = M * (1 / F) - 1.
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > for f from 0 to 31 do print(round(1 / (1 + f * 2^-5), SG, RN));
+LIBC_INLINE_VAR constexpr float ONE_OVER_F_F[32] = {
+    0x1p+0f,        0x1.f07c2p-1f,  0x1.e1e1e2p-1f, 0x1.d41d42p-1f,
+    0x1.c71c72p-1f, 0x1.bacf92p-1f, 0x1.af286cp-1f, 0x1.a41a42p-1f,
+    0x1.99999ap-1f, 0x1.8f9c18p-1f, 0x1.861862p-1f, 0x1.7d05f4p-1f,
+    0x1.745d18p-1f, 0x1.6c16c2p-1f, 0x1.642c86p-1f, 0x1.5c9882p-1f,
+    0x1.555556p-1f, 0x1.4e5e0ap-1f, 0x1.47ae14p-1f, 0x1.414142p-1f,
+    0x1.3b13b2p-1f, 0x1.3521dp-1f,  0x1.2f684cp-1f, 0x1.29e412p-1f,
+    0x1.24924ap-1f, 0x1.1f7048p-1f, 0x1.1a7b96p-1f, 0x1.15b1e6p-1f,
+    0x1.111112p-1f, 0x1.0c9714p-1f, 0x1.08421p-1f,  0x1.041042p-1f,
+};
+
+// Generated by Sollya with the following commands:
+//   > display = hexadecimal;
+//   > round(log10(2), SG, RN);
+LIBC_INLINE_VAR constexpr float LOG10F_2 = 0x1.344136p-2f;
+
+} // namespace log10p1bf16_internal
+
+LIBC_INLINE bfloat16 log10p1bf16(bfloat16 x) {
+  using namespace log10p1bf16_internal;
+  using FPBits = fputil::FPBits<bfloat16>;
+  FPBits x_bits(x);
+  uint16_t x_u = x_bits.uintval();
+  uint16_t x_abs = x_u & 0x7fffU;
+
+  // Handle special values; for |x| <= 0.125 (bits <= 0x3e00), use a direct
+  // polynomial to avoid cancellation in the range-reduction residual.
+  if (LIBC_UNLIKELY(x_abs <= 0x3e00U || x_abs >= 0x7f80U)) {
+    // log10p1(NaN) = NaN
+    if (x_bits.is_nan()) {
+      if (x_bits.is_signaling_nan()) {
+        fputil::raise_except_if_required(FE_INVALID);
+        return FPBits::quiet_nan().get_val();
+      }
+      return x;
+    }
+
+    // log10p1(+/-0) = +/-0
+    if (x_abs == 0U)
+      return x;
+
+    // log10p1(+inf) = +inf
+    if (x_u == 0x7f80U)
+      return x;
+
+    // log10p1(-inf) = NaN
+    if (x_abs >= 0x7f80U) {
----------------
Sukumarsawant wrote:

this is `>= inf` should be `==-inf` by your logic 

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