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

[via libc-commits]

================
@@ -0,0 +1,153 @@
+//===-- 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/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 3 - 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 commands:
+// > 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
+LIBC_INLINE_VAR 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) {
----------------
lntue wrote:

Remove `constexpr` for now.

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