[libc-commits] [libc] [llvm] [libc][math] Impl bfloat16 lgamma function. (PR #199312)
Kiriti Ponduri via libc-commits
libc-commits at lists.llvm.org
Mon Jun 1 23:53:14 PDT 2026
https://github.com/udaykiriti updated https://github.com/llvm/llvm-project/pull/199312
>From 8df441e478ca4336f383e1fb22786f94d2cbfc86 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Sat, 23 May 2026 07:56:37 +0530
Subject: [PATCH 01/10] [libc][math] Impl bfloat16 lgamma function.
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/config/baremetal/aarch64/entrypoints.txt | 1 +
libc/config/baremetal/arm/entrypoints.txt | 1 +
libc/config/baremetal/riscv/entrypoints.txt | 1 +
libc/config/darwin/aarch64/entrypoints.txt | 1 +
libc/config/darwin/x86_64/entrypoints.txt | 1 +
libc/config/gpu/amdgpu/entrypoints.txt | 1 +
libc/config/gpu/nvptx/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/arm/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/config/windows/entrypoints.txt | 1 +
libc/shared/math.h | 1 +
libc/shared/math/lgammabf16.h | 23 +++
libc/src/__support/math/CMakeLists.txt | 16 ++
libc/src/__support/math/lgammabf16.h | 193 ++++++++++++++++++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/generic/CMakeLists.txt | 12 ++
libc/src/math/generic/lgammabf16.cpp | 18 ++
libc/src/math/lgammabf16.h | 21 ++
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 3 +-
libc/test/src/math/CMakeLists.txt | 12 ++
libc/test/src/math/lgammabf16_test.cpp | 26 +++
libc/test/src/math/smoke/CMakeLists.txt | 12 ++
libc/test/src/math/smoke/lgammabf16_test.cpp | 32 +++
libc/utils/MPFRWrapper/MPCommon.cpp | 7 +
libc/utils/MPFRWrapper/MPCommon.h | 1 +
libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 +
libc/utils/MPFRWrapper/MPFRUtils.h | 1 +
.../llvm-project-overlay/libc/BUILD.bazel | 23 +++
31 files changed, 416 insertions(+), 1 deletion(-)
create mode 100644 libc/shared/math/lgammabf16.h
create mode 100644 libc/src/__support/math/lgammabf16.h
create mode 100644 libc/src/math/generic/lgammabf16.cpp
create mode 100644 libc/src/math/lgammabf16.h
create mode 100644 libc/test/src/math/lgammabf16_test.cpp
create mode 100644 libc/test/src/math/smoke/lgammabf16_test.cpp
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index dcb50135232e2..2c27769fe0786 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -835,6 +835,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index fac62bac939cc..78b0ee578f472 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -846,6 +846,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index a3b96225ff09d..9893ee44d4e18 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -843,6 +843,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index 914d2b7918da1..56abccb09dce7 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -656,6 +656,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/darwin/x86_64/entrypoints.txt b/libc/config/darwin/x86_64/entrypoints.txt
index b941c968255e8..bc164e4fcd3f6 100644
--- a/libc/config/darwin/x86_64/entrypoints.txt
+++ b/libc/config/darwin/x86_64/entrypoints.txt
@@ -278,6 +278,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index c76900a8371f7..913db74f49f12 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -684,6 +684,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index 6538732f785f5..1fdcd32ff10ee 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -686,6 +686,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index a50923fddec54..e7d42623cd56f 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -948,6 +948,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 49b30ef1830f3..8bac9164441d6 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -526,6 +526,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 51c769d853a52..5374674401b65 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -969,6 +969,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5545790fecd85..c3d284b758254 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -1034,6 +1034,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 94d1d00e676d9..5f14967e3c8a7 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -368,6 +368,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.iscanonicalbf16
libc.src.math.issignalingbf16
libc.src.math.ldexpbf16
+ libc.src.math.lgammabf16
libc.src.math.llogbbf16
libc.src.math.llrintbf16
libc.src.math.llroundbf16
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 2baeb07294a3b..9482efb1dc1b7 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -305,6 +305,7 @@
#include "math/ldexpf128.h"
#include "math/ldexpf16.h"
#include "math/ldexpl.h"
+#include "math/lgammabf16.h"
#include "math/llogb.h"
#include "math/llogbbf16.h"
#include "math/llogbf.h"
diff --git a/libc/shared/math/lgammabf16.h b/libc/shared/math/lgammabf16.h
new file mode 100644
index 0000000000000..85fdd03ff60ab
--- /dev/null
+++ b/libc/shared/math/lgammabf16.h
@@ -0,0 +1,23 @@
+//===-- Shared header for 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_SHARED_MATH_LGAMMABF16_H
+#define LLVM_LIBC_SHARED_MATH_LGAMMABF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/math/lgammabf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::lgammabf16;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_LGAMMABF16_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 9e3ec26cdc881..189d1d004abac 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -46,6 +46,22 @@ add_header_library(
libc.src.__support.macros.properties.types
)
+add_header_library(
+ lgammabf16
+ HDRS
+ lgammabf16.h
+ DEPENDS
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.bfloat16
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+)
+
add_header_library(
acosh_float_constants
HDRS
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
new file mode 100644
index 0000000000000..70694544fcc33
--- /dev/null
+++ b/libc/src/__support/math/lgammabf16.h
@@ -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)
+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},
+};
+
+// Compute lgamma for x >= 1 using polynomial or Stirling
+LIBC_INLINE float lgamma_positive(float x) {
+ if (LIBC_UNLIKELY(x == 1.0f || x == 2.0f))
+ return 0.0f;
+ // For x >= 8: Stirling approximation
+ // lgamma(x) ~ (x-0.5)*ln(x) - x + 0.5*ln(2*pi)
+ if (x >= 8.0f) {
+ constexpr float HALF_LN_2PI = 0x1.6840e6p-1f;
+ float lx = lgamma_logf(x);
+ // using FMA: (x - 0.5)*ln(x) - x + 0.5*ln(2*pi)
+ return fputil::multiply_add(x - 0.5f, lx, -x + HALF_LN_2PI);
+ }
+
+ // For x in [1, 8): use piecewise polynomial
+ int n = static_cast<int>(x); // floor for x >= 1
+ if (n >= 7)
+ n = 7;
+ float t = x - (static_cast<float>(n) + 0.5f);
+ const float *c = LGAMMA_POLY[n - 1];
+ // Co-efficients stored as {c0, c1, c2, c3, c4}
+ return fputil::polyeval(t, c[0], c[1], c[2], c[3], c[4]);
+}
+
+LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
+ using FPBits = fputil::FPBits<bfloat16>;
+ FPBits x_bits(x);
+
+ // Handle NaN
+ if (LIBC_UNLIKELY(x_bits.is_nan())) {
+ if (x_bits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ uint16_t x_u = x_bits.uintval();
+ uint16_t x_abs = x_u & 0x7fffU;
+
+ // +Inf or -Inf -> +Inf
+ if (LIBC_UNLIKELY(x_abs == 0x7f80U))
+ return FPBits::inf(Sign::POS).get_val();
+
+ // +-0 -> +Inf (pole error)
+ if (LIBC_UNLIKELY(x_abs == 0U)) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return FPBits::inf(Sign::POS).get_val();
+ }
+
+ float xf = static_cast<float>(x);
+
+ // Negative integers -> +Inf (pole error)
+ if (LIBC_UNLIKELY(x_bits.is_neg())) {
+ int biased_exp = x_abs >> FPBits::FRACTION_LEN;
+ if (biased_exp >= FPBits::EXP_BIAS) {
+ int e = biased_exp - FPBits::EXP_BIAS;
+ if (e >= FPBits::FRACTION_LEN ||
+ (x_bits.get_mantissa() &
+ static_cast<uint16_t>((1U << (FPBits::FRACTION_LEN - e)) - 1U)) ==
+ 0U) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return FPBits::inf(Sign::POS).get_val();
+ }
+ }
+
+ // Negative non-integer: reflection formula
+ // lgamma(x) = ln(pi) - ln|sin(pi*x)| - lgamma(1-x)
+ constexpr float LN_PI = 0x1.250d06p+0f; // ln(pi)
+ float ax = -xf;
+ float frac = ax - static_cast<float>(static_cast<int>(ax));
+
+ constexpr float A1 = 0x1.55553cp-2f; // pi^2/6 approx
+ constexpr float A2 = 0x1.11100ap-7f; // pi^4/120 approx
+ constexpr float PI = 0x1.921fb6p+1f;
+ float frac2 = frac * frac;
+ // sin(pi*frac) = PI * frac * polyeval(frac^2, 1, -A1, A2)
+ float sin_pi_frac = PI * frac * fputil::polyeval(frac2, 1.0f, -A1, A2);
+ if (sin_pi_frac < 0.0f)
+ sin_pi_frac = -sin_pi_frac;
+
+ float log_sin = lgamma_logf(sin_pi_frac);
+ float result = LN_PI - log_sin - lgamma_positive(1.0f + ax);
+ return fputil::cast<bfloat16>(result);
+ }
+
+ // Positive x in (0, 1): use recurrence lgamma(x) = lgamma(x+1) - ln(x)
+ if (xf < 1.0f) {
+ float lnx = lgamma_logf(xf);
+ return fputil::cast<bfloat16>(lgamma_positive(xf + 1.0f) - lnx);
+ }
+
+ // Positive x >= 1
+ return fputil::cast<bfloat16>(lgamma_positive(xf));
+}
+
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index b53817e2a1729..b50b848df999e 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -596,6 +596,7 @@ add_math_entrypoint_object(tanpif16)
add_math_entrypoint_object(tgamma)
add_math_entrypoint_object(tgammaf)
add_math_entrypoint_object(lgamma)
+add_math_entrypoint_object(lgammabf16)
add_math_entrypoint_object(lgamma_r)
add_math_entrypoint_object(totalorder)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 7ccbddba07b8d..c03843f2fcfb3 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4869,6 +4869,18 @@ add_entrypoint_object(
libc.src.__support.math.bf16fmal
)
+
+add_entrypoint_object(
+ lgammabf16
+ SRCS
+ lgammabf16.cpp
+ HDRS
+ ../lgammabf16.h
+ DEPENDS
+ libc.src.__support.math.lgammabf16
+ libc.src.__support.FPUtil.bfloat16
+)
+
add_entrypoint_object(
bf16fmaf128
SRCS
diff --git a/libc/src/math/generic/lgammabf16.cpp b/libc/src/math/generic/lgammabf16.cpp
new file mode 100644
index 0000000000000..a26e5fe437a90
--- /dev/null
+++ b/libc/src/math/generic/lgammabf16.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of lgammabf16 function -----------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/lgammabf16.h"
+#include "src/__support/math/lgammabf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(bfloat16, lgammabf16, (bfloat16 x)) {
+ return math::lgammabf16(x);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/lgammabf16.h b/libc/src/math/lgammabf16.h
new file mode 100644
index 0000000000000..965c8c881eb03
--- /dev/null
+++ b/libc/src/math/lgammabf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for 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_MATH_LGAMMABF16_H
+#define LLVM_LIBC_SRC_MATH_LGAMMABF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+bfloat16 lgammabf16(bfloat16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_LGAMMABF16_H
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index 84a8d4bf79b3d..237a7b1e209f1 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -297,6 +297,7 @@ add_fp_unittest(
libc.src.__support.math.ldexpbf16
libc.src.__support.math.ldexpf
libc.src.__support.math.ldexpl
+ libc.src.__support.math.lgammabf16
libc.src.__support.math.llogb
libc.src.__support.math.llogbbf16
libc.src.__support.math.llrint
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 634778380dc8e..c42c2ffd12f53 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -785,7 +785,8 @@ TEST(LlvmLibcSharedMathTest, AllBFloat16) {
EXPECT_FP_EQ(bfloat16(5.0),
LIBC_NAMESPACE::shared::hypotbf16(bfloat16(4.0), bfloat16(3.0)));
-
+ EXPECT_FP_EQ(bfloat16(0.0f),
+ LIBC_NAMESPACE::shared::lgammabf16(bfloat16(1.0f)));
EXPECT_FP_EQ(bfloat16(0.0), LIBC_NAMESPACE::shared::logbbf16(bfloat16(1.0f)));
bfloat16 setpayloadbf16_res = bfloat16(0.0);
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 4213c11eca515..ac7a1963ad6ce 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -30,6 +30,18 @@ add_fp_unittest(
FMA_OPT__ONLY
)
+add_fp_unittest(
+ lgammabf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ lgammabf16_test.cpp
+ DEPENDS
+ libc.src.math.lgammabf16
+ libc.src.__support.FPUtil.bfloat16
+)
+
add_fp_unittest(
cos_test
NEED_MPFR
diff --git a/libc/test/src/math/lgammabf16_test.cpp b/libc/test/src/math/lgammabf16_test.cpp
new file mode 100644
index 0000000000000..509df108158d9
--- /dev/null
+++ b/libc/test/src/math/lgammabf16_test.cpp
@@ -0,0 +1,26 @@
+//===-- Exhaustive test for lgammabf16 ------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/math/lgammabf16.h"
+#include "test/src/math/exhaustive/exhaustive_test.h"
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LlvmLibcLgammabf16ExhaustiveTest =
+ LlvmLibcUnaryOpExhaustiveMathTest<LIBC_NAMESPACE::fputil::BFloat16,
+ mpfr::Operation::Lgamma,
+ LIBC_NAMESPACE::lgammabf16>;
+
+TEST_F(LlvmLibcLgammabf16ExhaustiveTest, PositiveRange) {
+ test_full_range_all_roundings(0x0000U, 0x7F80U);
+}
+
+TEST_F(LlvmLibcLgammabf16ExhaustiveTest, NegativeRange) {
+ test_full_range_all_roundings(0x8000U, 0xFF80U);
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 28b85b1a25bbd..ad7387a52d41c 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -36,6 +36,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ lgammabf16_test
+ UNIT_TEST_ONLY
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ lgammabf16_test.cpp
+ DEPENDS
+ libc.src.math.lgammabf16
+ libc.src.__support.FPUtil.bfloat16
+)
+
add_fp_unittest(
cospif16_test
SUITE
diff --git a/libc/test/src/math/smoke/lgammabf16_test.cpp b/libc/test/src/math/smoke/lgammabf16_test.cpp
new file mode 100644
index 0000000000000..393cae08fd10c
--- /dev/null
+++ b/libc/test/src/math/smoke/lgammabf16_test.cpp
@@ -0,0 +1,32 @@
+//===-- Unittests for lgammabf16 ------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/math/lgammabf16.h"
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+class LlvmLibcLgammaBf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+ DECLARE_SPECIAL_CONSTANTS(bfloat16)
+public:
+ void test_special_numbers() {
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::lgammabf16(aNaN));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(zero));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(neg_zero));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(inf));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(neg_inf));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(1.0f)));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(2.0f)));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-1.0f)));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-2.0f)));
+ }
+};
+
+TEST_F(LlvmLibcLgammaBf16Test, SpecialNumbers) { test_special_numbers(); }
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index bdd8286148641..2d95148cc2a0b 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -379,6 +379,13 @@ MPFRNumber MPFRNumber::log1p() const {
return result;
}
+MPFRNumber MPFRNumber::lgamma() const {
+ MPFRNumber result(*this);
+ int signp;
+ mpfr_lgamma(result.value, &signp, value, mpfr_rounding);
+ return result;
+}
+
MPFRNumber MPFRNumber::pow(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_pow(result.value, value, b.value, mpfr_rounding);
diff --git a/libc/utils/MPFRWrapper/MPCommon.h b/libc/utils/MPFRWrapper/MPCommon.h
index 935a2614968a2..47bebffef9fec 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -218,6 +218,7 @@ class MPFRNumber {
MPFRNumber log10() const;
MPFRNumber log10p1() const;
MPFRNumber log1p() const;
+ MPFRNumber lgamma() const;
MPFRNumber pow(const MPFRNumber &b);
MPFRNumber remquo(const MPFRNumber &divisor, int "ient);
MPFRNumber round() const;
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 356764302bda8..2330ec05003b6 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -87,6 +87,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.log10p1();
case Operation::Log1p:
return mpfrInput.log1p();
+ case Operation::Lgamma:
+ return mpfrInput.lgamma();
case Operation::Mod2PI:
return mpfrInput.mod_2pi();
case Operation::ModPIOver2:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 6c1b15598b0f2..a9029b83c1e2f 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -54,6 +54,7 @@ enum class Operation : int {
Log10,
Log10p1,
Log1p,
+ Lgamma,
Mod2PI,
ModPIOver2,
ModPIOver4,
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 87963e4524b39..297f0fa2e4720 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -6233,6 +6233,22 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_lgammabf16",
+ hdrs = ["src/__support/math/lgammabf16.h"],
+ deps = [
+ ":__support_cpp_bit",
+ ":__support_fputil_bfloat16",
+ ":__support_fputil_cast",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ],
+)
+
libc_support_library(
name = "__support_math_llogbbf16",
hdrs = ["src/__support/math/llogbbf16.h"],
@@ -11397,6 +11413,13 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "lgammabf16",
+ additional_deps = [
+ ":__support_math_lgammabf16",
+ ],
+)
+
libc_math_function(
name = "llogb",
additional_deps = [
>From b4f798b71d6891e9778f26d85e70f34bb9f0fa54 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Mon, 25 May 2026 13:56:47 +0530
Subject: [PATCH 02/10] [libc][math] Implement lgammabf16 function
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 2 +
libc/test/src/math/lgammabf16_test.cpp | 73 ++++++++++++++++---
.../llvm-project-overlay/libc/BUILD.bazel | 2 +
3 files changed, 67 insertions(+), 10 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 189d1d004abac..792c77bea06ab 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -56,6 +56,8 @@ add_header_library(
libc.src.__support.FPUtil.bfloat16
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.PolyEval
+ libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.config
libc.src.__support.macros.optimization
libc.hdr.errno_macros
diff --git a/libc/test/src/math/lgammabf16_test.cpp b/libc/test/src/math/lgammabf16_test.cpp
index 509df108158d9..a137257381e84 100644
--- a/libc/test/src/math/lgammabf16_test.cpp
+++ b/libc/test/src/math/lgammabf16_test.cpp
@@ -1,4 +1,4 @@
-//===-- Exhaustive test for lgammabf16 ------------------------------------===//
+//===-- Unittests for lgammabf16 ------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
@@ -8,19 +8,72 @@
#include "src/__support/FPUtil/bfloat16.h"
#include "src/math/lgammabf16.h"
-#include "test/src/math/exhaustive/exhaustive_test.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLgammabf16Test = LIBC_NAMESPACE::testing::FPTest<bfloat16>;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
-using LlvmLibcLgammabf16ExhaustiveTest =
- LlvmLibcUnaryOpExhaustiveMathTest<LIBC_NAMESPACE::fputil::BFloat16,
- mpfr::Operation::Lgamma,
- LIBC_NAMESPACE::lgammabf16>;
+// Subnormal positive range
+static constexpr uint16_t SUBNORM_POS_START = 0x0001U;
+static constexpr uint16_t SUBNORM_POS_STOP = 0x007FU;
+
+// Normal positive range
+static constexpr uint16_t NORMAL_POS_START = 0x0080U;
+static constexpr uint16_t NORMAL_POS_STOP = 0x7F7FU;
+
+// Subnormal negative range
+static constexpr uint16_t SUBNORM_NEG_START = 0x8001U;
+static constexpr uint16_t SUBNORM_NEG_STOP = 0x807FU;
-TEST_F(LlvmLibcLgammabf16ExhaustiveTest, PositiveRange) {
- test_full_range_all_roundings(0x0000U, 0x7F80U);
+// Normal negative range
+static constexpr uint16_t NORMAL_NEG_START = 0x8080U;
+static constexpr uint16_t NORMAL_NEG_STOP = 0xFF7FU;
+
+TEST_F(LlvmLibcLgammabf16Test, SpecialNumbers) {
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::lgammabf16(aNaN));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(zero));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(neg_zero));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(inf));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(neg_inf));
+ // lgamma(1) = lgamma(2) = 0
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(1.0f)));
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(2.0f)));
+ // Negative integers are poles -> +inf
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-1.0f)));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-2.0f)));
}
-TEST_F(LlvmLibcLgammabf16ExhaustiveTest, NegativeRange) {
- test_full_range_all_roundings(0x8000U, 0xFF80U);
+TEST_F(LlvmLibcLgammabf16Test, SubnormalPositiveRange) {
+ for (uint16_t v = SUBNORM_POS_START; v <= SUBNORM_POS_STOP; ++v) {
+ bfloat16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammabf16(x), 0.5);
+ }
}
+
+TEST_F(LlvmLibcLgammabf16Test, NormalPositiveRange) {
+ for (uint16_t v = NORMAL_POS_START; v <= NORMAL_POS_STOP; ++v) {
+ bfloat16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammabf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcLgammabf16Test, SubnormalNegativeRange) {
+ for (uint16_t v = SUBNORM_NEG_START; v <= SUBNORM_NEG_STOP; ++v) {
+ bfloat16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammabf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcLgammabf16Test, NormalNegativeRange) {
+ for (uint16_t v = NORMAL_NEG_START; v <= NORMAL_NEG_STOP; ++v) {
+ bfloat16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammabf16(x), 0.5);
+ }
+}
\ No newline at end of file
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 297f0fa2e4720..55daa7e103e42 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -6242,6 +6242,8 @@ libc_support_library(
":__support_fputil_cast",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_polyeval",
":__support_macros_config",
":__support_macros_optimization",
":hdr_errno_macros",
>From 17978529af5438896b6c80c75b1492f893ad70ef Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Thu, 28 May 2026 11:09:06 +0530
Subject: [PATCH 03/10] [libc] Fix double-rounding in lgammabf16 for small
positive inputs
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 2 +-
libc/src/__support/math/lgammabf16.h | 159 ++++++++++---------
libc/src/math/generic/CMakeLists.txt | 2 -
libc/test/src/math/lgammabf16_test.cpp | 2 +-
libc/test/src/math/smoke/CMakeLists.txt | 3 +-
libc/test/src/math/smoke/lgammabf16_test.cpp | 52 ++++--
6 files changed, 129 insertions(+), 91 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 792c77bea06ab..4613f94eee304 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -56,7 +56,7 @@ add_header_library(
libc.src.__support.FPUtil.bfloat16
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
- libc.src.__support.FPUtil.PolyEval
+ libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.config
libc.src.__support.macros.optimization
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 70694544fcc33..8cf5b1bfd37ed 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -23,49 +23,11 @@
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)
+// (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,
@@ -90,27 +52,50 @@ LIBC_INLINE_VAR constexpr float LGAMMA_POLY[7][5] = {
0x1.f9d4bep-13f},
};
-// Compute lgamma for x >= 1 using polynomial or Stirling
-LIBC_INLINE float lgamma_positive(float x) {
- if (LIBC_UNLIKELY(x == 1.0f || x == 2.0f))
- return 0.0f;
- // For x >= 8: Stirling approximation
- // lgamma(x) ~ (x-0.5)*ln(x) - x + 0.5*ln(2*pi)
- if (x >= 8.0f) {
- constexpr float HALF_LN_2PI = 0x1.6840e6p-1f;
- float lx = lgamma_logf(x);
- // using FMA: (x - 0.5)*ln(x) - x + 0.5*ln(2*pi)
- return fputil::multiply_add(x - 0.5f, lx, -x + HALF_LN_2PI);
+// 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.
+LIBC_INLINE double lgamma_positive_d(double x) {
+ 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 = __builtin_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 [1, 8): use piecewise polynomial
- int n = static_cast<int>(x); // floor for x >= 1
+ // For x in (0, 4): apply recurrence until reaching [4, 8).
+ // Using the [4,5) polynomial (max_rel_err=6.61e-7) avoids the large errors
+ // near the edges of the [2,3) polynomial (max_rel_err=9.98e-6).
+ double log_product = 0.0;
+ double xs = x;
+ while (xs < 4.0) {
+ log_product += __builtin_log(xs);
+ xs += 1.0;
+ }
+ // xs is now in [4, 8); cast to float for polynomial evaluation.
+ float xf = static_cast<float>(xs);
+ int n = static_cast<int>(xf);
if (n >= 7)
n = 7;
- float t = x - (static_cast<float>(n) + 0.5f);
+ float t = xf - (static_cast<float>(n) + 0.5f);
const float *c = LGAMMA_POLY[n - 1];
- // Co-efficients stored as {c0, c1, c2, c3, c4}
- return fputil::polyeval(t, c[0], c[1], c[2], c[3], c[4]);
+ double lgamma_xs =
+ static_cast<double>(fputil::polyeval(t, c[0], c[1], c[2], c[3], c[4]));
+ return lgamma_xs - log_product;
}
LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
@@ -159,35 +144,55 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
// Negative non-integer: reflection formula
// lgamma(x) = ln(pi) - ln|sin(pi*x)| - lgamma(1-x)
- constexpr float LN_PI = 0x1.250d06p+0f; // ln(pi)
+ constexpr double LN_PI_D = 0x1.250d048e7a1bdp+0;
float ax = -xf;
float frac = ax - static_cast<float>(static_cast<int>(ax));
- constexpr float A1 = 0x1.55553cp-2f; // pi^2/6 approx
- constexpr float A2 = 0x1.11100ap-7f; // pi^4/120 approx
- constexpr float PI = 0x1.921fb6p+1f;
- float frac2 = frac * frac;
- // sin(pi*frac) = PI * frac * polyeval(frac^2, 1, -A1, A2)
- float sin_pi_frac = PI * frac * fputil::polyeval(frac2, 1.0f, -A1, A2);
- if (sin_pi_frac < 0.0f)
- sin_pi_frac = -sin_pi_frac;
-
- float log_sin = lgamma_logf(sin_pi_frac);
- float result = LN_PI - log_sin - lgamma_positive(1.0f + ax);
- return fputil::cast<bfloat16>(result);
- }
-
- // Positive x in (0, 1): use recurrence lgamma(x) = lgamma(x+1) - ln(x)
- if (xf < 1.0f) {
- float lnx = lgamma_logf(xf);
- return fputil::cast<bfloat16>(lgamma_positive(xf + 1.0f) - lnx);
+ // Map frac to [0, 0.5] to guarantee sin is positive
+ if (frac > 0.5f)
+ frac = 1.0f - frac;
+
+ // sin(pi*frac) in double via Taylor series for sin(x)/x:
+ // 1 - x^2/6 + x^4/120 - x^6/5040 + x^8/362880
+ constexpr double PI_D = 0x1.921fb54442d18p+1;
+ double frac_d = static_cast<double>(frac);
+ double x_pi_d = PI_D * frac_d;
+ double x_pi2_d = x_pi_d * x_pi_d;
+ constexpr double DC1 = -0x1.5555555555555p-3; // -1/6
+ constexpr double DC2 = 0x1.1111111111111p-7; // 1/120
+ constexpr double DC3 = -0x1.a01a01a01a01ap-13; // -1/5040
+ constexpr double DC4 = 0x1.71de3a556c734p-19; // 1/362880
+ double sin_pi_frac_d =
+ x_pi_d *
+ (1.0 +
+ x_pi2_d * (DC1 + x_pi2_d * (DC2 + x_pi2_d * (DC3 + x_pi2_d * DC4))));
+
+ double log_sin_d =
+ (sin_pi_frac_d == 1.0) ? 0.0 : __builtin_log(sin_pi_frac_d);
+ // Use double addition: 1.0 + double(ax) preserves tiny ax values that
+ // would be lost by 1.0f + ax in float (e.g. ax=2e-5 rounds to 1.0f).
+ double lgp_d = lgamma_positive_d(1.0 + static_cast<double>(ax));
+ double result_d = LN_PI_D - log_sin_d - lgp_d;
+ // Cast directly from double to bfloat16, bypassing float.
+ // A double->float->bfloat16 chain can double-round: the intermediate
+ // float result may land exactly on a bfloat16 tie point and round the
+ // wrong way, while the original double value was clearly on one side.
+ return fputil::cast<bfloat16>(result_d);
}
- // Positive x >= 1
- return fputil::cast<bfloat16>(lgamma_positive(xf));
+ // Positive x: cast directly from double to bfloat16 to avoid double-rounding.
+ //
+ // Example of the failure without this fix (x = 0x35E5 ≈ 1.706e-6):
+ // lgamma_positive_d returns 13.281250434... (double)
+ // -> static_cast<float> -> 13.28125 exactly (0x41548000)
+ // -> fputil::cast<bfloat16>: bottom 16 bits = 0x8000 (exact tie)
+ // tie-break rounds to even -> 0x4154 = 13.25 (WRONG)
+ // Direct double->bfloat16: 13.281250434 > midpoint 13.28125 -> 0x4155
+ // = 13.3125 (correct)
+ return fputil::cast<bfloat16>(lgamma_positive_d(static_cast<double>(xf)));
}
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
-#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
\ No newline at end of file
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index c03843f2fcfb3..c6d85ed090371 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4869,7 +4869,6 @@ add_entrypoint_object(
libc.src.__support.math.bf16fmal
)
-
add_entrypoint_object(
lgammabf16
SRCS
@@ -4878,7 +4877,6 @@ add_entrypoint_object(
../lgammabf16.h
DEPENDS
libc.src.__support.math.lgammabf16
- libc.src.__support.FPUtil.bfloat16
)
add_entrypoint_object(
diff --git a/libc/test/src/math/lgammabf16_test.cpp b/libc/test/src/math/lgammabf16_test.cpp
index a137257381e84..c35758e6ffd1b 100644
--- a/libc/test/src/math/lgammabf16_test.cpp
+++ b/libc/test/src/math/lgammabf16_test.cpp
@@ -76,4 +76,4 @@ TEST_F(LlvmLibcLgammabf16Test, NormalNegativeRange) {
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
LIBC_NAMESPACE::lgammabf16(x), 0.5);
}
-}
\ No newline at end of file
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index ad7387a52d41c..2d9dbaf48ee4a 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -38,7 +38,6 @@ add_fp_unittest(
add_fp_unittest(
lgammabf16_test
- UNIT_TEST_ONLY
SUITE
libc-math-smoke-tests
SRCS
@@ -46,6 +45,8 @@ add_fp_unittest(
DEPENDS
libc.src.math.lgammabf16
libc.src.__support.FPUtil.bfloat16
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
)
add_fp_unittest(
diff --git a/libc/test/src/math/smoke/lgammabf16_test.cpp b/libc/test/src/math/smoke/lgammabf16_test.cpp
index 393cae08fd10c..b165a0569a600 100644
--- a/libc/test/src/math/smoke/lgammabf16_test.cpp
+++ b/libc/test/src/math/smoke/lgammabf16_test.cpp
@@ -5,6 +5,7 @@
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
+
#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.h"
#include "src/__support/FPUtil/bfloat16.h"
@@ -15,17 +16,50 @@
class LlvmLibcLgammaBf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
DECLARE_SPECIAL_CONSTANTS(bfloat16)
+
public:
void test_special_numbers() {
- EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::lgammabf16(aNaN));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(zero));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(neg_zero));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(inf));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(neg_inf));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(1.0f)));
- EXPECT_FP_EQ(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(2.0f)));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-1.0f)));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-2.0f)));
+ // aNaN -> aNaN, no exception
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::lgammabf16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ // sNaN -> aNaN, FE_INVALID
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::lgammabf16(sNaN),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ // +Inf -> +Inf
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::lgammabf16(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ // -Inf -> +Inf
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::lgammabf16(neg_inf));
+ EXPECT_MATH_ERRNO(0);
+
+ // +-0 -> +Inf, pole error
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(zero), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(neg_zero), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ // lgamma(1) = lgamma(2) = 0
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(1.0f)));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::lgammabf16(bfloat16(2.0f)));
+ EXPECT_MATH_ERRNO(0);
+
+ // Negative integers -> +Inf, pole error
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-1.0f)), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-2.0f)), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
}
};
>From 9f074b0759fb4ff92b07080f45af700f3109e9dc Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Thu, 28 May 2026 11:13:24 +0530
Subject: [PATCH 04/10] added empty line EOF
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/lgammabf16.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 8cf5b1bfd37ed..1e69fab3f379b 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -195,4 +195,4 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
-#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
\ No newline at end of file
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
>From c8f0b1e01ecb19d8f3ca774afde7acb295bb217c Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Thu, 28 May 2026 11:30:23 +0530
Subject: [PATCH 05/10] use math::log instead of unqualified log to avoid
resolving to ::log
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 21 +++++++++++++++++++++
libc/src/__support/math/lgammabf16.h | 8 ++++----
2 files changed, 25 insertions(+), 4 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 4613f94eee304..f492271b9f7a8 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -60,8 +60,10 @@ add_header_library(
libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.config
libc.src.__support.macros.optimization
+ libc.src.__support.math.log
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
)
add_header_library(
@@ -419,6 +421,7 @@ add_header_library(
atanpif16.h
DEPENDS
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
@@ -467,6 +470,7 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.polyeval
@@ -1445,6 +1449,7 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.double_double
libc.src.__support.macros.config
)
@@ -2806,6 +2811,7 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -3066,6 +3072,7 @@ add_header_library(
.expf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.CPP.array
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
@@ -3985,6 +3992,7 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.common
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
@@ -4587,6 +4595,7 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -4601,6 +4610,7 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -4834,6 +4844,7 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
@@ -4907,6 +4918,7 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -4926,6 +4938,7 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -4944,6 +4957,7 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -5022,6 +5036,7 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -5144,6 +5159,7 @@ add_header_library(
.exp_constants
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fenv_impl
@@ -5229,6 +5245,7 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
@@ -5260,6 +5277,7 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.except_value_utils
@@ -5380,6 +5398,7 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -5411,6 +5430,7 @@ add_header_library(
DEPENDS
.expxf16_utils
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.src.__support.CPP.array
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -5450,6 +5470,7 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 1e69fab3f379b..e3eadb0af846f 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -19,6 +19,7 @@
#include "src/__support/FPUtil/multiply_add.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 {
@@ -70,7 +71,7 @@ LIBC_INLINE double lgamma_positive_d(double x) {
if (x >= 8.0) {
// Stirling series; 0.5*ln(2*pi)
constexpr double HALF_LN_2PI = 0x1.d67f1c864beb4p-1;
- double lx = __builtin_log(x);
+ 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);
@@ -83,7 +84,7 @@ LIBC_INLINE double lgamma_positive_d(double x) {
double log_product = 0.0;
double xs = x;
while (xs < 4.0) {
- log_product += __builtin_log(xs);
+ log_product += math::log(xs);
xs += 1.0;
}
// xs is now in [4, 8); cast to float for polynomial evaluation.
@@ -167,8 +168,7 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
(1.0 +
x_pi2_d * (DC1 + x_pi2_d * (DC2 + x_pi2_d * (DC3 + x_pi2_d * DC4))));
- double log_sin_d =
- (sin_pi_frac_d == 1.0) ? 0.0 : __builtin_log(sin_pi_frac_d);
+ double log_sin_d = (sin_pi_frac_d == 1.0) ? 0.0 : math::log(sin_pi_frac_d);
// Use double addition: 1.0 + double(ax) preserves tiny ax values that
// would be lost by 1.0f + ax in float (e.g. ax=2e-5 rounds to 1.0f).
double lgp_d = lgamma_positive_d(1.0 + static_cast<double>(ax));
>From bb8a9c58ea2520f175a06d039f79636b1025f63a Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Thu, 28 May 2026 11:48:16 +0530
Subject: [PATCH 06/10] [libc][math] Fix lgammabf16 deps, remove unused
includes
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 1 -
libc/src/__support/math/lgammabf16.h | 2 --
utils/bazel/llvm-project-overlay/libc/BUILD.bazel | 3 +--
3 files changed, 1 insertion(+), 5 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index f492271b9f7a8..3daecf594f75e 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -63,7 +63,6 @@ add_header_library(
libc.src.__support.math.log
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
)
add_header_library(
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index e3eadb0af846f..50ea78d684f2e 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -10,13 +10,11 @@
#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"
#include "src/__support/math/log.h"
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 55daa7e103e42..d2d1b7e2677c5 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -6237,15 +6237,14 @@ libc_support_library(
name = "__support_math_lgammabf16",
hdrs = ["src/__support/math/lgammabf16.h"],
deps = [
- ":__support_cpp_bit",
":__support_fputil_bfloat16",
":__support_fputil_cast",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
- ":__support_fputil_multiply_add",
":__support_fputil_polyeval",
":__support_macros_config",
":__support_macros_optimization",
+ ":__support_math_log",
":hdr_errno_macros",
":hdr_fenv_macros",
],
>From ddca4aa40e050a329ecc0f2e4e4275e6b34d763c Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Thu, 28 May 2026 11:53:07 +0530
Subject: [PATCH 07/10] removed unused dependencies in cmake
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 2 --
1 file changed, 2 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 3daecf594f75e..61b5207431bc9 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -51,13 +51,11 @@ add_header_library(
HDRS
lgammabf16.h
DEPENDS
- libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.bfloat16
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.polyeval
- libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.config
libc.src.__support.macros.optimization
libc.src.__support.math.log
>From 090c7e33fa93fdb00c3e9fbaa717b8a34b2998e7 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Fri, 29 May 2026 23:14:12 +0530
Subject: [PATCH 08/10] removed cmakelist targets unnecesary codes
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 20 --------------------
1 file changed, 20 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 61b5207431bc9..23c28df9fd427 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -58,7 +58,6 @@ add_header_library(
libc.src.__support.FPUtil.polyeval
libc.src.__support.macros.config
libc.src.__support.macros.optimization
- libc.src.__support.math.log
libc.hdr.errno_macros
libc.hdr.fenv_macros
)
@@ -418,7 +417,6 @@ add_header_library(
atanpif16.h
DEPENDS
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
@@ -467,7 +465,6 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.polyeval
@@ -1446,7 +1443,6 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.double_double
libc.src.__support.macros.config
)
@@ -2808,7 +2804,6 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -3069,7 +3064,6 @@ add_header_library(
.expf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.CPP.array
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
@@ -3989,7 +3983,6 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.common
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
@@ -4592,7 +4585,6 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -4607,7 +4599,6 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -4841,7 +4832,6 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
@@ -4915,7 +4905,6 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -4935,7 +4924,6 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -4954,7 +4942,6 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -5033,7 +5020,6 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
@@ -5156,7 +5142,6 @@ add_header_library(
.exp_constants
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fenv_impl
@@ -5242,7 +5227,6 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
@@ -5274,7 +5258,6 @@ add_header_library(
.expxf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.except_value_utils
@@ -5395,7 +5378,6 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -5427,7 +5409,6 @@ add_header_library(
DEPENDS
.expxf16_utils
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.src.__support.CPP.array
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -5467,7 +5448,6 @@ add_header_library(
.sincosf16_utils
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.math.log
libc.include.llvm-libc-macros.float16_macros
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
>From c803311c0a7ccd4f508dfefacf35d91fe6411a1e Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Tue, 2 Jun 2026 11:57:13 +0530
Subject: [PATCH 09/10] [libc] lgammabf16: move LGAMMA_POLY inside function,
fixed include guard
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/lgammabf16.h | 72 ++++++++++++++--------------
1 file changed, 36 insertions(+), 36 deletions(-)
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 50ea78d684f2e..4441164cc3815 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -5,8 +5,8 @@
// 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
+#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"
@@ -22,35 +22,6 @@
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.
@@ -58,11 +29,40 @@ LIBC_INLINE_VAR constexpr float LGAMMA_POLY[7][5] = {
// 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.
+// +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 fputil::polyeval(t, c[0]..c[4])
+ // where t = x - (n + 0.5), n = 1..7.
+ // Generated by gen_lgamma_bf16_coeffs.py (numpy.polyfit, degree 4).
+ // Maximum relative errors per interval are listed below.
+ static 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;
@@ -193,4 +193,4 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
-#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
+#endif // LLVM_LIBC_SRC__SUPPORT_MATH_LGAMMABF16_H
>From 9c8427f2900e59e5d1a8ac607a4b4d5d1a4a2ad7 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Tue, 2 Jun 2026 12:22:37 +0530
Subject: [PATCH 10/10] fixed issues and comments
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/lgammabf16.h | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 4441164cc3815..2de3253891564 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -5,8 +5,9 @@
// 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
+
+#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"
@@ -37,8 +38,7 @@ 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 fputil::polyeval(t, c[0]..c[4])
// where t = x - (n + 0.5), n = 1..7.
- // Generated by gen_lgamma_bf16_coeffs.py (numpy.polyfit, degree 4).
- // Maximum relative errors per interval are listed below.
+ // Maximum relative errors per interval.
static 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,
@@ -193,4 +193,4 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
} // namespace math
} // namespace LIBC_NAMESPACE_DECL
-#endif // LLVM_LIBC_SRC__SUPPORT_MATH_LGAMMABF16_H
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
More information about the libc-commits
mailing list