[libc-commits] [libc] [llvm] [libc][math][c23] Add exp10bf16 math function (PR #193299)

[via libc-commits]

================
@@ -0,0 +1,147 @@
+//===-- Implementation header for exp10bf16 ---------------------*- 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_EXP10BF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10BF16_H
+
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/ManipulationFunctions.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/optimization.h"
+
+// This is an algorithm for exp10(x) in bfloat16 which is correctly rounded for
+// all rounding modes, based on the implementation of exp10bf16(x)
+// from the RLIBM project at: https://people.cs.rutgers.edu/~sn349/rlibm
+
+// Step 1 - Range reduction of edge cases:
+// Find the largest value x for which exp10bf16(x) = +0.0
+// x < log10(min subnormal bfloat16)
+// Find the smallest value x for which exp10bf16(x) = +INF
+// x > log10(max normal bfloat16)
+
+// Step 2 - Range reduction to subnormal range:
+// exp10(x) = exp2(x * log2(10)) and further decompose x * log2(10) into
+// sum of two parts - closest integer i and a fraction f
+// exp10(x) = exp2(x * log2(10)) = exp2(i + f) = exp2(i) * exp2(f)
+
+// Step 4 - Polynomial approximation:
+// Approximate the range reduced function exp2(f) on a domain (-0.5, 0.5)
+// using a fourth degree polynomial P
+// Generated by Sollya with the following command:
+// > display = hexadecimal;
+// > exp2 = 2^x;
+// > I = [-0.5, 0.5];
+// > P = fpminimax(exp2, 4, [|1, SG...|], I);
+// > E = infnorm(exp2 - P, I);
+
+// Step 4 - Output compensation:
+// exp10(x) = exp2(i) * P(f)
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+// Generated by Sollya with the following command:
+// round(log2(10), SG, RN);
+LIBC_INLINE_VAR constexpr float LOG2F_10 = 0x1.a934fp1;
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+constexpr fputil::ExceptValues<bfloat16, 4> EXP10BF16_EXCEPTS{
+    {// {inputs, RZ output, RU offset, RD offset, RN offset}
+     // x = -0x1.6ep-7, exp10bf16(x) = 0x1.f2p-1 (RN)
+     {0xBC37U, 0x3F79U, 1U, 0U, 0U},
+     // x = -0x1.2ap-6, exp10bf16(x) = 0x1.eap-1 (RN)
+     {0xBC95U, 0x3F75U, 1U, 0U, 0U},
+     // x = -0x1.74p-3, exp10bf16(x) = 0x1.5p-1 (RN)
+     {0xBE3AU, 0x3F28U, 1U, 0U, 0U},
+     // x = -0x1.44p-1, exp10bf16(x) = 0x1.dcp-3 (RN)
+     {0xBF22U, 0x3E6EU, 1U, 0U, 0U}}};
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+LIBC_INLINE constexpr bfloat16 exp10bf16(bfloat16 x) {
+  using FPBits = typename fputil::FPBits<bfloat16>;
+  FPBits x_bits(x);
+
+  uint16_t x_u = x_bits.uintval();
+
+  // special case where exp10bf16(NaN) = NaN
+  if (LIBC_UNLIKELY(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;
+  }
+
+  // exp10bf16(x) = +0.0
+  if (x_bits.is_neg() && x_u >= 0xc222U) { // x <= -0x1.44p+5
+    return FPBits::zero().get_val();
+  }
+
+  // exp10bf16(x) = +INF
+  if (x_bits.is_pos() && x_u >= 0x421bU) { // x >= 0x1.36p+5
+    if (x_bits.is_inf()) {
+      return FPBits::inf().get_val();
----------------
ProfessionalMenace wrote:

Thank you for your feedback! Fixed

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