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

[via libc-commits]

================
@@ -0,0 +1,196 @@
+//===-- 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/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/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/log.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+// 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 interval)
+LIBC_INLINE_VAR constexpr float LGAMMA_POLY[7][5] = {
+    // [1,2), center=1.5, max_relative_err=1.29e-04
+    {-0x1.eeb280p-4f, 0x1.2f128ap-5f, 0x1.de1488p-2f, -0x1.2d373cp-3f,
+     0x1.08d8aep-4f},
+    // [2,3), center=2.5, max_rel_err=9.98e-06
+    {0x1.2383fcp-2f, 0x1.6809bep-1f, 0x1.f61322p-3f, -0x1.48c82ap-5f,
+     0x1.3b3cccp-7f},
+    // [3,4), center=3.5, max_rel_err=2.06e-06
+    {0x1.337302p+0f, 0x1.1a6912p+0f, 0x1.52475cp-3f, -0x1.2a2876p-6f,
+     0x1.859be8p-9f},
+    // [4,5), center=4.5, max_rel_err=6.61e-07
+    {0x1.3a140ap+1f, 0x1.638d3ep+0f, 0x1.fd62a0p-4f, -0x1.51efeep-7f,
+     0x1.4e023cp-10f},
+    // [5,6), center=5.5, max_rel_err=2.72e-07
+    {0x1.fa99a6p+1f, 0x1.9c70aep+0f, 0x1.984080p-4f, -0x1.b21838p-8f,
+     0x1.58a5b2p-11f},
+    // [6,7), center=6.5, max_rel_err=1.32e-07
+    {0x1.6a676ap+2f, 0x1.cafc46p+0f, 0x1.548cdap-4f, -0x1.2e0b0cp-8f,
+     0x1.909642p-12f},
+    // [7,8), center=7.5, max_rel_err=7.10e-08
+    {0x1.e23306p+2f, 0x1.f25eb8p+0f, 0x1.2413bep-4f, -0x1.bc5850p-9f,
+     0x1.f9d4bep-13f},
+};
+
+// lgamma_positive_d: compute lgamma(x) for x > 0, returning double.
+//
+// Takes double so callers can pass (1.0 + ax) without float precision loss.
+//
+// For x < 8, applies the recurrence lgamma(x) = lgamma(x+1) - ln(x) until
+// x reaches [4, 8), then evaluates the polynomial. This is critical because
+// the [2,3) polynomial has max_rel_err=9.98e-6 which, near its edges (t near
+// ±0.5), causes ~6e-6 absolute error. After subtracting ln(x), the result
+// near x=1 or x=2 can be as small as ~0.002, giving ~0.45 ULP error -- which
+// fails the directed-rounding tolerance test. Polynomials for [4,5) and above
+// have max_rel_err <= 6.61e-7, keeping final error well under 0.1 ULP.
----------------
Sukumarsawant wrote:

Haven't addressed previous comment regarding this .

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