[libc-commits] [libc] [llvm] [libc][math] Add tgammabf16 math support for bfloat16. (PR #208689)

[Vedant Neve via libc-commits]

================
@@ -0,0 +1,179 @@
+//===-- Implementation header for tgammabf16 --------------------*- 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_TGAMMABF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TGAMMABF16_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/bfloat16.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp.h"
+#include "src/__support/math/log.h"
+#include "src/__support/math/sin.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+namespace tgammabf16_internal {
+
+LIBC_INLINE_VAR constexpr double PI = 0x1.921fb54442d18p+1;
+LIBC_INLINE_VAR constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
+
+// Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
+LIBC_INLINE_VAR constexpr double LANCZOS_COEFFS[9] = {
+    0.99999999999980993,  676.5203681218851,     -1259.1392167224028,
+    771.32342877765313,   -176.61502916214059,   12.507343278224757,
+    -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
+
+LIBC_INLINE bool is_negative_integer(double x) {
+  if (x > -1.0)
+    return false;
+  if (x <= -128.0)
+    return true;
+  int n = static_cast<int>(x);
+  return x == static_cast<double>(n);
+}
+
+} // namespace tgammabf16_internal
----------------
0bVdnt wrote:

I have shifted the constants inside the function, since they are not needed anywhere else, but the Lanczos evaluation requires double intermediates so we can retain sufficient precision before final hardware rounding down to bfloat16. If we drop the intermediate math to float, the accumulated rounding errors during the polynomial evaluation cause us to lose precision and fail the strict 1.0 ULP bounds.

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