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

[via libc-commits]

================
@@ -0,0 +1,193 @@
+//===-- 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/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"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+// 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
+// Maximum relative error: ~1.71e-05 (well within bfloat16 precision)
+LIBC_INLINE float lgamma_logf(float x) {
+  constexpr float LN_COEFFS[6] = {
+      0x1.9f309cp-2f,  // c0
+      0x1.555720p-1f,  // c1
+      -0x1.c615ecp-3f, // c2
+      0x1.929d98p-4f,  // c3
+      -0x1.c2e22ap-5f, // c4
+      0x1.ef2630p-6f,  // c5
+  };
+  constexpr float LN2 = 0x1.62e430p-1f;
+
+  uint32_t bits = cpp::bit_cast<uint32_t>(x);
+  int e_tmp = (bits >> 23) & 0xFFU;
+
+  // Normalize subnormal float32 values
+  if (e_tmp == 0) {
+    x *= 0x1.0p23f; // Multiply by 2^23
+    bits = cpp::bit_cast<uint32_t>(x);
+    e_tmp = ((bits >> 23) & 0xFFU) - 23;
+  }
+
+  int e = e_tmp - 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;
+  return static_cast<float>(e) * LN2 +
+         fputil::polyeval(t, LN_COEFFS[0], LN_COEFFS[1], LN_COEFFS[2],
+                          LN_COEFFS[3], LN_COEFFS[4], LN_COEFFS[5]);
+}
+
+// Coefficients for lgamma on [n, n+1), centered at n+0.5
+// Each row: {c0, c1, c2, c3, c4} for fputil::polyeval in (x - (n+0.5))
+// Generated by numpy.polyfit with degree 4
+// Intervals: n = 1..7
+// (Maximum relative errors per intrvel)
----------------
Sukumarsawant wrote:

```suggestion
// (Maximum relative errors per intervel)
```
nit.


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