[libc-commits] [libc] [llvm] [libc][math] Impl bfloat16 lgamma function. (PR #199312)

[Kiriti Ponduri via libc-commits]

================
@@ -0,0 +1,195 @@
+//===-- Implementation of lgammabf16 ----------------------------*- 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_LGAMMABF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/bit.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/macros/config.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+// Evaluate degree-4 polynomial using Horner's method:
+// p(t) = c0 + t*(c1 + t*(c2 + t*(c3 + t*c4)))
+// coeffs = {c4, c3, c2, c1, c0}
+static inline float poly4(float t, const float c[5]) {
+  return c[4] + t * (c[3] + t * (c[2] + t * (c[1] + t * c[0])));
+}
+
+// Compute natural log of positive float x using range reduction.
+// x = 2^e * m, m in [1, 2)
+// ln(x) = e * ln(2) + ln(m)
+// ln(m) approximated by degree-5 poly centered at 1.5
+static inline float lgamma_logf(float x) {
+  // Coefficients for ln(m) on [1,2), centered at 1.5
+  // Generated by numpy.polyfit, max_err = 1.71e-05 (well within bf16 precision)
+  constexpr float LN_COEFFS[6] = {
+      0x1.ef26300000000p-6f,  // c5
+      -0x1.c2e22a0000000p-5f, // c4
+      0x1.929d980000000p-4f,  // c3
+      -0x1.c615ec0000000p-3f, // c2
+      0x1.5557200000000p-1f,  // c1
+      0x1.9f309c0000000p-2f,  // c0
+  };
+  constexpr float LN2 = 0x1.62e430p-1f;
+
+  // Extract exponent and mantissa
+  uint32_t bits = cpp::bit_cast<uint32_t>(x);
+  int e = static_cast<int>((bits >> 23) & 0xFFU) - 127;
+  // Set exponent to 127 (i.e. value in [1,2))
+  bits = (bits & 0x807FFFFFU) | (127U << 23);
+  float m = cpp::bit_cast<float>(bits);
+
+  float t = m - 1.5f;
+  float lnm =
+      LN_COEFFS[5] +
+      t * (LN_COEFFS[4] +
+           t * (LN_COEFFS[3] +
+                t * (LN_COEFFS[2] + t * (LN_COEFFS[1] + t * LN_COEFFS[0]))));
+  return static_cast<float>(e) * LN2 + lnm;
+}
+
+// Coefficients for lgamma on [n, n+1), centered at n+0.5
+// Each row: {c4, c3, c2, c1, c0} for poly in (x - (n+0.5))
+// Generated by numpy.polyfit with degree 4
+// Intervals: n = 1..7
+static constexpr float LGAMMA_POLY[7][5] = {
----------------
udaykiriti wrote:

can i procedd with constexpr or LIBC_CONSTEXPR

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