[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:47:25 PDT 2026
https://github.com/udaykiriti updated https://github.com/llvm/llvm-project/pull/199312
>From fe433888b4c2a41e43c9dc34b2a091b32e35dbf2 Mon Sep 17 00:00:00 2001
From: udaykiriti <udaykiriti9 at gmail.com>
Date: Fri, 24 Jul 2026 15:17:10 +0530
Subject: [PATCH] [libc][math] Implement lgammabf16 function
Adds bfloat16 lgamma implementation, fixes double-rounding for
small positive inputs, cleans up cmake/bazel deps, and addresses
review comments.
---
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 | 20 +-
libc/src/__support/math/lgammabf16.h | 260 ++++++++++++++++++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/generic/CMakeLists.txt | 10 +
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 | 79 ++++++
libc/test/src/math/smoke/CMakeLists.txt | 13 +
libc/test/src/math/smoke/lgammabf16_test.cpp | 106 +++++++
.../llvm-project-overlay/libc/BUILD.bazel | 25 ++
27 files changed, 602 insertions(+), 3 deletions(-)
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..e6dbc89cfd2ae 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -61,6 +61,24 @@ add_header_library(
libc.src.__support.macros.properties.types
)
+add_header_library(
+ lgammabf16
+ HDRS
+ lgammabf16.h
+ DEPENDS
+ 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.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
+)
+
add_header_library(
acosh_float_constants
HDRS
@@ -4671,8 +4689,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
diff --git a/libc/src/__support/math/lgammabf16.h b/libc/src/__support/math/lgammabf16.h
new file mode 100644
index 0000000000000..a4fb097e1c945
--- /dev/null
+++ b/libc/src/__support/math/lgammabf16.h
@@ -0,0 +1,260 @@
+//===-- 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___SUPPORT_MATH_LGAMMABF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMABF16_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/log.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace lgammabf16_internal {
+
+// lgamma_positive_d: compute lgamma(x) for x > 0, returning double.
+//
+// Takes double so callers can pass (1.0 + ax) without float precision loss.
+// The double return is necessary to avoid double-rounding: the final
+// fputil::cast<bfloat16> needs the full double precision so it can correctly
+// break ties (see the 0x35E5 example in lgammabf16 below).
+//
+// For x < 4, applies the recurrence lgamma(x) = lgamma(x+1) - ln(x) until
+// x reaches [4, 8), then evaluates the polynomial. This is critical because
+// the [2,3) polynomial has max_rel_err=9.98e-6 which, near its edges (t near
+// +0.5 or -0.5), causes ~6e-6 absolute error. After subtracting ln(x), the
+// result near x=1 or x=2 can be as small as ~0.002, giving ~0.45 ULP error --
+// which fails the directed-rounding tolerance test. Polynomials for [4,5) and
+// above have max_rel_err <= 6.61e-7, keeping final error well under 0.1 ULP.
+//
+LIBC_INLINE double lgamma_positive_d(double x) {
+ // Coefficients for lgamma on [n, n+1), centered at n+0.5.
+ // Each row: {c0, c1, c2, c3, c4} for Estrin evaluation
+ // where t = x - (n + 0.5), n = 1..7.
+ // 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.
+ constexpr float LGAMMA_POLY[7][5] = {
+ // [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))
+ return 0.0;
+
+ if (x >= 8.0) {
+ // Stirling series; 0.5*ln(2*pi)
+ // > 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;
+ result += 1.0 / (12.0 * x) - 1.0 / (360.0 * x * x2);
+ return result;
+ }
+
+ // For x in (0, 4): apply recurrence relation
+ // lgamma(x) = lgamma(x+n) - ln(x*(x+1)*...*(x+n-1))
+ // to shift x into the stable [4, 8) range for polynomial evaluation.
+ // x already in [4, 8) needs no shift -- evaluate the polynomial directly.
+ double log_product, xs, product;
+
+ 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) {
+ product = x * (x + 1.0);
+ log_product = math::log(product);
+ xs = x + 2.0;
+ } else if (x >= 1.0) {
+ product = x * (x + 1.0);
+ product = product * (x + 2.0);
+ log_product = math::log(product);
+ xs = x + 3.0;
+ } else {
+ product = x * (x + 1.0);
+ product = product * (x + 2.0);
+ product = product * (x + 3.0);
+ log_product = math::log(product);
+ xs = x + 4.0;
+ }
+
+ // xs in [4, 8); select polynomial interval.
+ float xf = static_cast<float>(xs);
+ int n = static_cast<int>(xf);
+ if (n >= 7)
+ n = 7;
+ float t = xf - (static_cast<float>(n) + 0.5f);
+
+ float c0 = LGAMMA_POLY[n - 1][0];
+ float c1 = LGAMMA_POLY[n - 1][1];
+ float c2 = LGAMMA_POLY[n - 1][2];
+ float c3 = LGAMMA_POLY[n - 1][3];
+ float c4 = LGAMMA_POLY[n - 1][4];
+
+ // Estrin's scheme for p(t) = c0 + c1*t + c2*t^2 + c3*t^3 + c4*t^4:
+ // p(t) = (c0 + c1*t) + t^2 * ((c2 + c3*t) + t^2 * c4)
+ float t2 = t * t;
+ float p01 = fputil::multiply_add(t, c1, c0);
+ float p23 = fputil::multiply_add(t, c3, c2);
+ float p234 = fputil::multiply_add(t2, c4, p23);
+ float lgamma_xs_f = fputil::multiply_add(t2, p234, p01);
+ double lgamma_xs = static_cast<double>(lgamma_xs_f);
+
+ return lgamma_xs - log_product;
+}
+} // namespace lgammabf16_internal
+
+LIBC_INLINE bfloat16 lgammabf16(bfloat16 x) {
+ using FPBits = fputil::FPBits<bfloat16>;
+ FPBits x_bits(x);
+
+ // Handles 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)
+ // > mpmath.mp.dps = 50; hex(mpmath.log(mpmath.pi))
+ constexpr double LN_PI_D = 0x1.250d048e7a1bdp+0;
+ float ax = -xf;
+
+ // 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) 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;
+
+ // 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!
+ -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 =
+ 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->bfloat16 (not via float) to avoid
+ // double-rounding at bfloat16 tie points.
+ return fputil::cast<bfloat16>(result_d);
+ }
+
+ // Cast directly from double->bfloat16 for the same reason as above.
+ return fputil::cast<bfloat16>(
+ lgammabf16_internal::lgamma_positive_d(static_cast<double>(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..31930682eded1 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4942,6 +4942,16 @@ add_entrypoint_object(
libc.src.__support.math.bf16fmal
)
+add_entrypoint_object(
+ lgammabf16
+ SRCS
+ lgammabf16.cpp
+ HDRS
+ ../lgammabf16.h
+ DEPENDS
+ libc.src.__support.math.lgammabf16
+)
+
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..c35758e6ffd1b
--- /dev/null
+++ b/libc/test/src/math/lgammabf16_test.cpp
@@ -0,0 +1,79 @@
+//===-- 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 "src/__support/FPUtil/bfloat16.h"
+#include "src/math/lgammabf16.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;
+
+// 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;
+
+// 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(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);
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 81f76fa7681fd..996c2d7b419a3 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -36,6 +36,19 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ lgammabf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ lgammabf16_test.cpp
+ DEPENDS
+ libc.src.math.lgammabf16
+ libc.src.__support.FPUtil.bfloat16
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+)
+
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..193f66eafa197
--- /dev/null
+++ b/libc/test/src/math/smoke/lgammabf16_test.cpp
@@ -0,0 +1,106 @@
+//===-- 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() {
+ // 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);
+ }
+ 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/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index bee6f7aa5766c..5ece4758c41a2 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7094,6 +7094,24 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_lgammabf16",
+ hdrs = ["src/__support/math/lgammabf16.h"],
+ deps = [
+ ":__support_fputil_bfloat16",
+ ":__support_fputil_cast",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_log",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ],
+)
+
libc_support_library(
name = "__support_math_llogbbf16",
hdrs = ["src/__support/math/llogbbf16.h"],
@@ -12369,6 +12387,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 = [
More information about the libc-commits
mailing list