[libc-commits] [libc] [llvm] [libc][math] Impl bfloat16 lgamma function. (PR #199312)
Uday Kiriti via libc-commits
libc-commits at lists.llvm.org
Fri Jul 24 02:43:17 PDT 2026
https://github.com/udaykiriti updated https://github.com/llvm/llvm-project/pull/199312
>From 05637bd73a805fd51c711b98a81ae5bf97d4abe3 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/15] [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 93808e4c471fc..7f0db7091070d 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/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/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 71ecbe438bf7b..b7fc7f2592aaf 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -857,6 +857,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 b38f60ab9d5f5..45c95a0d0a6e2 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -854,6 +854,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 15adad72ab459..42f901a5555a6 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -666,6 +666,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 18184e8136a9e..0a19c9b7b5a9c 100644
--- a/libc/config/darwin/x86_64/entrypoints.txt
+++ b/libc/config/darwin/x86_64/entrypoints.txt
@@ -280,6 +280,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 f94a6edd69759..2f48ebdd3e9a1 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -693,6 +693,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 3c52dd8437daa..7ba94c9bd6134 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -709,6 +709,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 2e7b9276ab068..58979d49deb7b 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -979,6 +979,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 1986d6a5347dc..43cf36f10e7dc 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -553,6 +553,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 8a29e77ffe0fd..701fbb9b4ab60 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -1071,6 +1071,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 0a4fb747c2940..aab3cca9061ff 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -1076,6 +1076,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 5333bc4042070..6e4e2225e8299 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -370,6 +370,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 b86fcbc0abcea..ba389a1c2ce0c 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -310,6 +310,7 @@
#include "math/ldexpf128.h"
#include "math/ldexpf16.h"
#include "math/ldexpl.h"
+#include "math/lgammabf16.h"
#include "math/lgammaf16.h"
#include "math/llogb.h"
#include "math/llogbbf16.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 9f4624682eafe..858ae4fae496e 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -61,6 +61,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 8c9a184ca517f..e59439afbbea7 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -604,6 +604,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 d5661f6f55e91..84f49390a46e6 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4942,6 +4942,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 9208dbcc5ade6..d4059f90406b0 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -312,6 +312,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 9da54744b5307..268464703bc78 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -791,7 +791,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 240689c0de02f..adf64dd7f1565 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 81f76fa7681fd..fc6aa3c861724 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 2422bcf45222f..61150c06f3a5c 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -386,6 +386,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 38fd15fcc956c..8a674ef46ff3e 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -219,6 +219,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 d585baa2e0d2c..86c240b8fd4a4 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -89,6 +89,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 84e59674295f5..51d078b689e16 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -55,6 +55,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 bee6f7aa5766c..6233ba0d42548 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7094,6 +7094,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"],
@@ -12369,6 +12385,13 @@ libc_math_function(
additional_deps = [":__support_math_lgammaf16"],
)
+libc_math_function(
+ name = "lgammabf16",
+ additional_deps = [
+ ":__support_math_lgammabf16",
+ ],
+)
+
libc_math_function(
name = "llogb",
additional_deps = [
>From 2804eb9cbb81afab3f6b9d4acde747996f93d95b 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/15] [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 858ae4fae496e..6461c1b1fc96b 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -71,6 +71,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 6233ba0d42548..f1697a1e84432 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7103,6 +7103,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 8712262b0a271530b568e2658525a56ce6dad99f 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/15] [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 6461c1b1fc96b..2af11733b54ea 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -71,7 +71,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 84f49390a46e6..31930682eded1 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4942,7 +4942,6 @@ add_entrypoint_object(
libc.src.__support.math.bf16fmal
)
-
add_entrypoint_object(
lgammabf16
SRCS
@@ -4951,7 +4950,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 fc6aa3c861724..996c2d7b419a3 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 ab1a9a5bb527c268ba5b351e2bc8437086f1faa0 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/15] 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 3ef41d1d0821bf20d74b4a08331ed830ca091e0f 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/15] 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 2af11733b54ea..06f050f602ce3 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -75,8 +75,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(
@@ -434,6 +436,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
@@ -482,6 +485,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
@@ -1474,6 +1478,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
)
@@ -2855,6 +2860,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
@@ -3115,6 +3121,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
@@ -4065,6 +4072,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
@@ -4673,6 +4681,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
@@ -4682,6 +4691,7 @@ add_header_library(
libc.src.__support.macros.properties.types
)
+
add_header_library(
rsqrtf16
HDRS
@@ -4926,6 +4936,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
@@ -4999,6 +5010,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
@@ -5018,6 +5030,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
@@ -5036,6 +5049,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
@@ -5114,6 +5128,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
@@ -5236,6 +5251,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
@@ -5321,6 +5337,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
@@ -5352,6 +5369,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
@@ -5489,6 +5507,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
@@ -5520,6 +5539,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
@@ -5559,6 +5579,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 167d7feb6087c1129e6f0a77fa01b3efe1012d68 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/15] [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 06f050f602ce3..5acf34942ec7b 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -78,7 +78,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 f1697a1e84432..e2198d2b5b1fd 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7098,15 +7098,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 14b1fdc6d959002bebe7298e175017f20b753aa9 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/15] 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 5acf34942ec7b..ad553c4d7ca9d 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -66,13 +66,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 99eca46d727a26f7f3a629c1a42a32e048c966d8 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/15] removed cmakelist targets unnecesary codes
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 22 ----------------------
1 file changed, 22 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index ad553c4d7ca9d..ef1a577cb1983 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -73,7 +73,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
)
@@ -433,7 +432,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
@@ -482,7 +480,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
@@ -1475,7 +1472,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
)
@@ -2857,7 +2853,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
@@ -3118,7 +3113,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
@@ -4069,7 +4063,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
@@ -4678,7 +4671,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
@@ -4688,7 +4680,6 @@ add_header_library(
libc.src.__support.macros.properties.types
)
-
add_header_library(
rsqrtf16
HDRS
@@ -4696,8 +4687,6 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.include.llvm-libc-macros.float16_macros
- libc.src.__support.CPP.bit
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
@@ -4933,7 +4922,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
@@ -5007,7 +4995,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
@@ -5027,7 +5014,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
@@ -5046,7 +5032,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
@@ -5125,7 +5110,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
@@ -5248,7 +5232,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
@@ -5334,7 +5317,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
@@ -5366,7 +5348,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
@@ -5504,7 +5485,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
@@ -5536,7 +5516,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
@@ -5576,7 +5555,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 29e4d8787346906db14df166100ccdae927033a3 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/15] [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 f278698d09ba66c4e329fb4c7ec7e61564c913c2 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/15] 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
>From c190a41f5bb5015be4bb42f4636db8771426d896 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Sat, 4 Jul 2026 21:16:57 +0530
Subject: [PATCH 11/15] [libc][math] lgammabf16
- loop unrolled
- used estrin's scheme..
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/lgammabf16.h | 105 +++++++++++++++++++--------
1 file changed, 76 insertions(+), 29 deletions(-)
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 2de3253891564..0b6219f670ca6 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -16,6 +16,8 @@
#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"
@@ -26,17 +28,22 @@ namespace math {
// 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 < 8, applies the recurrence lgamma(x) = lgamma(x+1) - ln(x) until
+// 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 fputil::polyeval(t, c[0]..c[4])
+ // Each row: {c0, c1, c2, c3, c4} for Estrin evaluation
// where t = x - (n + 0.5), n = 1..7.
// Maximum relative errors per interval.
static constexpr float LGAMMA_POLY[7][5] = {
@@ -76,24 +83,49 @@ LIBC_INLINE double lgamma_positive_d(double x) {
return result;
}
- // 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 += math::log(xs);
- xs += 1.0;
+ // For x in (0, 4): apply recurrence lgamma(x) = lgamma(x+n) -
+ // ln(x*(x+1)*...*(x+n-1)) until xs reaches [4, 8). We unroll the loop
+ // explicitly for efficiency. 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, xs;
+
+ if (x >= 3.0) {
+ log_product = math::log(x);
+ xs = x + 1.0;
+ } else if (x >= 2.0) {
+ log_product = math::log(x * (x + 1.0));
+ xs = x + 2.0;
+ } else if (x >= 1.0) {
+ log_product = math::log(x * (x + 1.0) * (x + 2.0));
+ xs = x + 3.0;
+ } else {
+ log_product = math::log(x * (x + 1.0) * (x + 2.0) * (x + 3.0));
+ xs = x + 4.0;
}
- // xs is now in [4, 8); cast to float for polynomial evaluation.
+
+ // 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);
- const float *c = LGAMMA_POLY[n - 1];
- double lgamma_xs =
- static_cast<double>(fputil::polyeval(t, c[0], c[1], c[2], c[3], c[4]));
+
+ 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;
}
@@ -101,7 +133,7 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
using FPBits = fputil::FPBits<bfloat16>;
FPBits x_bits(x);
- // Handle NaN
+ // Handles NaN
if (LIBC_UNLIKELY(x_bits.is_nan())) {
if (x_bits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
@@ -145,32 +177,47 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
// lgamma(x) = ln(pi) - ln|sin(pi*x)| - lgamma(1-x)
constexpr double LN_PI_D = 0x1.250d048e7a1bdp+0;
float ax = -xf;
- float frac = ax - static_cast<float>(static_cast<int>(ax));
- // Map frac to [0, 0.5] to guarantee sin is positive
+ // nearest_integer avoids truncation-toward-zero of static_cast<int>
+ float frac = fputil::abs(ax - fputil::nearest_integer(ax));
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
+ // sin(pi*frac) via degree-4 Taylor series in double.
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 : math::log(sin_pi_frac_d);
+
+ // DC[k] = (-1)^(k+1) / (2k+2)! for k=0..3
+ constexpr double DC[4] = {
+ -0x1.5555555555555p-3, // -1/6 = -1/3!
+ 0x1.1111111111111p-7, // 1/120 = 1/5!
+ -0x1.a01a01a01a01ap-13, // -1/5040 = -1/7!
+ 0x1.71de3a556c734p-19, // 1/362880 = 1/9!
+ };
+
+ // Estrin's scheme for q(u) = 1 + DC[0]*u + DC[1]*u^2 + DC[2]*u^3 +
+ // DC[3]*u^4 where u = x_pi2_d:
+ // q(u) = (1 + DC[0]*u) + u^2 * ((DC[1] + DC[2]*u) + u^2 * DC[3])
+ double u2 = x_pi2_d * x_pi2_d;
+ double q01 = fputil::multiply_add(x_pi2_d, DC[0], 1.0);
+ double q12 = fputil::multiply_add(x_pi2_d, DC[2], DC[1]);
+ double q123 = fputil::multiply_add(u2, DC[3], q12);
+ double poly = fputil::multiply_add(u2, q123, q01);
+ double sin_pi_frac_d = x_pi_d * poly;
+
+ // A fast (not correctly-rounded) log suffices here: the final result is
+ // cast to bfloat16, so we only need ~8 bits of accuracy in log_sin_d.
+ double log_sin_d = (sin_pi_frac_d == 1.0)
+ ? 0.0
+ : math::log(static_cast<float>(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
>From 5771a43eb691f8773064d524a68c0d1690d82d1b Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Sat, 11 Jul 2026 18:35:50 +0530
Subject: [PATCH 12/15] [libc][math] bfloat16 lgamma
- added missing dependnecies...
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 5 +-
libc/src/__support/math/lgammabf16.h | 46 +++++++++++++------
.../llvm-project-overlay/libc/BUILD.bazel | 2 +
3 files changed, 38 insertions(+), 15 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index ef1a577cb1983..877e9bf0d735c 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -70,9 +70,12 @@ 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.FPUtil.nearest_integer
libc.src.__support.macros.config
libc.src.__support.macros.optimization
+ libc.src.__support.math.log
libc.hdr.errno_macros
libc.hdr.fenv_macros
)
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index 0b6219f670ca6..71e8494aa1527 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -25,6 +25,8 @@
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.
@@ -40,13 +42,21 @@ namespace math {
// 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.
- static constexpr float LGAMMA_POLY[7][5] = {
+ // 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},
@@ -76,31 +86,35 @@ 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 = math::log(x);
+ double lx = math::log_internal::log_finite(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 lgamma(x) = lgamma(x+n) -
- // ln(x*(x+1)*...*(x+n-1)) until xs reaches [4, 8). We unroll the loop
- // explicitly for efficiency. 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, xs;
+ // 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) {
- log_product = math::log(x * (x + 1.0));
+ product = x * (x + 1.0);
+ log_product = math::log(product);
xs = x + 2.0;
} else if (x >= 1.0) {
- log_product = math::log(x * (x + 1.0) * (x + 2.0));
+ product = x * (x + 1.0);
+ product = product * (x + 2.0);
+ log_product = math::log(product);
xs = x + 3.0;
} else {
- log_product = math::log(x * (x + 1.0) * (x + 2.0) * (x + 3.0));
+ product = x * (x + 1.0);
+ product = product * (x + 2.0);
+ product = product * (x + 3.0);
+ log_product = math::log(product);
xs = x + 4.0;
}
@@ -129,6 +143,8 @@ LIBC_INLINE double lgamma_positive_d(double x) {
return lgamma_xs - log_product;
}
+} // namespace lgammabf16_internal
+
LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
using FPBits = fputil::FPBits<bfloat16>;
FPBits x_bits(x);
@@ -215,7 +231,8 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
// 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 lgp_d =
+ lgammabf16_internal::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.
@@ -234,7 +251,8 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
// 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)));
+ return fputil::cast<bfloat16>(
+ lgammabf16_internal::lgamma_positive_d(static_cast<double>(xf)));
}
} // namespace math
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index e2198d2b5b1fd..d562f45e34474 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7102,6 +7102,8 @@ libc_support_library(
":__support_fputil_cast",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer",
":__support_fputil_polyeval",
":__support_macros_config",
":__support_macros_optimization",
>From 4538155c75cab4b19471a5c74464575b20c12b3c Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Sat, 11 Jul 2026 18:41:09 +0530
Subject: [PATCH 13/15] ....
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 71e8494aa1527..d9804f745dfda 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -86,7 +86,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 = math::log_internal::log_finite(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);
>From 92f0f9bfdf949d4e6b19948d12e93661192187e9 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Fri, 24 Jul 2026 15:04:08 +0530
Subject: [PATCH 14/15] updated comments for consts
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
libc/src/__support/math/CMakeLists.txt | 1 -
libc/src/__support/math/lgammabf16.h | 93 ++++++++++----------
libc/test/src/math/smoke/lgammabf16_test.cpp | 40 +++++++++
libc/utils/MPFRWrapper/MPCommon.cpp | 7 --
libc/utils/MPFRWrapper/MPCommon.h | 1 -
libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 -
libc/utils/MPFRWrapper/MPFRUtils.h | 1 -
7 files changed, 86 insertions(+), 59 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 877e9bf0d735c..e6dbc89cfd2ae 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -70,7 +70,6 @@ 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.FPUtil.nearest_integer
libc.src.__support.macros.config
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
index d9804f745dfda..a4fb097e1c945 100644
--- a/libc/src/__support/math/lgammabf16.h
+++ b/libc/src/__support/math/lgammabf16.h
@@ -1,4 +1,4 @@
-//===-- Implementation of lgammabf16 ----------------------------*- C++ -*-===//
+//===-- 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.
@@ -13,7 +13,6 @@
#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"
@@ -24,7 +23,6 @@
namespace LIBC_NAMESPACE_DECL {
namespace math {
-
namespace lgammabf16_internal {
// lgamma_positive_d: compute lgamma(x) for x > 0, returning double.
@@ -46,38 +44,38 @@ 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):
+ // P is a degree-4 fit to lgamma(t + n+0.5) on t in [-0.5, 0.5].
+ // For n = 3..7: mpmath.chebyfit (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)
+ // For n = 1, 2: chebyfit's relative error is too poor near the
+ // zeros of lgamma at x=1 and x=2 (see fit_irls_relative_error()),
+ // so these two rows instead use an IRLS (Lawson-iteration) fit that
+ // directly minimizes relative error rather than plain interpolation.
// 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},
+ // [1,2), center=1.5, max_rel_err=1.95e-03 (IRLS relative-error fit)
+ {-0x1.ee9a58p-4f, 0x1.36b582p-5f, 0x1.dd697ep-2f, -0x1.36ad56p-3f,
+ 0x1.130d8p-4f},
+ // [2,3), center=2.5, max_rel_err=4.32e-05 (IRLS relative-error fit)
+ {0x1.2386f2p-2f, 0x1.6803b6p-1f, 0x1.f53b64p-3f, -0x1.46e712p-5f,
+ 0x1.736f46p-7f},
+ // [3,4), center=3.5, max_rel_err=1.79e-06 (chebyfit)
+ {0x1.337302p+0f, 0x1.1a6936p+0f, 0x1.52480cp-3f, -0x1.2a67eep-6f,
+ 0x1.84fbbep-9f},
+ // [4,5), center=4.5, max_rel_err=2.22e-07 (chebyfit)
+ {0x1.3a140ap+1f, 0x1.638d48p+0f, 0x1.fd62f6p-4f, -0x1.52194ep-7f,
+ 0x1.4db306p-10f},
+ // [5,6), center=5.5, max_rel_err=5.49e-08 (chebyfit)
+ {0x1.fa99a6p+1f, 0x1.9c70b4p+0f, 0x1.98409cp-4f, -0x1.b23a9cp-8f,
+ 0x1.5870b8p-11f},
+ // [6,7), center=6.5, max_rel_err=5.71e-08 (chebyfit)
+ {0x1.6a676ap+2f, 0x1.cafc4ap+0f, 0x1.548ce6p-4f, -0x1.2e1bc6p-8f,
+ 0x1.906b32p-12f},
+ // [7,8), center=7.5, max_rel_err=1.04e-08 (chebyfit)
+ {0x1.e23306p+2f, 0x1.f25ebap+0f, 0x1.2413c4p-4f, -0x1.bc6a74p-9f,
+ 0x1.f9ac9ep-13f},
};
if (LIBC_UNLIKELY(x == 1.0 || x == 2.0))
@@ -85,7 +83,8 @@ 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;
+ // > mpmath.mp.dps = 50; hex(0.5 * mpmath.log(2 * mpmath.pi))
+ constexpr double HALF_LN_2PI = 0x1.d67f1c864beb5p-1;
double lx = math::log(x);
double x2 = x * x;
double result = (x - 0.5) * lx - x + HALF_LN_2PI;
@@ -96,9 +95,13 @@ LIBC_INLINE double lgamma_positive_d(double x) {
// 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.
+ // x already in [4, 8) needs no shift -- evaluate the polynomial directly.
double log_product, xs, product;
- if (x >= 3.0) {
+ if (x >= 4.0) {
+ log_product = 0.0;
+ xs = x;
+ } else if (x >= 3.0) {
log_product = math::log(x);
xs = x + 1.0;
} else if (x >= 2.0) {
@@ -142,7 +145,6 @@ LIBC_INLINE double lgamma_positive_d(double x) {
return lgamma_xs - log_product;
}
-
} // namespace lgammabf16_internal
LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
@@ -191,6 +193,7 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
// Negative non-integer: reflection formula
// lgamma(x) = ln(pi) - ln|sin(pi*x)| - lgamma(1-x)
+ // > mpmath.mp.dps = 50; hex(mpmath.log(mpmath.pi))
constexpr double LN_PI_D = 0x1.250d048e7a1bdp+0;
float ax = -xf;
@@ -200,12 +203,18 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
frac = 1.0f - frac;
// sin(pi*frac) via degree-4 Taylor series in double.
+ // > mpmath.mp.dps = 50; hex(mpmath.pi)
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;
- // DC[k] = (-1)^(k+1) / (2k+2)! for k=0..3
+ // Taylor coefficients for sin(y)/y = 1 - y^2/3! + y^4/5! - y^6/7! + y^8/9!,
+ // evaluated at u = y^2 (so q(u) = 1 + DC[0]*u + DC[1]*u^2 + DC[2]*u^3 +
+ // DC[3]*u^4). DC[k] = (-1)^(k+1) / (2k+3)! for k = 0..3, generated via:
+ // > mpmath.mp.dps = 50
+ // > [hex((-1)**(k+1) * mpmath.mpf(1) / mpmath.factorial(2*k+3)) for k in
+ // range(4)]
constexpr double DC[4] = {
-0x1.5555555555555p-3, // -1/6 = -1/3!
0x1.1111111111111p-7, // 1/120 = 1/5!
@@ -235,22 +244,12 @@ LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
lgammabf16_internal::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.
+ // Cast directly from double->bfloat16 (not via float) to avoid
+ // double-rounding at bfloat16 tie points.
return fputil::cast<bfloat16>(result_d);
}
- // 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)
+ // Cast directly from double->bfloat16 for the same reason as above.
return fputil::cast<bfloat16>(
lgammabf16_internal::lgamma_positive_d(static_cast<double>(xf)));
}
diff --git a/libc/test/src/math/smoke/lgammabf16_test.cpp b/libc/test/src/math/smoke/lgammabf16_test.cpp
index b165a0569a600..193f66eafa197 100644
--- a/libc/test/src/math/smoke/lgammabf16_test.cpp
+++ b/libc/test/src/math/smoke/lgammabf16_test.cpp
@@ -61,6 +61,46 @@ class LlvmLibcLgammaBf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-2.0f)), FE_DIVBYZERO);
EXPECT_MATH_ERRNO(ERANGE);
}
+ void test_negative_integers() {
+ // More negative integer poles -> +Inf, pole error
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-3.0f)), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-4.0f)), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-10.0f)), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::lgammabf16(bfloat16(-100.0f)), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
+
+ void test_exact_values() {
+ // lgamma(3) = ln(2) ~ 0.693147
+ EXPECT_FP_EQ(bfloat16(0x1.62p-1f),
+ LIBC_NAMESPACE::lgammabf16(bfloat16(3.0f)));
+ EXPECT_MATH_ERRNO(0);
+ // lgamma(4) = ln(6) ~ 1.791759
+ EXPECT_FP_EQ(bfloat16(0x1.cap+0f),
+ LIBC_NAMESPACE::lgammabf16(bfloat16(4.0f)));
+ EXPECT_MATH_ERRNO(0);
+ // lgamma(0.5) = ln(sqrt(pi)) ~ 0.572365
+ EXPECT_FP_EQ(bfloat16(0x1.26p-1f),
+ LIBC_NAMESPACE::lgammabf16(bfloat16(0.5f)));
+ EXPECT_MATH_ERRNO(0);
+ }
+
+ void test_overflow() {
+ // Large x -> lgamma overflows to +Inf
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammabf16(max_normal),
+ FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
};
TEST_F(LlvmLibcLgammaBf16Test, SpecialNumbers) { test_special_numbers(); }
+TEST_F(LlvmLibcLgammaBf16Test, NegativeIntegers) { test_negative_integers(); }
+TEST_F(LlvmLibcLgammaBf16Test, ExactValues) { test_exact_values(); }
+TEST_F(LlvmLibcLgammaBf16Test, Overflow) { test_overflow(); }
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index 61150c06f3a5c..2422bcf45222f 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -386,13 +386,6 @@ 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 8a674ef46ff3e..38fd15fcc956c 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -219,7 +219,6 @@ 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 86c240b8fd4a4..d585baa2e0d2c 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -89,8 +89,6 @@ 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 51d078b689e16..84e59674295f5 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -55,7 +55,6 @@ enum class Operation : int {
Log10,
Log10p1,
Log1p,
- Lgamma,
Mod2PI,
ModPIOver2,
ModPIOver4,
>From 6a3f724913763f848374bf47b7853c4afd48b483 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Fri, 24 Jul 2026 15:12:48 +0530
Subject: [PATCH 15/15] removed unused dependency from bazel
Signed-off-by: udaykiriti <udaykiriti9 at gmail.com>
---
utils/bazel/llvm-project-overlay/libc/BUILD.bazel | 1 -
1 file changed, 1 deletion(-)
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index d562f45e34474..5ece4758c41a2 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7104,7 +7104,6 @@ libc_support_library(
":__support_fputil_fp_bits",
":__support_fputil_multiply_add",
":__support_fputil_nearest_integer",
- ":__support_fputil_polyeval",
":__support_macros_config",
":__support_macros_optimization",
":__support_math_log",
More information about the libc-commits
mailing list