[libc-commits] [libc] [llvm] [libc][math] Impl bfloat16 lgamma function. (PR #199312)
via libc-commits
libc-commits at lists.llvm.org
Wed Jul 22 07:45:56 PDT 2026
================
@@ -0,0 +1,261 @@
+//===-- 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/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.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 {
+
+namespace lgammabf16_internal {
+
+// lgamma_positive_d: compute lgamma(x) for x > 0, returning double.
+//
+// Takes double so callers can pass (1.0 + ax) without float precision loss.
+// The double return is necessary to avoid double-rounding: the final
+// fputil::cast<bfloat16> needs the full double precision so it can correctly
+// break ties (see the 0x35E5 example in lgammabf16 below).
+//
+// For x < 4, 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 or -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.
+//
+LIBC_INLINE double lgamma_positive_d(double x) {
+ // Coefficients for lgamma on [n, n+1), centered at n+0.5.
+ // Each row: {c0, c1, c2, c3, c4} for Estrin evaluation
+ // where t = x - (n + 0.5), n = 1..7.
+ // Maximum relative errors per interval.
+ // lgamma(x) on [n, n+1), n = 1..7, as P(t) where t = x - (n+0.5).
+ // P is a degree-4 fit to lgamma(t + n+0.5) on t in [-0.5, 0.5],
+ // generated via mpmath (no native lgamma in Sollya):
+ // > mpmath.mp.dps = 50
+ // > f = lambda t: mpmath.loggamma(t + (n + 0.5))
+ // > coeffs = mpmath.chebyfit(f, [-0.5, 0.5], 5)
+ // Reversed to ascending [c0..c4] and rounded to float32 for Estrin.
+ // Max relative error: ~1.29e-4 (n=1) to ~7.10e-8 (n=7), skipping
+ // |lgamma(x)| < 0.01 where error explodes near the zeros at x=1,2.
+ 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},
+ };
+
+ if (LIBC_UNLIKELY(x == 1.0 || x == 2.0))
+ return 0.0;
+
+ if (x >= 8.0) {
+ // Stirling series; 0.5*ln(2*pi)
+ constexpr double HALF_LN_2PI = 0x1.d67f1c864beb4p-1;
+ double lx = math::log(x);
+ double x2 = x * x;
+ double result = (x - 0.5) * lx - x + HALF_LN_2PI;
+ result += 1.0 / (12.0 * x) - 1.0 / (360.0 * x * x2);
+ return result;
+ }
+
+ // For x in (0, 4): apply recurrence relation
+ // lgamma(x) = lgamma(x+n) - ln(x*(x+1)*...*(x+n-1))
+ // to shift x into the stable [4, 8) range for polynomial evaluation.
+ double log_product, xs, product;
+
+ if (x >= 3.0) {
+ log_product = math::log(x);
+ xs = x + 1.0;
+ } else if (x >= 2.0) {
+ product = x * (x + 1.0);
+ log_product = math::log(product);
+ xs = x + 2.0;
+ } else if (x >= 1.0) {
+ product = x * (x + 1.0);
+ product = product * (x + 2.0);
+ log_product = math::log(product);
+ xs = x + 3.0;
+ } else {
+ product = x * (x + 1.0);
+ product = product * (x + 2.0);
+ product = product * (x + 3.0);
+ log_product = math::log(product);
+ xs = x + 4.0;
+ }
+
+ // xs in [4, 8); select polynomial interval.
+ float xf = static_cast<float>(xs);
+ int n = static_cast<int>(xf);
+ if (n >= 7)
+ n = 7;
+ float t = xf - (static_cast<float>(n) + 0.5f);
+
+ float c0 = LGAMMA_POLY[n - 1][0];
+ float c1 = LGAMMA_POLY[n - 1][1];
+ float c2 = LGAMMA_POLY[n - 1][2];
+ float c3 = LGAMMA_POLY[n - 1][3];
+ float c4 = LGAMMA_POLY[n - 1][4];
+
+ // Estrin's scheme for p(t) = c0 + c1*t + c2*t^2 + c3*t^3 + c4*t^4:
+ // p(t) = (c0 + c1*t) + t^2 * ((c2 + c3*t) + t^2 * c4)
+ float t2 = t * t;
+ float p01 = fputil::multiply_add(t, c1, c0);
+ float p23 = fputil::multiply_add(t, c3, c2);
+ float p234 = fputil::multiply_add(t2, c4, p23);
+ float lgamma_xs_f = fputil::multiply_add(t2, p234, p01);
+ double lgamma_xs = static_cast<double>(lgamma_xs_f);
+
+ return lgamma_xs - log_product;
+}
+
----------------
Sukumarsawant wrote:
```suggestion
```
https://github.com/llvm/llvm-project/pull/199312
More information about the libc-commits
mailing list