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

[via libc-commits]

================
@@ -0,0 +1,240 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Implementation of the bfloat16 tgamma function.
+///
+//===----------------------------------------------------------------------===//
+
+#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/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/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/expf.h"
+#include "src/__support/math/sinpif.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+// Sets ERANGE, raises FE_OVERFLOW, and returns the overflow value
+// (+inf or max_normal) according to the sign and active IEEE 754 rounding mode.
+LIBC_INLINE bfloat16 tgammabf16_overflow(Sign sign) {
+  fputil::set_errno_if_required(ERANGE);
+  fputil::raise_except_if_required(FE_OVERFLOW);
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+  return fputil::FPBits<bfloat16>::inf(sign).get_val();
+#else
+  switch (fputil::quick_get_round()) {
+  case FE_TONEAREST:
+    return fputil::FPBits<bfloat16>::inf(sign).get_val();
+  case FE_UPWARD:
+    return sign == Sign::POS
+               ? fputil::FPBits<bfloat16>::inf(Sign::POS).get_val()
+               : fputil::FPBits<bfloat16>::max_normal(Sign::NEG).get_val();
+  case FE_DOWNWARD:
+    return sign == Sign::POS
+               ? fputil::FPBits<bfloat16>::max_normal(Sign::POS).get_val()
+               : fputil::FPBits<bfloat16>::inf(Sign::NEG).get_val();
+  case FE_TOWARDZERO:
+    return fputil::FPBits<bfloat16>::max_normal(sign).get_val();
+  default:
+    return fputil::FPBits<bfloat16>::max_normal(sign).get_val();
+  }
+#endif
+}
+
+// This is only called with x in [7.5, 135.5] or with the positive Lanczos
+// coefficient sum. Range reduction into [sqrt(2)/2, sqrt(2)] confines
+// r = (m - 1) / (m + 1) to [-0.172, 0.172], making the degree-13 atanh Taylor
+// series accurate to ~2^-41 (sufficient for bfloat16 requirements).
+LIBC_INLINE double tgammabf16_log(double x) {
+  using DoubleBits = fputil::FPBits<double>;
+  // The hexadecimal literals were generated with Sollya:
+  // > display = hexadecimal;
+  // > round(log(2), D, RN);
+  // > round(sqrt(2), D, RN);
+  constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
+  constexpr double SQRT_2 = 0x1.6a09e667f3bcdp0;
+
+  DoubleBits xbits(x);
+  int exponent =
+      static_cast<int>(xbits.get_biased_exponent()) - DoubleBits::EXP_BIAS;
+  xbits.set_biased_exponent(DoubleBits::EXP_BIAS);
+  double mantissa = xbits.get_val();
+  if (mantissa > SQRT_2) {
+    mantissa *= 0.5;
+    ++exponent;
+  }
+
+  double r = (mantissa - 1.0) / (mantissa + 1.0);
+  double r2 = r * r;
+  // Maclaurin series: log(mantissa) = log((1 + r) / (1 - r)) = 2 * atanh(r)
+  //                                 = 2 * r * (1 + r^2/3 + r^4/5 + ...).
+  // The coefficients 1/(2k + 1) for k = 1..6 were generated with Sollya:
+  // > for k from 1 to 6 do print(round(1 / (2 * k + 1), D, RN));
+  double p = fputil::polyeval(r2, 0x1.5555555555555p-2, 0x1.999999999999ap-3,
+                              0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4,
+                              0x1.745d1745d1746p-4, 0x1.3b13b13b13b14p-4);
+  return fputil::multiply_add(static_cast<double>(exponent), LOG_2,
+                              2.0 * r * fputil::multiply_add(r2, p, 1.0));
+}
+
+LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
+  using FPBits = fputil::FPBits<bfloat16>;
+  // The hexadecimal literals were generated with Sollya:
+  // > display = hexadecimal;
+  // > round(pi, D, RN);
+  // > round(log(sqrt(2 * pi)), D, RN);
+  constexpr double PI = 0x1.921fb54442d18p+1;
+  constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
+  // Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
+  // Reference: "A note on the computation of the convergent Lanczos complex
+  // Gamma approximation" by Paul Godfrey (2001).
+  // Original:
+  // https://web.archive.org/web/20060915161115/http://my.fit.edu/~gabdo/gamma.txt
+  // Mirror: http://www.mrob.com/pub/ries/lanczos-gamma.html
+  constexpr double LANCZOS_COEFFS[9] = {
----------------
lntue wrote:

Move the constexpr array out of the function body.

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