[libc-commits] [libc] [llvm] [libc][math] Add tgammabf16 math support for bfloat16. (PR #208689)
Vedant Neve via libc-commits
libc-commits at lists.llvm.org
Thu Jul 16 03:28:19 PDT 2026
https://github.com/0bVdnt updated https://github.com/llvm/llvm-project/pull/208689
>From 7d19fcec11b40bddb1a203f2c02f7011d58bb1bf Mon Sep 17 00:00:00 2001
From: Vedant Neve <vedantneve13 at gmail.com>
Date: Tue, 7 Jul 2026 11:42:52 +0530
Subject: [PATCH 1/4] [libc][math] Add bfloat16 tgamma function
This patch implements the `tgammabf16` function for LLVM-libc.
Fixes #207495
Signed-off-by: Vedant Neve <vedantneve13 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/freebsd/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 | 2 +-
libc/shared/math.h | 1 +
libc/shared/math/tgammabf16.h | 23 +++
libc/src/__support/math/CMakeLists.txt | 19 ++
libc/src/__support/math/tgammabf16.h | 179 ++++++++++++++++++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/generic/CMakeLists.txt | 10 +
libc/src/math/generic/tgammabf16.cpp | 18 ++
libc/src/math/tgammabf16.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/smoke/CMakeLists.txt | 13 ++
libc/test/src/math/smoke/tgammabf16_test.cpp | 79 ++++++++
libc/test/src/math/tgammabf16_test.cpp | 42 ++++
libc/utils/MPFRWrapper/MPCommon.cpp | 6 +
libc/utils/MPFRWrapper/MPCommon.h | 1 +
libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 +
libc/utils/MPFRWrapper/MPFRUtils.h | 1 +
.../llvm-project-overlay/libc/BUILD.bazel | 26 +++
32 files changed, 471 insertions(+), 1 deletion(-)
create mode 100644 libc/shared/math/tgammabf16.h
create mode 100644 libc/src/__support/math/tgammabf16.h
create mode 100644 libc/src/math/generic/tgammabf16.cpp
create mode 100644 libc/src/math/tgammabf16.h
create mode 100644 libc/test/src/math/smoke/tgammabf16_test.cpp
create mode 100644 libc/test/src/math/tgammabf16_test.cpp
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index 149c790bb2448..5db00c527ef52 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -869,6 +869,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 8a8ab3d6a78fd..e9e34e59d3688 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -880,6 +880,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index d0b8f091cd6e2..b22952210b0d7 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -877,6 +877,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index eeb9ec244b415..1b6174a569d2d 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -690,6 +690,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/darwin/x86_64/entrypoints.txt b/libc/config/darwin/x86_64/entrypoints.txt
index 18184e8136a9e..e561cdae925f6 100644
--- a/libc/config/darwin/x86_64/entrypoints.txt
+++ b/libc/config/darwin/x86_64/entrypoints.txt
@@ -306,6 +306,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/freebsd/x86_64/entrypoints.txt b/libc/config/freebsd/x86_64/entrypoints.txt
index df6c0f907725b..8ea696469cc3c 100644
--- a/libc/config/freebsd/x86_64/entrypoints.txt
+++ b/libc/config/freebsd/x86_64/entrypoints.txt
@@ -629,6 +629,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index ea5dff846499f..17e797e9c7ae2 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -713,6 +713,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index 8ccb39cfd6605..4bff935042427 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -715,6 +715,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index c4db45ad3c50e..12200581643ca 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -991,6 +991,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 20763eec40e8b..556caef00465d 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -568,6 +568,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index bef3bca77df0e..0613f407ec0de 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -1012,6 +1012,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index e0bf4ff52d13a..fe902b29db717 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -1083,6 +1083,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 5333bc4042070..9a71f82a5dc22 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -395,6 +395,7 @@ list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.setpayloadsigbf16
libc.src.math.sqrtbf16
libc.src.math.tanbf16
+ libc.src.math.tgammabf16
libc.src.math.truncbf16
libc.src.math.totalorderbf16
libc.src.math.totalordermagbf16
@@ -406,4 +407,3 @@ set(TARGET_LLVMLIBC_ENTRYPOINTS
${TARGET_LIBC_ENTRYPOINTS}
${TARGET_LIBM_ENTRYPOINTS}
)
-
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 8e01154c5ab56..edae1d5f77c55 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -482,6 +482,7 @@
#include "math/tanhf16.h"
#include "math/tanpif.h"
#include "math/tanpif16.h"
+#include "math/tgammabf16.h"
#include "math/totalorder.h"
#include "math/totalorderbf16.h"
#include "math/totalorderf.h"
diff --git a/libc/shared/math/tgammabf16.h b/libc/shared/math/tgammabf16.h
new file mode 100644
index 0000000000000..1d1db83182300
--- /dev/null
+++ b/libc/shared/math/tgammabf16.h
@@ -0,0 +1,23 @@
+//===-- Shared tgammabf16 function ------------------------------*- 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_TGAMMABF16_H
+#define LLVM_LIBC_SHARED_MATH_TGAMMABF16_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/tgammabf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::tgammabf16;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_TGAMMABF16_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index deafc85b487a6..a076a4f219992 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -5438,6 +5438,25 @@ add_header_library(
libc.src.__support.FPUtil.bfloat16
)
+add_header_library(
+ tgammabf16
+ HDRS
+ tgammabf16.h
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.bfloat16
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.math.exp
+ libc.src.__support.math.log
+ libc.src.__support.math.sin
+)
+
add_header_library(
tanf
HDRS
diff --git a/libc/src/__support/math/tgammabf16.h b/libc/src/__support/math/tgammabf16.h
new file mode 100644
index 0000000000000..275f9fc15745b
--- /dev/null
+++ b/libc/src/__support/math/tgammabf16.h
@@ -0,0 +1,179 @@
+//===-- Implementation header for tgammabf16 --------------------*- 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_TGAMMABF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TGAMMABF16_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/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/exp.h"
+#include "src/__support/math/log.h"
+#include "src/__support/math/sin.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+namespace tgammabf16_internal {
+
+LIBC_INLINE_VAR constexpr double PI = 0x1.921fb54442d18p+1;
+LIBC_INLINE_VAR constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
+
+// Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
+LIBC_INLINE_VAR constexpr double LANCZOS_COEFFS[9] = {
+ 0.99999999999980993, 676.5203681218851, -1259.1392167224028,
+ 771.32342877765313, -176.61502916214059, 12.507343278224757,
+ -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
+
+LIBC_INLINE bool is_negative_integer(double x) {
+ if (x > -1.0)
+ return false;
+ if (x <= -128.0)
+ return true;
+ int n = static_cast<int>(x);
+ return x == static_cast<double>(n);
+}
+
+} // namespace tgammabf16_internal
+
+LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
+ using FPBits = fputil::FPBits<bfloat16>;
+ using namespace tgammabf16_internal;
+
+ FPBits xbits(x);
+
+ if (LIBC_UNLIKELY(xbits.is_nan())) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ if (LIBC_UNLIKELY(xbits.is_inf())) {
+ if (xbits.is_pos())
+ return x;
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ if (LIBC_UNLIKELY(xbits.is_zero())) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return FPBits::inf(xbits.sign()).get_val();
+ }
+
+ double xd = static_cast<double>(static_cast<float>(x));
+ if (LIBC_UNLIKELY(is_negative_integer(xd))) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // Fast path for exact positive integers
+ if (xd > 0.0 && xd <= 35.0) {
+ int n = static_cast<int>(xd);
+ if (xd == static_cast<double>(n)) {
+ double res = 1.0;
+ for (int i = 1; i < n; ++i)
+ res *= i;
+ return fputil::cast<bfloat16>(res);
+ }
+ }
+
+ bool reflection = false;
+ bool divide_by_x = false;
+ double x_eval = xd;
+ double res;
+
+ // Fast path for tiny positive inputs to prevent exact-boundary overshoots
+ // For tiny x, Gamma(x) ~= 1/x - gamma
+ if (LIBC_UNLIKELY(xd > 0.0 && xd < 0x1.0p-8)) {
+ res = (1.0 / xd) - 0.577215664901532860606;
+ } else {
+ if (xd < 0.0) {
+ reflection = true;
+ x_eval = 1.0 - xd;
+ } else if (xd < 1.0) {
+ divide_by_x = true;
+ x_eval = xd + 1.0;
+ }
+
+ double z = x_eval - 1.0;
+ double a = LANCZOS_COEFFS[0];
+ for (int i = 1; i < 9; ++i) {
+ a += LANCZOS_COEFFS[i] / (z + static_cast<double>(i));
+ }
+
+ double t = z + 7.5;
+
+ // ILP Optimization: FMA for the exponent argument
+ double log_t = math::log(t);
+ double base_term = (LOG_SQRT_2_PI - t) + math::log(a);
+ double exp_arg = fputil::multiply_add(z + 0.5, log_t, base_term);
+
+ res = math::exp(exp_arg);
+
+ if (reflection) {
+ double sin_pi_x = math::sin(PI * xd);
+ double denom = sin_pi_x * res;
+
+ if (LIBC_UNLIKELY(denom == 0.0)) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_UNDERFLOW);
+ return FPBits::zero(Sign::POS).get_val();
+ }
+ res = PI / denom;
+ } else if (divide_by_x) {
+ res = res / xd;
+ }
+ }
+
+ double abs_res = res < 0.0 ? -res : res;
+
+ // 0x1.ffp127 is max normal + 0.5 ULP in bfloat16, the exact RTN overflow
+ // threshold
+ if (LIBC_UNLIKELY(abs_res > 0x1.ffp+127)) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return FPBits::inf(res > 0.0 ? Sign::POS : Sign::NEG).get_val();
+#else
+ Sign sign = res > 0.0 ? Sign::POS : Sign::NEG;
+ switch (fputil::quick_get_round()) {
+ case FE_TONEAREST:
+ return FPBits::inf(sign).get_val();
+ case FE_UPWARD:
+ return sign == Sign::POS ? FPBits::inf(Sign::POS).get_val()
+ : FPBits::max_normal(Sign::NEG).get_val();
+ case FE_DOWNWARD:
+ return sign == Sign::POS ? FPBits::max_normal(Sign::POS).get_val()
+ : FPBits::inf(Sign::NEG).get_val();
+ case FE_TOWARDZERO:
+ return FPBits::max_normal(sign).get_val();
+ default:
+ return FPBits::max_normal(sign).get_val();
+ }
+#endif
+ }
+
+ return fputil::cast<bfloat16>(res);
+}
+
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TGAMMABF16_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index fd768eaed26af..f3e86d9439f4a 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -600,6 +600,7 @@ add_math_entrypoint_object(tanpif)
add_math_entrypoint_object(tanpif16)
add_math_entrypoint_object(tgamma)
+add_math_entrypoint_object(tgammabf16)
add_math_entrypoint_object(tgammaf)
add_math_entrypoint_object(lgamma)
add_math_entrypoint_object(lgamma_r)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 28b708e0c8f7a..64da6c553eb25 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -425,6 +425,16 @@ add_entrypoint_object(
libc.src.__support.math.tanbf16
)
+add_entrypoint_object(
+ tgammabf16
+ SRCS
+ tgammabf16.cpp
+ HDRS
+ ../tgammabf16.h
+ DEPENDS
+ libc.src.__support.math.tgammabf16
+)
+
add_entrypoint_object(
tanf
SRCS
diff --git a/libc/src/math/generic/tgammabf16.cpp b/libc/src/math/generic/tgammabf16.cpp
new file mode 100644
index 0000000000000..c5772f00c154f
--- /dev/null
+++ b/libc/src/math/generic/tgammabf16.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of tgammabf16 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/__support/math/tgammabf16.h"
+#include "src/math/tgammabf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(bfloat16, tgammabf16, (bfloat16 x)) {
+ return math::tgammabf16(x);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/tgammabf16.h b/libc/src/math/tgammabf16.h
new file mode 100644
index 0000000000000..1aa93bb10c103
--- /dev/null
+++ b/libc/src/math/tgammabf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for tgammabf16 --------------------*- 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_TGAMMABF16_H
+#define LLVM_LIBC_SRC_MATH_TGAMMABF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+bfloat16 tgammabf16(bfloat16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_TGAMMABF16_H
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index 5156e7e5b32d0..1685a9234c055 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -487,6 +487,7 @@ add_fp_unittest(
libc.src.__support.math.tanhf16
libc.src.__support.math.tanpif
libc.src.__support.math.tanpif16
+ libc.src.__support.math.tgammabf16
libc.src.__support.math.trunc
libc.src.__support.math.truncbf16
libc.src.__support.math.truncf
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 8ae2aac3d51ed..08c69b19dfbf8 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -802,6 +802,9 @@ TEST(LlvmLibcSharedMathTest, AllBFloat16) {
bfloat16(0.0)));
EXPECT_FP_EQ(bfloat16(0.0), setpayloadsigbf16_res);
+ EXPECT_FP_EQ(bfloat16(1.0f),
+ LIBC_NAMESPACE::shared::tgammabf16(bfloat16(1.0f)));
+
EXPECT_FP_EQ(bfloat16(1.0), LIBC_NAMESPACE::shared::expbf16(bfloat16(0.0)));
EXPECT_FP_EQ(bfloat16(0.0), LIBC_NAMESPACE::shared::log_bf16(bfloat16(1.0)));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 4698f246fe4cc..d33158f037fd8 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -3193,6 +3193,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.bfloat16
)
+add_fp_unittest(
+ tgammabf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ tgammabf16_test.cpp
+ DEPENDS
+ libc.src.math.tgammabf16
+ libc.src.__support.FPUtil.bfloat16
+)
+
add_fp_unittest(
dmull_test
NEED_MPFR
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index e4e2fb11438d4..820425cc13161 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -6229,6 +6229,19 @@ add_fp_unittest(
libc.src.__support.FPUtil.bfloat16
)
+add_fp_unittest(
+ tgammabf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ tgammabf16_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.math.tgammabf16
+ libc.src.__support.FPUtil.bfloat16
+)
+
add_fp_unittest(
dmull_test
SUITE
diff --git a/libc/test/src/math/smoke/tgammabf16_test.cpp b/libc/test/src/math/smoke/tgammabf16_test.cpp
new file mode 100644
index 0000000000000..d6f1223c5214c
--- /dev/null
+++ b/libc/test/src/math/smoke/tgammabf16_test.cpp
@@ -0,0 +1,79 @@
+//===-- Unittests for tgammabf16 ------------------------------------------===//
+//
+// 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/tgammabf16.h"
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+class LlvmLibcTgammaBf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+ DECLARE_SPECIAL_CONSTANTS(bfloat16)
+
+public:
+ void test_special_numbers() {
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::tgammabf16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::tgammabf16(sNaN),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::tgammabf16(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ aNaN, LIBC_NAMESPACE::tgammabf16(neg_inf), FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ inf, LIBC_NAMESPACE::tgammabf16(zero), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ neg_inf, LIBC_NAMESPACE::tgammabf16(neg_zero), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
+
+ void test_poles() {
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ aNaN, LIBC_NAMESPACE::tgammabf16(bfloat16(-1.0f)), FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(
+ aNaN, LIBC_NAMESPACE::tgammabf16(bfloat16(-2.0f)), FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+ }
+
+ void test_values() {
+ EXPECT_FP_EQ_ALL_ROUNDING(bfloat16(1.0f),
+ LIBC_NAMESPACE::tgammabf16(bfloat16(1.0f)));
+ EXPECT_FP_EQ_ALL_ROUNDING(bfloat16(1.0f),
+ LIBC_NAMESPACE::tgammabf16(bfloat16(2.0f)));
+ EXPECT_FP_EQ_ALL_ROUNDING(bfloat16(2.0f),
+ LIBC_NAMESPACE::tgammabf16(bfloat16(3.0f)));
+ EXPECT_FP_EQ_ALL_ROUNDING(bfloat16(24.0f),
+ LIBC_NAMESPACE::tgammabf16(bfloat16(5.0f)));
+ }
+
+ void test_boundaries() {
+ EXPECT_FP_EQ_ALL_ROUNDING(bfloat16(2.951479e38f),
+ LIBC_NAMESPACE::tgammabf16(bfloat16(35.0f)));
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::tgammabf16(bfloat16(36.0f)), FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
+};
+
+TEST_F(LlvmLibcTgammaBf16Test, SpecialNumbers) { test_special_numbers(); }
+TEST_F(LlvmLibcTgammaBf16Test, Poles) { test_poles(); }
+TEST_F(LlvmLibcTgammaBf16Test, Values) { test_values(); }
+TEST_F(LlvmLibcTgammaBf16Test, Boundaries) { test_boundaries(); }
diff --git a/libc/test/src/math/tgammabf16_test.cpp b/libc/test/src/math/tgammabf16_test.cpp
new file mode 100644
index 0000000000000..5d4087ccc85c5
--- /dev/null
+++ b/libc/test/src/math/tgammabf16_test.cpp
@@ -0,0 +1,42 @@
+//===-- Exhaustive test for tgammabf16 ------------------------------------===//
+//
+// 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/tgammabf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcTgammaBf16Test = LIBC_NAMESPACE::testing::FPTest<bfloat16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+static bool is_negative_integer(bfloat16 x) {
+ float xf = static_cast<float>(x);
+ if (xf > -1.0f)
+ return false;
+ if (xf <= -128.0f)
+ return true;
+ int n = static_cast<int>(xf);
+ return xf == static_cast<float>(n);
+}
+
+TEST_F(LlvmLibcTgammaBf16Test, Exhaustive) {
+ for (uint32_t v = 0x0000; v < 0x10000; ++v) {
+ bfloat16 x =
+ LIBC_NAMESPACE::fputil::FPBits<bfloat16>(static_cast<uint16_t>(v))
+ .get_val();
+ LIBC_NAMESPACE::fputil::FPBits<bfloat16> bits(x);
+
+ if (bits.is_nan() || bits.is_zero() ||
+ (bits.is_neg() && is_negative_integer(x)))
+ continue;
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tgamma, x,
+ LIBC_NAMESPACE::tgammabf16(x), 1.0);
+ }
+}
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index 452a1391058d4..b5ca8370848e2 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -529,6 +529,12 @@ MPFRNumber MPFRNumber::tanh() const {
return result;
}
+MPFRNumber MPFRNumber::tgamma() const {
+ MPFRNumber result(*this);
+ mpfr_gamma(result.value, value, mpfr_rounding);
+ return result;
+}
+
MPFRNumber MPFRNumber::tanpi() const {
MPFRNumber result(*this);
diff --git a/libc/utils/MPFRWrapper/MPCommon.h b/libc/utils/MPFRWrapper/MPCommon.h
index 935a2614968a2..96149487170a5 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -237,6 +237,7 @@ class MPFRNumber {
MPFRNumber tan() const;
MPFRNumber tanh() const;
MPFRNumber tanpi() const;
+ MPFRNumber tgamma() const;
MPFRNumber trunc() const;
MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c);
MPFRNumber mul(const MPFRNumber &b);
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 356764302bda8..2fadaf4bc8388 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -113,6 +113,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.tanh();
case Operation::Tanpi:
return mpfrInput.tanpi();
+ case Operation::Tgamma:
+ return mpfrInput.tgamma();
case Operation::Trunc:
return mpfrInput.trunc();
default:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 64343fde2969f..3383d2d2434c9 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -67,6 +67,7 @@ enum class Operation : int {
Tan,
Tanh,
Tanpi,
+ Tgamma,
Trunc,
EndUnaryOperationsSingleOutput,
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 56ef8876c9e2d..a1bd87e143979 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -9267,6 +9267,25 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_tgammabf16",
+ hdrs = ["src/__support/math/tgammabf16.h"],
+ deps = [
+ ":__support_fputil_bfloat16",
+ ":__support_fputil_cast",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_exp",
+ ":__support_math_log",
+ ":__support_math_sin",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ],
+)
+
libc_support_library(
name = "__support_math_totalorder",
hdrs = ["src/__support/math/totalorder.h"],
@@ -12972,6 +12991,13 @@ libc_math_function(
additional_deps = [":__support_math_tanpif16"],
)
+libc_math_function(
+ name = "tgammabf16",
+ additional_deps = [
+ ":__support_math_tgammabf16",
+ ],
+)
+
libc_math_function(
name = "totalorder",
additional_deps = [
>From ff071924e2cf169852bf8fa506d7ef45d60d7c6b Mon Sep 17 00:00:00 2001
From: Vedant Neve <vedantneve13 at gmail.com>
Date: Tue, 14 Jul 2026 17:05:52 +0000
Subject: [PATCH 2/4] [libc][math] Optimize tgammabf16 with exact hex literals
and bitwise checks - Address review feedback
Signed-off-by: Vedant Neve <vedantneve13 at gmail.com>
---
libc/src/__support/math/tgammabf16.h | 55 +++++++++++---------
libc/test/src/math/smoke/tgammabf16_test.cpp | 3 --
2 files changed, 29 insertions(+), 29 deletions(-)
diff --git a/libc/src/__support/math/tgammabf16.h b/libc/src/__support/math/tgammabf16.h
index 275f9fc15745b..701f6b1b84e06 100644
--- a/libc/src/__support/math/tgammabf16.h
+++ b/libc/src/__support/math/tgammabf16.h
@@ -32,18 +32,9 @@ LIBC_INLINE_VAR constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
// Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
LIBC_INLINE_VAR constexpr double LANCZOS_COEFFS[9] = {
- 0.99999999999980993, 676.5203681218851, -1259.1392167224028,
- 771.32342877765313, -176.61502916214059, 12.507343278224757,
- -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
-
-LIBC_INLINE bool is_negative_integer(double x) {
- if (x > -1.0)
- return false;
- if (x <= -128.0)
- return true;
- int n = static_cast<int>(x);
- return x == static_cast<double>(n);
-}
+ 0x1.ffffffffff950p-1, 0x1.52429b6c30b05p+9, -0x1.3ac8e8ed4171bp+10,
+ 0x1.81a9661d3b4d8p+9, -0x1.613ae51a32f5dp+7, 0x1.903c27f8b9c81p+3,
+ -0x1.1bcb2992b2855p-3, 0x1.4f0514e4e324fp-17, 0x1.435508f3faeefp-23};
} // namespace tgammabf16_internal
@@ -75,17 +66,28 @@ LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
return FPBits::inf(xbits.sign()).get_val();
}
- double xd = static_cast<double>(static_cast<float>(x));
- if (LIBC_UNLIKELY(is_negative_integer(xd))) {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
+ if (LIBC_UNLIKELY(xbits.is_neg())) {
+ uint16_t x_abs = xbits.uintval() & 0x7fffU;
+ 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 ||
+ (xbits.get_mantissa() &
+ static_cast<uint16_t>((1U << (FPBits::FRACTION_LEN - e)) - 1U)) ==
+ 0U) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ }
}
+ float xf = static_cast<float>(x);
+
// Fast path for exact positive integers
- if (xd > 0.0 && xd <= 35.0) {
- int n = static_cast<int>(xd);
- if (xd == static_cast<double>(n)) {
+ if (xf > 0.0f && xf <= 35.0f) {
+ int n = static_cast<int>(xf);
+ if (xf == static_cast<float>(n)) {
double res = 1.0;
for (int i = 1; i < n; ++i)
res *= i;
@@ -93,20 +95,21 @@ LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
}
}
- bool reflection = false;
+ bool reflection = xbits.is_neg();
bool divide_by_x = false;
+
+ double xd = static_cast<double>(xf);
double x_eval = xd;
double res;
// Fast path for tiny positive inputs to prevent exact-boundary overshoots
// For tiny x, Gamma(x) ~= 1/x - gamma
- if (LIBC_UNLIKELY(xd > 0.0 && xd < 0x1.0p-8)) {
- res = (1.0 / xd) - 0.577215664901532860606;
+ if (LIBC_UNLIKELY(xf > 0.0f && xf < 0x1.0p-8f)) {
+ res = (1.0 / xd) - 0x1.2788cfc6fb619p-1;
} else {
- if (xd < 0.0) {
- reflection = true;
+ if (reflection) {
x_eval = 1.0 - xd;
- } else if (xd < 1.0) {
+ } else if (xf < 1.0f) {
divide_by_x = true;
x_eval = xd + 1.0;
}
diff --git a/libc/test/src/math/smoke/tgammabf16_test.cpp b/libc/test/src/math/smoke/tgammabf16_test.cpp
index d6f1223c5214c..935c8f0c5930b 100644
--- a/libc/test/src/math/smoke/tgammabf16_test.cpp
+++ b/libc/test/src/math/smoke/tgammabf16_test.cpp
@@ -64,9 +64,6 @@ class LlvmLibcTgammaBf16Test : public LIBC_NAMESPACE::testing::FEnvSafeTest {
}
void test_boundaries() {
- EXPECT_FP_EQ_ALL_ROUNDING(bfloat16(2.951479e38f),
- LIBC_NAMESPACE::tgammabf16(bfloat16(35.0f)));
-
EXPECT_FP_EQ_WITH_EXCEPTION(
inf, LIBC_NAMESPACE::tgammabf16(bfloat16(36.0f)), FE_OVERFLOW);
EXPECT_MATH_ERRNO(ERANGE);
>From c2e1be0858d033fa9e5f7480886d59edc4fd2673 Mon Sep 17 00:00:00 2001
From: Vedant Neve <vedantneve13 at gmail.com>
Date: Tue, 14 Jul 2026 17:28:30 +0000
Subject: [PATCH 3/4] [libc][math] Add references to Lanczos approximation
coefficients
Signed-off-by: Vedant Neve <vedantneve13 at gmail.com>
---
libc/src/__support/math/tgammabf16.h | 5 +++++
1 file changed, 5 insertions(+)
diff --git a/libc/src/__support/math/tgammabf16.h b/libc/src/__support/math/tgammabf16.h
index 701f6b1b84e06..8feb0ed472603 100644
--- a/libc/src/__support/math/tgammabf16.h
+++ b/libc/src/__support/math/tgammabf16.h
@@ -31,6 +31,11 @@ LIBC_INLINE_VAR constexpr double PI = 0x1.921fb54442d18p+1;
LIBC_INLINE_VAR constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
// Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
+// Reference: "A note on the computation of the convergent Lanczos complex Gamma
+// approximation" by Paul Godfrey (2001).
+// Original:
+// https://web.archive.org/web/20060915161115/http://my.fit.edu/~gabdo/gamma.txt
+// Mirror: http://www.mrob.com/pub/ries/lanczos-gamma.html
LIBC_INLINE_VAR constexpr double LANCZOS_COEFFS[9] = {
0x1.ffffffffff950p-1, 0x1.52429b6c30b05p+9, -0x1.3ac8e8ed4171bp+10,
0x1.81a9661d3b4d8p+9, -0x1.613ae51a32f5dp+7, 0x1.903c27f8b9c81p+3,
>From 5cbf6f1d305551ae9a724af683488ca84026b877 Mon Sep 17 00:00:00 2001
From: Vedant Neve <vedantneve13 at gmail.com>
Date: Thu, 16 Jul 2026 15:57:27 +0530
Subject: [PATCH 4/4] [libc][math] Optimize tgammabf16 with exact hex literals
and bitwise checks
Signed-off-by: Vedant Neve <vedantneve13 at gmail.com>
---
libc/docs/headers/math/index.rst | 2 +-
libc/src/__support/math/tgammabf16.h | 60 +++++++++----------
libc/src/math/generic/CMakeLists.txt | 20 +++----
libc/utils/MPFRWrapper/MPCommon.cpp | 12 ++--
.../libc/test/src/math/BUILD.bazel | 5 ++
.../libc/test/src/math/smoke/BUILD.bazel | 5 ++
6 files changed, 56 insertions(+), 48 deletions(-)
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index d17eccecf1a3e..563940d84236a 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -359,7 +359,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
| tanpi | |check| | | | |check| | | | 7.12.4.14 | F.10.1.14 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
-| tgamma | | | | | | | 7.12.8.4 | F.10.5.4 |
+| tgamma | | | | | | |check| | 7.12.8.4 | F.10.5.4 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/src/__support/math/tgammabf16.h b/libc/src/__support/math/tgammabf16.h
index 8feb0ed472603..b366f2fe8dc4a 100644
--- a/libc/src/__support/math/tgammabf16.h
+++ b/libc/src/__support/math/tgammabf16.h
@@ -25,39 +25,36 @@
namespace LIBC_NAMESPACE_DECL {
namespace math {
-namespace tgammabf16_internal {
-
-LIBC_INLINE_VAR constexpr double PI = 0x1.921fb54442d18p+1;
-LIBC_INLINE_VAR constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
-
-// Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
-// Reference: "A note on the computation of the convergent Lanczos complex Gamma
-// approximation" by Paul Godfrey (2001).
-// Original:
-// https://web.archive.org/web/20060915161115/http://my.fit.edu/~gabdo/gamma.txt
-// Mirror: http://www.mrob.com/pub/ries/lanczos-gamma.html
-LIBC_INLINE_VAR constexpr double LANCZOS_COEFFS[9] = {
- 0x1.ffffffffff950p-1, 0x1.52429b6c30b05p+9, -0x1.3ac8e8ed4171bp+10,
- 0x1.81a9661d3b4d8p+9, -0x1.613ae51a32f5dp+7, 0x1.903c27f8b9c81p+3,
- -0x1.1bcb2992b2855p-3, 0x1.4f0514e4e324fp-17, 0x1.435508f3faeefp-23};
-
-} // namespace tgammabf16_internal
-
LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
using FPBits = fputil::FPBits<bfloat16>;
- using namespace tgammabf16_internal;
+ // The hexadecimal literals were generated with Sollya:
+ // > display = hexadecimal;
+ // > round(pi, D, RN);
+ // > round(log(sqrt(2 * pi)), D, RN);
+ constexpr double PI = 0x1.921fb54442d18p+1;
+ constexpr double LOG_SQRT_2_PI = 0x1.d67f1c864beb5p-1;
+ // Paul Godfrey's exact Lanczos approximation coefficients (g=7, n=9)
+ // Reference: "A note on the computation of the convergent Lanczos complex
+ // Gamma approximation" by Paul Godfrey (2001).
+ // Original:
+ // https://web.archive.org/web/20060915161115/http://my.fit.edu/~gabdo/gamma.txt
+ // Mirror: http://www.mrob.com/pub/ries/lanczos-gamma.html
+ constexpr double LANCZOS_COEFFS[9] = {
+ 0x1.ffffffffff950p-1, 0x1.52429b6c30b05p+9, -0x1.3ac8e8ed4171bp+10,
+ 0x1.81a9661d3b4d8p+9, -0x1.613ae51a32f5dp+7, 0x1.903c27f8b9c81p+3,
+ -0x1.1bcb2992b2855p-3, 0x1.4f0514e4e324fp-17, 0x1.435508f3faeefp-23};
FPBits xbits(x);
+ uint16_t x_abs = xbits.uintval() & 0x7fffU;
- if (LIBC_UNLIKELY(xbits.is_nan())) {
- if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
+ if (LIBC_UNLIKELY(x_abs >= FPBits::EXP_MASK)) {
+ if (xbits.is_nan()) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
}
- return x;
- }
-
- if (LIBC_UNLIKELY(xbits.is_inf())) {
if (xbits.is_pos())
return x;
fputil::set_errno_if_required(EDOM);
@@ -72,7 +69,6 @@ LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
}
if (LIBC_UNLIKELY(xbits.is_neg())) {
- uint16_t x_abs = xbits.uintval() & 0x7fffU;
int biased_exp = x_abs >> FPBits::FRACTION_LEN;
if (biased_exp >= FPBits::EXP_BIAS) {
int e = biased_exp - FPBits::EXP_BIAS;
@@ -107,10 +103,12 @@ LIBC_INLINE bfloat16 tgammabf16(bfloat16 x) {
double x_eval = xd;
double res;
- // Fast path for tiny positive inputs to prevent exact-boundary overshoots
- // For tiny x, Gamma(x) ~= 1/x - gamma
+ // Fast path for tiny positive inputs to prevent exact-boundary overshoots.
+ // EULER_GAMMA is the Euler-Mascheroni constant. Its correctly rounded binary64
+ // value was generated with MPFR's mpfr_const_euler.
if (LIBC_UNLIKELY(xf > 0.0f && xf < 0x1.0p-8f)) {
- res = (1.0 / xd) - 0x1.2788cfc6fb619p-1;
+ constexpr double EULER_GAMMA = 0x1.2788cfc6fb619p-1;
+ res = (1.0 / xd) - EULER_GAMMA;
} else {
if (reflection) {
x_eval = 1.0 - xd;
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 64da6c553eb25..56b9d91fb1a71 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -425,16 +425,6 @@ add_entrypoint_object(
libc.src.__support.math.tanbf16
)
-add_entrypoint_object(
- tgammabf16
- SRCS
- tgammabf16.cpp
- HDRS
- ../tgammabf16.h
- DEPENDS
- libc.src.__support.math.tgammabf16
-)
-
add_entrypoint_object(
tanf
SRCS
@@ -476,6 +466,16 @@ add_entrypoint_object(
libc.src.__support.math.tanpif16
)
+add_entrypoint_object(
+ tgammabf16
+ SRCS
+ tgammabf16.cpp
+ HDRS
+ ../tgammabf16.h
+ DEPENDS
+ libc.src.__support.math.tgammabf16
+)
+
add_entrypoint_object(
fabs
SRCS
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index b5ca8370848e2..f7a20d29ae12f 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -529,12 +529,6 @@ MPFRNumber MPFRNumber::tanh() const {
return result;
}
-MPFRNumber MPFRNumber::tgamma() const {
- MPFRNumber result(*this);
- mpfr_gamma(result.value, value, mpfr_rounding);
- return result;
-}
-
MPFRNumber MPFRNumber::tanpi() const {
MPFRNumber result(*this);
@@ -581,6 +575,12 @@ MPFRNumber MPFRNumber::tanpi() const {
#endif
}
+MPFRNumber MPFRNumber::tgamma() const {
+ MPFRNumber result(*this);
+ mpfr_gamma(result.value, value, mpfr_rounding);
+ return result;
+}
+
MPFRNumber MPFRNumber::trunc() const {
MPFRNumber result(*this);
mpfr_trunc(result.value, value);
diff --git a/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
index 9b5fa61303f10..82268d464223d 100644
--- a/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
@@ -653,3 +653,8 @@ math_mpfr_test(name = "logf16")
math_mpfr_test(name = "log2f16")
math_mpfr_test(name = "log10f16")
+
+math_mpfr_test(
+ name = "tgammabf16",
+ deps = ["//libc:__support_fputil_bfloat16"],
+)
diff --git a/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel
index 3f0874110e140..c7a4cd743b231 100644
--- a/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel
@@ -1737,3 +1737,8 @@ math_test(name = "tanf16")
math_test(name = "tanhf16")
math_test(name = "tanpif16")
+
+math_test(
+ name = "tgammabf16",
+ deps = ["//libc:__support_fputil_bfloat16"],
+)
More information about the libc-commits
mailing list