[libc-commits] [libc] [llvm] [libc][math] Implement single precision lgamma function (PR #205490)
via libc-commits
libc-commits at lists.llvm.org
Tue Jun 23 23:22:28 PDT 2026
https://github.com/AnonMiraj created https://github.com/llvm/llvm-project/pull/205490
The implementation reduces to `|x|` and selects a path
- **Tiny** ($|x| < 2^{-23}$): truncated Laurent series $\text{lgamma}(x) = -\log|x| - \gamma\,|x|$.
- **Small** ($|x| < 0.66$): degree-22 Chebyshev for $h(t) = \frac{\text{lgamma}(t) + \log(t)}{t}$, recovering $\text{lgamma}(t) = t \cdot h(t) - \log(t)$.
- **Medium** ($|x| \in [0.66, 3.37]$): degree-22 Chebyshev split into three sub-ranges that factor out the function's roots —
- $[0.66, 1)$: $\text{lgamma}(t) = (t-1)\,P(t)$
- $[1, 2)$: $\text{lgamma}(t) = (t-1)(t-2)\,P(t)$
- $[2, 3.37)$: $\text{lgamma}(t) = (t-2)\,P(t)$
- **Stirling** ($|x| > 3.37$): $(x-0.5)\log(x) - x + \frac{\log(2\pi)}{2}$ plus 1/2/4/8-term Bernoulli corrections selected by sub-range (the dominant term is evaluated in double-double for large $|x|$).
Negative $x$ uses the reflection formula $\text{lgamma}(x) = \log\pi - \log|\sin(\pi x)| - \text{lgamma}(|x|)$.
Exhaustive tests pass in all 4 rounding modes.
```
-- CORE-MATH reciprocal throughput --
Ntrial = 20 ; Min = 54.476 + 7.003 clc/call; Median-Min = 1.211 clc/call; Max = 67.906 clc/call;
-- System LIBC reciprocal throughput --
Ntrial = 20 ; Min = 70.125 + 8.204 clc/call; Median-Min = 4.416 clc/call; Max = 90.469 clc/call;
```
Depends on #192834 so lgammaf16 must be merged first.
>From 9a59d5aeb8bafe962e0d2d87f4f4088c0c5ca303 Mon Sep 17 00:00:00 2001
From: Anonmiraj <ezzibrahimx at gmail.com>
Date: Sun, 19 Apr 2026 08:32:58 +0200
Subject: [PATCH 1/3] [libc][math] Implement half precision lgamma function
---
libc/config/baremetal/aarch64/entrypoints.txt | 1 +
libc/config/baremetal/arm/entrypoints.txt | 1 +
libc/config/baremetal/riscv/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/include/math.yaml | 7 +
libc/shared/math.h | 1 +
libc/shared/math/lgammaf16.h | 27 ++
libc/src/__support/math/CMakeLists.txt | 18 +
libc/src/__support/math/lgammaf16.h | 311 ++++++++++++++++++
libc/src/math/CMakeLists.txt | 2 +
libc/src/math/generic/CMakeLists.txt | 13 +
libc/src/math/generic/lgammaf16.cpp | 20 ++
libc/src/math/lgammaf16.h | 21 ++
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 1 +
libc/test/src/math/CMakeLists.txt | 11 +
libc/test/src/math/lgammaf16_test.cpp | 42 +++
libc/test/src/math/smoke/CMakeLists.txt | 11 +
libc/test/src/math/smoke/lgammaf16_test.cpp | 61 ++++
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 ++
26 files changed, 587 insertions(+)
create mode 100644 libc/shared/math/lgammaf16.h
create mode 100644 libc/src/__support/math/lgammaf16.h
create mode 100644 libc/src/math/generic/lgammaf16.cpp
create mode 100644 libc/src/math/lgammaf16.h
create mode 100644 libc/test/src/math/lgammaf16_test.cpp
create mode 100644 libc/test/src/math/smoke/lgammaf16_test.cpp
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index dcb50135232e2..c2d9e860388f9 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -655,6 +655,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index fac62bac939cc..6c80f317917e2 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -667,6 +667,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index a3b96225ff09d..3dfa9013e6623 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -663,6 +663,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index c7458ef36c941..9bb3cfa06923f 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -595,6 +595,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf16
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index e3336420cbbcc..89b931faf29c2 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -604,6 +604,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf16
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 3c3b0e835429f..b9e73d95ec244 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -664,6 +664,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf16
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index 3984a665b234b..67d19399b99ae 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -1682,6 +1682,13 @@ functions:
arguments:
- type: long double
- type: int *
+ - name: lgammaf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: llogb
standards:
- stdc
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 2baeb07294a3b..b9060dc5f6245 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -305,6 +305,7 @@
#include "math/ldexpf128.h"
#include "math/ldexpf16.h"
#include "math/ldexpl.h"
+#include "math/lgammaf16.h"
#include "math/llogb.h"
#include "math/llogbbf16.h"
#include "math/llogbf.h"
diff --git a/libc/shared/math/lgammaf16.h b/libc/shared/math/lgammaf16.h
new file mode 100644
index 0000000000000..4ae8d36e5ba53
--- /dev/null
+++ b/libc/shared/math/lgammaf16.h
@@ -0,0 +1,27 @@
+//===-- Shared lgammaf16 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_LGAMMAF16_H
+#define LLVM_LIBC_SHARED_MATH_LGAMMAF16_H
+
+#ifdef LIBC_TYPES_HAS_FLOAT16
+
+#include "shared/libc_common.h"
+#include "src/__support/math/lgammaf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::lgammaf16;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_TYPES_HAS_FLOAT16
+
+#endif // LLVM_LIBC_SHARED_MATH_LGAMMAF16_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 9e3ec26cdc881..a7ff2cec529b5 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3881,6 +3881,24 @@ add_header_library(
libc.include.llvm-libc-types.float128
)
+add_header_library(
+ lgammaf16
+ HDRS
+ lgammaf16.h
+ DEPENDS
+ 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
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer_operations
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_header_library(
llogbf16
HDRS
diff --git a/libc/src/__support/math/lgammaf16.h b/libc/src/__support/math/lgammaf16.h
new file mode 100644
index 0000000000000..2e6424fc4826c
--- /dev/null
+++ b/libc/src/__support/math/lgammaf16.h
@@ -0,0 +1,311 @@
+//===-- Implementation header for lgammaf16 ---------------------*- 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_LGAMMAF16_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H
+
+#include "include/llvm-libc-macros/float16-macros.h"
+
+#ifdef LIBC_TYPES_HAS_FLOAT16
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/NearestIntegerOperations.h"
+#include "src/__support/FPUtil/PolyEval.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/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace lgammaf16_internal {
+
+LIBC_INLINE constexpr bool is_integer(float16 x) {
+ using FPBits = fputil::FPBits<float16>;
+ FPBits xbits(x);
+ uint16_t x_u = xbits.uintval();
+ unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
+ unsigned lsb = static_cast<unsigned>(
+ cpp::countr_zero(static_cast<uint16_t>(x_u | FPBits::EXP_MASK)));
+ constexpr unsigned UNIT_EXPONENT =
+ static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
+ return x_e + lsb >= UNIT_EXPONENT;
+}
+
+// sin(pi*x) for x in (0, 1), written as (0.25 - u^2) * P(u^2) where u = x-0.5.
+// Coefficients for P are a degree-7 polynomial in u^2
+// approximating sin(pi*(u+0.5)) / (pi*(0.25-u^2))
+// with max error ~2^{-55}, generated by Sollya with:
+// > P = fpminimax(sin(pi*(x+0.5))/(pi*(0.25-x^2)),
+// [|0,2,4,6,8,10,12,14|], [|D...|], [-0.5, 0.5]);
+LIBC_INLINE constexpr double lg_sinpi(double x) {
+ constexpr double COEFFS[8] = {0x1p+2,
+ -0x1.de9e64df22ea4p+1,
+ 0x1.472be122401f8p+0,
+ -0x1.d4fcd82df91bp-3,
+ 0x1.9f05c97e0aab2p-6,
+ -0x1.f3091c427b611p-10,
+ 0x1.b22c9bfdca547p-14,
+ -0x1.15484325ef569p-18};
+ double u = x - 0.5;
+ double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
+ double p01 = fputil::multiply_add(u2, COEFFS[1], COEFFS[0]);
+ double p23 = fputil::multiply_add(u2, COEFFS[3], COEFFS[2]);
+ double p45 = fputil::multiply_add(u2, COEFFS[5], COEFFS[4]);
+ double p67 = fputil::multiply_add(u2, COEFFS[7], COEFFS[6]);
+ double p03 = fputil::multiply_add(u4, p23, p01);
+ double p47 = fputil::multiply_add(u4, p67, p45);
+ return (0.25 - u2) * fputil::multiply_add(u8, p47, p03);
+}
+
+// Natural logarithm of x (x > 0), using 16-entry table + degree-7 polynomial.
+// Range reduction: x = 2^e * m where m in [1, 2), decomposed as
+// m = (1+i/16)*(1+z) so log(x) = e*log(2) + log(1+i/16) + log(1+z)
+// = e*log(2) + IL[i] + z*P(z).
+// P approximates log(1+z)/z on z in [-1/16, 1/16] with max
+// error ~2^{-54}, generated by Sollya with:
+// > P = fpminimax(log(1+x)/x, [|0,1,2,3,4,5,6,7|], [|D...|], [-1/16, 1/16]);
+LIBC_INLINE double lg_ln(double x) {
+ using FPBits = fputil::FPBits<double>;
+ uint64_t u = FPBits(x).uintval();
+ int e = static_cast<int>(FPBits(x).get_biased_exponent()) - 0x3ff;
+
+ // Coefficients for log(1 + z)/z on z in [-1/16, 1/16]
+ constexpr double COEFFS[8] = {0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2,
+ 0x1.55555537802dep-2, -0x1.ffffeca81b866p-3,
+ 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3,
+ 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4};
+ // IL[i] = log(1 + i/16) for i = 0..15
+ constexpr double IL[16] = {
+ 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4,
+ 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2,
+ 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2,
+ 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1,
+ 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1,
+ 0x1.52a2d265bc5abp-1};
+ // IX[i] = 1 / (1 + i/16) for i = 0..15
+ constexpr double IX[16] = {
+ 0x000000000000001p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1,
+ 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1,
+ 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1,
+ 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1,
+ 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1,
+ 0x1.0842108421084p-1};
+
+ int i = static_cast<int>((u >> 48) & 0xf);
+ // Reduce to mantissa in [1, 2)
+ uint64_t mant_u = (u & (~uint64_t(0) >> 12)) | (uint64_t(0x3ff) << 52);
+ double mant = FPBits(mant_u).get_val();
+ double z = IX[i] * mant - 1.0, z2 = z * z, z4 = z2 * z2;
+ double q01 = fputil::multiply_add(z, COEFFS[1], COEFFS[0]);
+ double q23 = fputil::multiply_add(z, COEFFS[3], COEFFS[2]);
+ double q45 = fputil::multiply_add(z, COEFFS[5], COEFFS[4]);
+ double q67 = fputil::multiply_add(z, COEFFS[7], COEFFS[6]);
+ double q03 = fputil::multiply_add(z2, q23, q01);
+ double q47 = fputil::multiply_add(z2, q67, q45);
+ return e * 0x1.62e42fefa39efp-1 + IL[i] +
+ z * fputil::multiply_add(z4, q47, q03);
+}
+
+} // namespace lgammaf16_internal
+
+LIBC_INLINE float16 lgammaf16(float16 x) {
+ using namespace lgammaf16_internal;
+ using FPBits = fputil::FPBits<float16>;
+
+ FPBits xbits(x);
+ uint16_t x_abs = xbits.abs().uintval();
+
+ // Handle NaN and Inf
+ if (LIBC_UNLIKELY(x_abs >= 0x7c00u)) {
+ if (x_abs == 0x7c00u)
+ return FPBits::inf().get_val();
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ // +-0 pole error
+ if (LIBC_UNLIKELY(x_abs == 0)) {
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ fputil::set_errno_if_required(ERANGE);
+ return FPBits::inf().get_val();
+ }
+
+ if (LIBC_UNLIKELY(is_integer(x))) {
+ if (xbits.is_neg()) {
+ // pole -> +Inf
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ fputil::set_errno_if_required(ERANGE);
+ return FPBits::inf().get_val();
+ }
+ // lgamma(1) = lgamma(2) = 0
+ if (x_abs == 0x3c00u || x_abs == 0x4000u)
+ return FPBits::zero().get_val();
+ }
+
+ double xd = fputil::cast<double>(x);
+ double abs_xd = xd < 0.0 ? -xd : xd;
+ double lgamma_val;
+
+ if (LIBC_UNLIKELY(x_abs < 0x3948u)) {
+ // Small path: |x| < 0.66015625
+ // For t = |x|: lgamma(t) = t*P(u) - log(t), where P(u) approximates
+ // h(t) = (lgamma(t) + log(t)) / t, u = (t - MID_S) / HW_S in [-1, 1]
+ // For x < 0: lgamma(-t) = log(pi) - log(sin(pi*t)) - lgamma(t) - log(t)
+ // = log(pi) - lg_ln(lg_sinpi(t)) - t*P(u)
+ // Degree-12 Chebyshev approximation for h in u; max lgamma error ~1.5e-13
+ // Coefficients generated by DCT at 13 Chebyshev nodes on
+ // [1e-5, 0.66015625].
+ // Evaluated via Clenshaw: sum_{k=0}^{12} CHEB_H[k] * T_k(u)
+ constexpr double MID_S = 0x1.52014f8b588e3p-2;
+ constexpr double HW_S = 0x1.51feb074a771dp-2;
+ constexpr double CHEB_H[13] = {
+ -0x1.6ab0351d4e02fp-2, 0x1.ac55f136e49cep-3, -0x1.9fc5280cd41bfp-7,
+ 0x1.1825d37c81c1ep-10, -0x1.ae2f107c31c9ep-14, 0x1.6174bdf95e9abp-17,
+ -0x1.2e6abd22bc4ecp-20, 0x1.09b66d4f14b06p-23, -0x1.dbb938f189d8ap-27,
+ 0x1.afcb748e20e4ep-30, -0x1.8c221a9d89d8ap-33, 0x1.6e708e49b3aeep-36,
+ -0x1.50accec4ec4ecp-39};
+ // Clenshaw's algorithm: p = sum_k CHEB_H[k] * T_k(u_s)
+ double u_s = (abs_xd - MID_S) / HW_S, u_s2 = 2.0 * u_s;
+ double bh1 = CHEB_H[12], bh2 = 0.0;
+ for (int k = 11; k >= 1; k--) {
+ double bh0 = fputil::multiply_add(u_s2, bh1, CHEB_H[k] - bh2);
+ bh2 = bh1;
+ bh1 = bh0;
+ }
+ double poly_h = fputil::multiply_add(u_s, bh1, CHEB_H[0] - bh2);
+ // abs_xd * poly_h = lgamma(abs_xd) + log(abs_xd)
+ double poly_val = abs_xd * poly_h;
+ if (xbits.is_neg())
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lg_ln(lg_sinpi(abs_xd)) - poly_val;
+ else
+ lgamma_val = poly_val - lg_ln(abs_xd);
+ } else if (x_abs > 0x42BEu) {
+ double log_abs_xd = lg_ln(abs_xd);
+ // lgamma(x) = (x-0.5)*log(x) - x + log(2*pi)/2 + Bernoulli corrections
+ // Stirling base: (x-0.5)*(log(x)-1) + (log(2*pi)-1)/2
+ lgamma_val = (abs_xd - 0.5) * (log_abs_xd - 1.0) + 0x1.acfe390c97d69p-2;
+
+ double inv_x = 1.0 / abs_xd, inv_x2 = inv_x * inv_x;
+
+ if (x_abs > 0x64AEu) {
+ // |x| > 1198 1-term Stirling-Bernoulli: 1/(12*x)
+ lgamma_val += inv_x * 0x1.5555555555555p-4;
+ } else if (x_abs > 0x549Eu) {
+ // |x| in (73.875, 1198] 2-term correction
+ // Minimax polynomial in 1/x^2 for R(x) = lgamma(x) - stirling_base(x),
+ // where stirling_base(x) = (x-0.5)*log(x) - x + log(2*pi)/2
+ // Exact Bernoulli: B_2/(1*2) = 1/12, B_4/(3*4) = -1/360
+ // Generated by Sollya with:
+ // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2);
+ // }; > P = fpminimax(stir(x), [|1,3|], [|D...|], [73.875, 1198]);
+ constexpr double BERN2[2] = {0x1.555555547fbadp-4, -0x1.6c0fd270c465p-9};
+ lgamma_val += inv_x * fputil::polyeval(inv_x2, BERN2[0], BERN2[1]);
+ } else if (x_abs > 0x4955u) {
+ // |x| in (10.664, 73.875] 4-term correction
+ // Minimax polynomial in 1/x^2 for R(x); exact Bernoulli terms are
+ // 1/12, -1/360, 1/1260, -1/1680
+ // Generated by Sollya with:
+ // > display = hexadecimal;
+ // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2);
+ // }; > P = fpminimax(stir(x), [|1,3,5,7|], [|D...|], [10.664, 73.875]);
+ constexpr double BERN4[4] = {0x1.555555554de0bp-4, -0x1.6c16bdc45944fp-9,
+ 0x1.a0077f300ecb3p-11,
+ -0x1.2e9cfff3b29c2p-11};
+ double inv_x4 = inv_x2 * inv_x2;
+ lgamma_val +=
+ inv_x * fputil::multiply_add(
+ inv_x4, fputil::multiply_add(inv_x2, BERN4[3], BERN4[2]),
+ fputil::multiply_add(inv_x2, BERN4[1], BERN4[0]));
+ } else {
+ // |x| in (3.373, 10.664]: 8-term correction.
+ // Degree-7 minimax polynomial in 1/x^2 for R(x), generated by Sollya
+ // with:
+ // > stir = proc(x) { return lgamma(x)-((x-0.5)*log(x)-x+log(2*pi)/2);
+ // }; > P = fpminimax(stir(x), [|1,3,5,7,9,11,13,15|], [|D...|],
+ // [3.373, 10.664]);
+ constexpr double BERN8[8] = {
+ 0x1.5555555551286p-4, -0x1.6c16c0e7c4cf4p-9, 0x1.a0193267fe6f2p-11,
+ -0x1.37e87ec19cb45p-11, 0x1.b40011dfff081p-11, -0x1.c16c8946b19b6p-10,
+ 0x1.e9f47ace150d8p-9, -0x1.4f5843a71a338p-8};
+ double inv_x4 = inv_x2 * inv_x2, inv_x8 = inv_x4 * inv_x4;
+ double p01 = fputil::multiply_add(inv_x2, BERN8[1], BERN8[0]);
+ double p23 = fputil::multiply_add(inv_x2, BERN8[3], BERN8[2]);
+ double p45 = fputil::multiply_add(inv_x2, BERN8[5], BERN8[4]);
+ double p67 = fputil::multiply_add(inv_x2, BERN8[7], BERN8[6]);
+ double p03 = fputil::multiply_add(inv_x4, p23, p01);
+ double p47 = fputil::multiply_add(inv_x4, p67, p45);
+ lgamma_val += inv_x * fputil::multiply_add(inv_x8, p47, p03);
+ }
+
+ if (xbits.is_neg()) {
+ // Reflection formula
+ // lgamma(x) = log(pi) - log|x| - log|sin(pi*frac_x)| - lgamma|x|
+ double frac_x = xd - fputil::floor(xd);
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val - log_abs_xd;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x));
+ }
+ } else {
+ // Medium path: |x| in [0.66015625, 3.373046875]
+ // lgamma(|x|) = (|x|-1)*(|x|-2) * P(u), where P(u) approximates
+ // g(t) = lgamma(t) / ((t-1)*(t-2)), u = (t - MID_M) / HW_M in [-1, 1]
+ // Degree-20 Chebyshev approximation for g in u; max lgamma error ~1.4e-10
+ // Coefficients generated by DCT at 21 Chebyshev nodes on
+ // [0.66015625, 3.373046875]. Evaluated via Clenshaw: sum_{k=0}^{20}
+ // CHEB_G[k] * T_k(u)
+ constexpr double MID_M = 0x1.0220000000000p+1;
+ constexpr double HW_M = 0x1.5b40000000000p+0;
+ constexpr double CHEB_G[21] = {
+ 0x1.d690d4de8b355p-2, -0x1.53439726f7d8bp-3, 0x1.568fc6653c026p-5,
+ -0x1.8d074f73f36efp-7, 0x1.eefebf4ca8e86p-9, -0x1.4275913d85ff5p-10,
+ 0x1.b075f0d96e492p-12, -0x1.27ed35f357339p-13, 0x1.9b12cbe71a800p-15,
+ -0x1.20c949c2829b5p-16, 0x1.996eb95241555p-18, -0x1.2460d6c60f800p-19,
+ 0x1.a42297aae1862p-21, -0x1.2f6c7d28230a5p-22, 0x1.b8394aac30c31p-24,
+ -0x1.409726ce1110ep-25, 0x1.d484ac42db6dbp-27, -0x1.574dae55861abp-28,
+ 0x1.f782605555555p-30, -0x1.6cefb8f65aa2ep-31, 0x1.d9f268f3cf3cfp-33};
+ // Clenshaw's algorithm: poly_g = sum_k CHEB_G[k] * T_k(u_m)
+ double u_m = (abs_xd - MID_M) / HW_M, u_m2 = 2.0 * u_m;
+ double bg1 = CHEB_G[20], bg2 = 0.0;
+ for (int k = 19; k >= 1; k--) {
+ double bg0 = fputil::multiply_add(u_m2, bg1, CHEB_G[k] - bg2);
+ bg2 = bg1;
+ bg1 = bg0;
+ }
+ double poly_g = fputil::multiply_add(u_m, bg1, CHEB_G[0] - bg2);
+ lgamma_val = (abs_xd - 1.0) * (abs_xd - 2.0) * poly_g;
+
+ if (xbits.is_neg()) {
+ // Reflection
+ // (x) = log(pi) - lgamma|x| - log(|x|*|sin(pi*frac_x)|)
+ double frac_x = xd - fputil::floor(xd);
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd);
+ }
+ }
+
+ float16 result = fputil::cast<float16>(lgamma_val);
+ if (LIBC_UNLIKELY(fputil::FPBits<float16>(result).is_inf()))
+ fputil::set_errno_if_required(ERANGE);
+ return result;
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LIBC_TYPES_HAS_FLOAT16
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF16_H
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index b53817e2a1729..8d2a981679980 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -376,6 +376,8 @@ add_math_entrypoint_object(issignalingf16)
add_math_entrypoint_object(issignalingf128)
add_math_entrypoint_object(issignalingbf16)
+add_math_entrypoint_object(lgammaf16)
+
add_math_entrypoint_object(llogb)
add_math_entrypoint_object(llogbf)
add_math_entrypoint_object(llogbl)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 7ccbddba07b8d..49ca911ef3005 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1567,6 +1567,19 @@ add_entrypoint_object(
libc.src.__support.math.ilogbbf16
)
+add_entrypoint_object(
+ lgammaf16
+ SRCS
+ lgammaf16.cpp
+ HDRS
+ ../lgammaf16.h
+ DEPENDS
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.math.lgammaf16
+ libc.src.errno.errno
+)
+
add_entrypoint_object(
llogb
SRCS
diff --git a/libc/src/math/generic/lgammaf16.cpp b/libc/src/math/generic/lgammaf16.cpp
new file mode 100644
index 0000000000000..91b14aea3b01f
--- /dev/null
+++ b/libc/src/math/generic/lgammaf16.cpp
@@ -0,0 +1,20 @@
+//===-- Implementation of lgammaf16 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/lgammaf16.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/math/lgammaf16.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(float16, lgammaf16, (float16 x)) {
+ return math::lgammaf16(x);
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/lgammaf16.h b/libc/src/math/lgammaf16.h
new file mode 100644
index 0000000000000..6d233848ca1ea
--- /dev/null
+++ b/libc/src/math/lgammaf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for lgammaf16 ---------------------*- 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_LGAMMAF16_H
+#define LLVM_LIBC_SRC_MATH_LGAMMAF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 lgammaf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_LGAMMAF16_H
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index 84a8d4bf79b3d..68907f9561b6d 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -332,6 +332,7 @@ add_fp_unittest(
libc.src.__support.math.ldexpf
libc.src.__support.math.ldexpf128
libc.src.__support.math.ldexpf16
+ libc.src.__support.math.lgammaf16
libc.src.__support.math.llogbf
libc.src.__support.math.llogbf128
libc.src.__support.math.llogbf16
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 634778380dc8e..61f9bd547a001 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -58,6 +58,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::shared::log2f16(2.0f16));
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::shared::log2p1f16(1.0f16));
EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::shared::logbf16(1.0f16));
+ EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::shared::lgammaf16(1.0f16));
EXPECT_EQ(0L, LIBC_NAMESPACE::shared::llogbf16(1.0f16));
EXPECT_FP_EQ(0x1.921fb6p+0f16, LIBC_NAMESPACE::shared::acosf16(0.0f16));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 4213c11eca515..3278cf545c501 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -243,6 +243,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ lgammaf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ lgammaf16_test.cpp
+ DEPENDS
+ libc.src.math.lgammaf16
+)
+
add_fp_unittest(
tanpif16_test
NEED_MPFR
diff --git a/libc/test/src/math/lgammaf16_test.cpp b/libc/test/src/math/lgammaf16_test.cpp
new file mode 100644
index 0000000000000..5713b6449f276
--- /dev/null
+++ b/libc/test/src/math/lgammaf16_test.cpp
@@ -0,0 +1,42 @@
+//===-- Unittests for lgammaf16 -------------------------------------------===//
+//
+// 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/lgammaf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLgammaf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: ]0, +Inf[
+static constexpr uint16_t POS_START = 0x0001U;
+static constexpr uint16_t POS_STOP = 0x7bffU;
+
+// Range: ]-Inf, 0[
+static constexpr uint16_t NEG_START = 0x8001U;
+static constexpr uint16_t NEG_STOP = 0xfbffU;
+
+TEST_F(LlvmLibcLgammaf16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammaf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcLgammaf16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ if (FPBits(v).is_nan())
+ continue;
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammaf16(x), 0.5);
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 28b85b1a25bbd..3871798d00b7f 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -145,6 +145,17 @@ add_fp_unittest(
libc.src.math.tanpif
)
+add_fp_unittest(
+ lgammaf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ lgammaf16_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.src.math.lgammaf16
+)
+
add_fp_unittest(
tanpif16_test
SUITE
diff --git a/libc/test/src/math/smoke/lgammaf16_test.cpp b/libc/test/src/math/smoke/lgammaf16_test.cpp
new file mode 100644
index 0000000000000..8b302421b5ca5
--- /dev/null
+++ b/libc/test/src/math/smoke/lgammaf16_test.cpp
@@ -0,0 +1,61 @@
+//===-- Unittests for lgammaf16 -------------------------------------------===//
+//
+// 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 "src/math/lgammaf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcLgammaf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcLgammaf16Test, SpecialNumbers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::lgammaf16(sNaN),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::lgammaf16(aNaN));
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::lgammaf16(inf));
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::lgammaf16(neg_inf));
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf16(zero),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf16(neg_zero),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST_F(LlvmLibcLgammaf16Test, NegativeIntegers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf16(-1.0f16),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf16(-2.0f16),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf16(-3.0f16),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST_F(LlvmLibcLgammaf16Test, ExactValues) {
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::lgammaf16(1.0f16));
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::lgammaf16(2.0f16));
+}
+
+TEST_F(LlvmLibcLgammaf16Test, Overflow) {
+ // lgamma(8192) = 65623 > 65520 so it
+ // overflows to +Inf in round-to-nearest mode.
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf16(8192.0f16),
+ FE_OVERFLOW | FE_INEXACT);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index bdd8286148641..af69bd1421f95 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -323,6 +323,13 @@ MPFRNumber MPFRNumber::hypot(const MPFRNumber &b) {
return result;
}
+MPFRNumber MPFRNumber::lgamma() const {
+ MPFRNumber result(*this);
+ int signp;
+ mpfr_lgamma(result.value, &signp, value, mpfr_rounding);
+ return result;
+}
+
MPFRNumber MPFRNumber::log() const {
MPFRNumber result(*this);
mpfr_log(result.value, value, mpfr_rounding);
diff --git a/libc/utils/MPFRWrapper/MPCommon.h b/libc/utils/MPFRWrapper/MPCommon.h
index 935a2614968a2..38fd15fcc956c 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -212,6 +212,7 @@ class MPFRNumber {
MPFRNumber fmod(const MPFRNumber &b);
MPFRNumber frexp(int &exp);
MPFRNumber hypot(const MPFRNumber &b);
+ MPFRNumber lgamma() const;
MPFRNumber log() const;
MPFRNumber log2() const;
MPFRNumber log2p1() const;
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 356764302bda8..d585baa2e0d2c 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -75,6 +75,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.expm1();
case Operation::Floor:
return mpfrInput.floor();
+ case Operation::Lgamma:
+ return mpfrInput.lgamma();
case Operation::Log:
return mpfrInput.log();
case Operation::Log2:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index 6c1b15598b0f2..d490433685bcf 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -48,6 +48,7 @@ enum class Operation : int {
Exp10m1,
Expm1,
Floor,
+ Lgamma,
Log,
Log2,
Log2p1,
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index e6fd811ca50c2..2719832f7b268 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7594,6 +7594,24 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_lgammaf16",
+ hdrs = ["src/__support/math/lgammaf16.h"],
+ deps = [
+ ":__support_cpp_bit",
+ ":__support_fputil_cast",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer_operations",
+ ":__support_fputil_polyeval",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_types",
+ ":llvm_libc_macros_float16_macros",
+ ],
+)
+
libc_support_library(
name = "__support_math_llogbf16",
hdrs = ["src/__support/math/llogbf16.h"],
@@ -11279,6 +11297,11 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "lgammaf16",
+ additional_deps = [":__support_math_lgammaf16"],
+)
+
libc_math_function(
name = "llogb",
additional_deps = [
>From 954fb318700d56a59cbeba036f339aecd6e9fed2 Mon Sep 17 00:00:00 2001
From: Anonmiraj <ezzibrahimx at gmail.com>
Date: Thu, 7 May 2026 08:14:53 +0300
Subject: [PATCH 2/3] fix nits
---
libc/config/darwin/aarch64/entrypoints.txt | 1 +
libc/config/gpu/amdgpu/entrypoints.txt | 1 +
libc/config/gpu/nvptx/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 2 +-
libc/config/linux/riscv/entrypoints.txt | 2 +-
libc/config/linux/x86_64/entrypoints.txt | 2 +-
libc/shared/math/lgammaf16.h | 2 ++
libc/src/__support/math/lgammaf16.h | 12 ++++--------
8 files changed, 12 insertions(+), 11 deletions(-)
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index 914d2b7918da1..cdd04d8bbb74f 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -476,6 +476,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index c76900a8371f7..c1b29ff893891 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -597,6 +597,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index 6538732f785f5..73e2e62c9f3fa 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -599,6 +599,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 9bb3cfa06923f..2f52540adad92 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -595,7 +595,6 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
- libc.src.math.lgammaf16
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
@@ -770,6 +769,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 89b931faf29c2..12a44004cbe01 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -604,7 +604,6 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
- libc.src.math.lgammaf16
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
@@ -786,6 +785,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index b9e73d95ec244..801b0f55ccf7c 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -664,7 +664,6 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
- libc.src.math.lgammaf16
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
@@ -850,6 +849,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.iscanonicalf16
libc.src.math.issignalingf16
libc.src.math.ldexpf16
+ libc.src.math.lgammaf16
libc.src.math.llogbf16
libc.src.math.llrintf16
libc.src.math.llroundf16
diff --git a/libc/shared/math/lgammaf16.h b/libc/shared/math/lgammaf16.h
index 4ae8d36e5ba53..07de5ef6dd13d 100644
--- a/libc/shared/math/lgammaf16.h
+++ b/libc/shared/math/lgammaf16.h
@@ -9,6 +9,8 @@
#ifndef LLVM_LIBC_SHARED_MATH_LGAMMAF16_H
#define LLVM_LIBC_SHARED_MATH_LGAMMAF16_H
+#include "include/llvm-libc-macros/float16-macros.h"
+
#ifdef LIBC_TYPES_HAS_FLOAT16
#include "shared/libc_common.h"
diff --git a/libc/src/__support/math/lgammaf16.h b/libc/src/__support/math/lgammaf16.h
index 2e6424fc4826c..d93b801a4a38c 100644
--- a/libc/src/__support/math/lgammaf16.h
+++ b/libc/src/__support/math/lgammaf16.h
@@ -49,14 +49,10 @@ LIBC_INLINE constexpr bool is_integer(float16 x) {
// > P = fpminimax(sin(pi*(x+0.5))/(pi*(0.25-x^2)),
// [|0,2,4,6,8,10,12,14|], [|D...|], [-0.5, 0.5]);
LIBC_INLINE constexpr double lg_sinpi(double x) {
- constexpr double COEFFS[8] = {0x1p+2,
- -0x1.de9e64df22ea4p+1,
- 0x1.472be122401f8p+0,
- -0x1.d4fcd82df91bp-3,
- 0x1.9f05c97e0aab2p-6,
- -0x1.f3091c427b611p-10,
- 0x1.b22c9bfdca547p-14,
- -0x1.15484325ef569p-18};
+ constexpr double COEFFS[8] = {0x000000000000001p+2, -0x1.de9e64df22ea4p+1,
+ 0x1.472be122401f8p+0, -0x1.d4fcd82df91bp-3,
+ 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10,
+ 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18};
double u = x - 0.5;
double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
double p01 = fputil::multiply_add(u2, COEFFS[1], COEFFS[0]);
>From 0327facae226cf55b770c54df08224d951051ce3 Mon Sep 17 00:00:00 2001
From: Anonmiraj <ezzibrahimx at gmail.com>
Date: Wed, 24 Jun 2026 06:48:49 +0300
Subject: [PATCH 3/3] [libc][math] Implement single precision lgamma function
---
libc/config/baremetal/aarch64/entrypoints.txt | 1 +
libc/config/baremetal/arm/entrypoints.txt | 1 +
libc/config/baremetal/riscv/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/shared/math.h | 1 +
libc/shared/math/lgammaf.h | 23 ++
libc/src/__support/math/CMakeLists.txt | 34 +-
libc/src/__support/math/gamma_util.h | 116 ++++++
libc/src/__support/math/lgammaf.h | 374 ++++++++++++++++++
libc/src/__support/math/lgammaf16.h | 91 +----
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/generic/CMakeLists.txt | 13 +
libc/src/math/generic/lgammaf.cpp | 18 +
libc/src/math/lgammaf.h | 20 +
libc/test/src/math/CMakeLists.txt | 12 +
libc/test/src/math/exhaustive/CMakeLists.txt | 16 +
.../test/src/math/exhaustive/lgammaf_test.cpp | 35 ++
libc/test/src/math/lgammaf_test.cpp | 100 +++++
libc/test/src/math/smoke/CMakeLists.txt | 11 +
libc/test/src/math/smoke/lgammaf_test.cpp | 62 +++
.../llvm-project-overlay/libc/BUILD.bazel | 37 +-
23 files changed, 879 insertions(+), 91 deletions(-)
create mode 100644 libc/shared/math/lgammaf.h
create mode 100644 libc/src/__support/math/gamma_util.h
create mode 100644 libc/src/__support/math/lgammaf.h
create mode 100644 libc/src/math/generic/lgammaf.cpp
create mode 100644 libc/src/math/lgammaf.h
create mode 100644 libc/test/src/math/exhaustive/lgammaf_test.cpp
create mode 100644 libc/test/src/math/lgammaf_test.cpp
create mode 100644 libc/test/src/math/smoke/lgammaf_test.cpp
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index c2d9e860388f9..3f05a2f0ed004 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -478,6 +478,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 6c80f317917e2..bc70b5797af84 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -490,6 +490,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index 3dfa9013e6623..33fc0591bc84c 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -486,6 +486,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 2f52540adad92..64d87df1dc32f 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -595,6 +595,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 12a44004cbe01..df2441dff5056 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -604,6 +604,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 801b0f55ccf7c..8d49c2674de9c 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -664,6 +664,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.ldexp
libc.src.math.ldexpf
libc.src.math.ldexpl
+ libc.src.math.lgammaf
libc.src.math.llogb
libc.src.math.llogbf
libc.src.math.llogbl
diff --git a/libc/shared/math.h b/libc/shared/math.h
index b9060dc5f6245..cc95ab2b74234 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -305,6 +305,7 @@
#include "math/ldexpf128.h"
#include "math/ldexpf16.h"
#include "math/ldexpl.h"
+#include "math/lgammaf.h"
#include "math/lgammaf16.h"
#include "math/llogb.h"
#include "math/llogbbf16.h"
diff --git a/libc/shared/math/lgammaf.h b/libc/shared/math/lgammaf.h
new file mode 100644
index 0000000000000..6976109fe9832
--- /dev/null
+++ b/libc/shared/math/lgammaf.h
@@ -0,0 +1,23 @@
+//===-- Shared lgammaf 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_LGAMMAF_H
+#define LLVM_LIBC_SHARED_MATH_LGAMMAF_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/lgammaf.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::lgammaf;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_LGAMMAF_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index a7ff2cec529b5..b568ef151f85e 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3713,6 +3713,17 @@ add_header_library(
libc.src.__support.FPUtil.generic.sqrt
)
+add_header_library(
+ gamma_util
+ HDRS
+ gamma_util.h
+ DEPENDS
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.macros.config
+)
+
add_header_library(
getpayload
HDRS
@@ -3881,13 +3892,33 @@ add_header_library(
libc.include.llvm-libc-types.float128
)
+add_header_library(
+ lgammaf
+ HDRS
+ lgammaf.h
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.except_value_utils
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer_operations
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.math.gamma_util
+ libc.src.__support.math.log
+)
+
add_header_library(
lgammaf16
HDRS
lgammaf16.h
DEPENDS
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
@@ -3897,6 +3928,7 @@ add_header_library(
libc.src.__support.macros.config
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.types
+ libc.src.__support.math.gamma_util
)
add_header_library(
diff --git a/libc/src/__support/math/gamma_util.h b/libc/src/__support/math/gamma_util.h
new file mode 100644
index 0000000000000..09024a521ff77
--- /dev/null
+++ b/libc/src/__support/math/gamma_util.h
@@ -0,0 +1,116 @@
+//===-- Shared helpers for gamma family math functions ----------*- 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_GAMMA_UTIL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_GAMMA_UTIL_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace gamma_internal {
+
+template <typename T> LIBC_INLINE constexpr bool is_integer(T x) {
+ using FPBits = fputil::FPBits<T>;
+ using StorageType = typename FPBits::StorageType;
+ FPBits xbits(x);
+ StorageType x_u = xbits.uintval();
+ unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
+ unsigned lsb = static_cast<unsigned>(
+ cpp::countr_zero(static_cast<StorageType>(x_u | FPBits::EXP_MASK)));
+ constexpr unsigned UNIT_EXPONENT =
+ static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
+ return x_e + lsb >= UNIT_EXPONENT;
+}
+
+// sin(pi*x) for x in (0, 1), written as (0.25 - u^2) * P(u^2) where u = x-0.5.
+// Coefficients for P are a degree-7 polynomial in u^2
+// approximating sin(pi*(u+0.5)) / (pi*(0.25-u^2))
+// with max error ~2^{-55}, generated by Sollya with:
+// > P = fpminimax(sin(pi*(x+0.5))/(pi*(0.25-x^2)),
+// [|0,2,4,6,8,10,12,14|], [|D...|], [-0.5, 0.5]);
+LIBC_INLINE constexpr double lg_sinpi(double x) {
+ constexpr double COEFFS[8] = {0x000000000000001p+2, -0x1.de9e64df22ea4p+1,
+ 0x1.472be122401f8p+0, -0x1.d4fcd82df91bp-3,
+ 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10,
+ 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18};
+ double u = x - 0.5;
+ double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
+ double p01 = fputil::multiply_add(u2, COEFFS[1], COEFFS[0]);
+ double p23 = fputil::multiply_add(u2, COEFFS[3], COEFFS[2]);
+ double p45 = fputil::multiply_add(u2, COEFFS[5], COEFFS[4]);
+ double p67 = fputil::multiply_add(u2, COEFFS[7], COEFFS[6]);
+ double p03 = fputil::multiply_add(u4, p23, p01);
+ double p47 = fputil::multiply_add(u4, p67, p45);
+ // Compute (0.25 - u^2) = (0.5 - u) * (0.5 + u) to avoid ~10-digit
+ // catastrophic cancellation when |u| ~ 0.5 (i.e. x near 0 or 1).
+ return (0.5 - u) * (0.5 + u) * fputil::multiply_add(u8, p47, p03);
+}
+
+// Natural logarithm of x (x > 0), using 16-entry table + degree-7 polynomial.
+// Range reduction: x = 2^e * m where m in [1, 2), decomposed as
+// m = (1+i/16)*(1+z) so log(x) = e*log(2) + log(1+i/16) + log(1+z)
+// = e*log(2) + IL[i] + z*P(z).
+// P approximates log(1+z)/z on z in [-1/16, 1/16] with max
+// error ~2^{-54}, generated by Sollya with:
+// > P = fpminimax(log(1+x)/x, [|0,1,2,3,4,5,6,7|], [|D...|], [-1/16, 1/16]);
+LIBC_INLINE double lg_ln(double x) {
+ using FPBits = fputil::FPBits<double>;
+ uint64_t u = FPBits(x).uintval();
+ int e = static_cast<int>(FPBits(x).get_biased_exponent()) - 0x3ff;
+
+ // Coefficients for log(1 + z)/z on z in [-1/16, 1/16]
+ constexpr double COEFFS[8] = {0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2,
+ 0x1.55555537802dep-2, -0x1.ffffeca81b866p-3,
+ 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3,
+ 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4};
+ // IL[i] = log(1 + i/16) for i = 0..15
+ constexpr double IL[16] = {
+ 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4,
+ 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2,
+ 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2,
+ 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1,
+ 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1,
+ 0x1.52a2d265bc5abp-1};
+ // IX[i] = 1 / (1 + i/16) for i = 0..15
+ constexpr double IX[16] = {
+ 0x000000000000001p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1,
+ 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1,
+ 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1,
+ 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1,
+ 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1,
+ 0x1.0842108421084p-1};
+
+ int i = static_cast<int>((u >> 48) & 0xf);
+ // Reduce to mantissa in [1, 2)
+ uint64_t mant_u = (u & (~uint64_t(0) >> 12)) | (uint64_t(0x3ff) << 52);
+ double mant = FPBits(mant_u).get_val();
+ double z = IX[i] * mant - 1.0, z2 = z * z, z4 = z2 * z2;
+ double q01 = fputil::multiply_add(z, COEFFS[1], COEFFS[0]);
+ double q23 = fputil::multiply_add(z, COEFFS[3], COEFFS[2]);
+ double q45 = fputil::multiply_add(z, COEFFS[5], COEFFS[4]);
+ double q67 = fputil::multiply_add(z, COEFFS[7], COEFFS[6]);
+ double q03 = fputil::multiply_add(z2, q23, q01);
+ double q47 = fputil::multiply_add(z2, q67, q45);
+ return e * 0x1.62e42fefa39efp-1 + IL[i] +
+ z * fputil::multiply_add(z4, q47, q03);
+}
+
+} // namespace gamma_internal
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_GAMMA_UTIL_H
diff --git a/libc/src/__support/math/lgammaf.h b/libc/src/__support/math/lgammaf.h
new file mode 100644
index 0000000000000..9e4f4825d15c0
--- /dev/null
+++ b/libc/src/__support/math/lgammaf.h
@@ -0,0 +1,374 @@
+//===-- Implementation header for lgammaf -----------------------*- 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_LGAMMAF_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF_H
+
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/NearestIntegerOperations.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/gamma_util.h"
+#include "src/__support/math/log.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+LIBC_INLINE float lgammaf(float x) {
+ using namespace gamma_internal;
+ using FPBits = fputil::FPBits<float>;
+
+ FPBits xbits(x);
+ uint32_t x_abs = xbits.abs().uintval();
+
+ // NaN / Inf
+ if (LIBC_UNLIKELY(x_abs >= 0x7f800000u)) {
+ if (x_abs == 0x7f800000u)
+ return FPBits::inf().get_val();
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ // +/- 0 -> +Inf pole
+ if (LIBC_UNLIKELY(x_abs == 0)) {
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ fputil::set_errno_if_required(ERANGE);
+ return FPBits::inf().get_val();
+ }
+
+ // Negative integers and lgamma(1) = lgamma(2) = 0.
+ if (LIBC_UNLIKELY(is_integer(x))) {
+ if (xbits.is_neg()) {
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ fputil::set_errno_if_required(ERANGE);
+ return FPBits::inf().get_val();
+ }
+ if (x_abs == 0x3f800000u || x_abs == 0x40000000u)
+ return FPBits::zero().get_val();
+ }
+
+ constexpr fputil::ExceptValues<float, 11> LGAMMAF_EXCEPTS{{
+ // input, toward-zero result, RU, RD, RN
+ {0x3b7c53aau, 0x40b1d661u, 1, 0, 0},
+ {0x42468b59u, 0x430f25a7u, 1, 0, 0},
+ {0x50522f52u, 0x5292ee5du, 1, 0, 1},
+ {0x9b7679ffu, 0x4247c72cu, 1, 0, 1},
+ {0x9e88452du, 0x4236bd8bu, 1, 0, 1},
+ {0xa77a8e47u, 0x42052b94u, 1, 0, 1},
+ {0xb0e17820u, 0x41a1d37bu, 1, 0, 0},
+ {0xc02c060fu, 0xbdabc3e2u, 0, 1, 1},
+ {0xc134eb14u, 0xc1875615u, 0, 1, 0},
+ {0xc33139a3u, 0xc43991afu, 0, 1, 0},
+ {0xc6f7e151u, 0xc89116deu, 0, 1, 0},
+ }};
+ if (auto r = LGAMMAF_EXCEPTS.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
+ double xd = fputil::cast<double>(x);
+ double abs_xd = xd < 0.0 ? -xd : xd;
+ double lgamma_val;
+
+ // For very tiny |x| (< 2^-23), use the truncated Laurent series:
+ // lgamma(x) = -log(x) - gamma*x + O(x^2) for tiny x > 0
+ // lgamma(-y) = -log(y) + gamma*y + O(y^2) for tiny y > 0
+ // Even though gamma*|x| < 2^-25 is below 1 float ULP of |result|, it tips
+ // rounding at boundary cases. The x^2 term is at most gamma*2^-46 << 2^-50.
+ if (x_abs < 0x34000000u) { // |x| < 2^-23
+ constexpr double EULER_GAMMA = 0x1.2788cfc6fb619p-1;
+ double sign_corr = xbits.is_neg() ? EULER_GAMMA : -EULER_GAMMA;
+ return fputil::cast<float>(
+ fputil::multiply_add(sign_corr, abs_xd, -math::log(abs_xd)));
+ }
+
+ if (x_abs < 0x3f290000u) {
+ // Small: t = |x| < 0.66015625. Degree-22 Chebyshev (23 coeffs), Clenshaw.
+ // P(u) approximates h(t) = (lgamma(t) + log(t)) / t (smooth on [0, 0.66]).
+ constexpr double MID_S = 0x1.5200000000000p-2;
+ constexpr double HW_S = 0x1.5200000000000p-2;
+ constexpr double CHEB_S[23] = {
+ -0x1.6ab12743e97f5p-2, 0x1.ac5810962c7c9p-3, -0x1.9fc950c97fe08p-7,
+ 0x1.182a0ad81ff45p-10, -0x1.ae37b402356cbp-14, 0x1.617d9d1a7276fp-17,
+ -0x1.2e73d84dd6c00p-20, 0x1.09bfc3957e271p-23, -0x1.dbcc653579713p-27,
+ 0x1.afdf73641ecc6p-30, -0x1.8c3941844323bp-33, 0x1.6ea57fc6795d5p-36,
+ -0x1.55ad18a495c31p-39, 0x1.404bbf9a485bcp-42, -0x1.2dc653db9bc94p-45,
+ 0x1.1d92560a2236bp-48, -0x1.0f46c20905858p-51, 0x1.02925450b709dp-54,
+ -0x1.ee6a24e4dcbd9p-58, 0x1.d9f7b9780e3c5p-61, -0x1.c77a4c49dbb90p-64,
+ 0x1.b6972ac7c1946p-67, -0x1.a12b8ba3e9c77p-70};
+ double u = (abs_xd - MID_S) / HW_S;
+ double u2 = 2.0 * u;
+ double b1 = CHEB_S[22], b2 = 0.0;
+ for (int k = 21; k >= 1; --k) {
+ double b0 = fputil::multiply_add(u2, b1, CHEB_S[k] - b2);
+ b2 = b1;
+ b1 = b0;
+ }
+ double poly_h = fputil::multiply_add(u, b1, CHEB_S[0] - b2);
+ if (xbits.is_neg()) {
+ // log(pi) - log(sin(pi*x)) - poly_val, fused via FMA to cut roundings.
+ double neg_log_sin = -math::log(lg_sinpi(abs_xd));
+ lgamma_val = fputil::multiply_add(-abs_xd, poly_h,
+ 0x1.250d048e7a1bdp+0 + neg_log_sin);
+ } else {
+ // poly_val - log(abs_xd) with single rounding via FMA.
+ lgamma_val = fputil::multiply_add(abs_xd, poly_h, -math::log(abs_xd));
+ }
+ } else if (x_abs < 0x3f800000u) {
+ // M1: t in [0.66015625, 1.0). lgamma(t) = (t-1) * P_M1(u). Degree 22.
+ constexpr double MID_M1 = 0x1.a900000000000p-1;
+ constexpr double HW_M1 = 0x1.5c00000000000p-3;
+ constexpr double CHEB_M1[23] = {
+ -0x1.7a316a5bdc767p-1, 0x1.5b2bb19b764c1p-3, -0x1.1b0b304434f9bp-7,
+ 0x1.3ed224837f481p-11, -0x1.96032c59a8abfp-15, 0x1.1312af063c88ep-18,
+ -0x1.82beb75121dd4p-22, 0x1.16a1dfdaac960p-25, -0x1.98a1140d9aeeap-29,
+ 0x1.2faea3c424724p-32, -0x1.c8307861441f9p-36, 0x1.599dab4beb82dp-39,
+ -0x1.07bc40924e4c5p-42, 0x1.94f8259fea698p-46, -0x1.388d3c796c827p-49,
+ 0x1.e4a6b23ad71bcp-53, -0x1.7942bc85950eap-56, 0x1.26b4a27f3a966p-59,
+ -0x1.cde22532dd461p-63, 0x1.6af8963399134p-66, -0x1.1df7a1df10080p-69,
+ 0x1.c398dfa01e0e3p-73, -0x1.61e18974cb9b3p-76};
+ double u = (abs_xd - MID_M1) / HW_M1;
+ double u2 = 2.0 * u;
+ double b1 = CHEB_M1[22], b2 = 0.0;
+ for (int k = 21; k >= 1; --k) {
+ double b0 = fputil::multiply_add(u2, b1, CHEB_M1[k] - b2);
+ b2 = b1;
+ b1 = b0;
+ }
+ double poly = fputil::multiply_add(u, b1, CHEB_M1[0] - b2);
+ lgamma_val = (abs_xd - 1.0) * poly;
+ if (xbits.is_neg()) {
+ double frac_x = xd - fputil::floor(xd);
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd);
+ }
+ } else if (x_abs < 0x40000000u) {
+ // M2: t in [1.0, 2.0). lgamma(t) = (t-1)*(t-2) * P_M2(u). Degree 22 —
+ // hard cases near the lgamma minimum (~1.46) demand larger headroom.
+ constexpr double MID_M2 = 0x1.8000000000000p+0;
+ constexpr double HW_M2 = 0x1.0000000000000p-1;
+ constexpr double CHEB_M2[23] = {
+ 0x1.f73598c73fb27p-2, -0x1.37c37bb231109p-4, 0x1.144f77fa9d4e5p-7,
+ -0x1.1afb9900d693bp-10, 0x1.387dc81da3774p-13, -0x1.68eb0dd19f30ep-16,
+ 0x1.ad3824063f31fp-19, -0x1.047929385160cp-21, 0x1.40e5e09ac458ap-24,
+ -0x1.8fe38a574f8d3p-27, 0x1.f6dc9c81f711ap-30, -0x1.3e84fe8ed90fdp-32,
+ 0x1.96005b0ce149ap-35, -0x1.041c83c9d7060p-37, 0x1.4ecb51ecebfc7p-40,
+ -0x1.b09e091ec416cp-43, 0x1.187b631152304p-45, -0x1.6cd0f97db61f7p-48,
+ 0x1.dbd1c1a5b72f3p-51, -0x1.371247796da23p-53, 0x1.97a45af7c3ef7p-56,
+ -0x1.0b7547cbef874p-58, 0x1.5690ac55f96fcp-61};
+ double u = (abs_xd - MID_M2) / HW_M2;
+ double u2 = 2.0 * u;
+ double b1 = CHEB_M2[22], b2 = 0.0;
+ for (int k = 21; k >= 1; --k) {
+ double b0 = fputil::multiply_add(u2, b1, CHEB_M2[k] - b2);
+ b2 = b1;
+ b1 = b0;
+ }
+ double poly = fputil::multiply_add(u, b1, CHEB_M2[0] - b2);
+ lgamma_val = (abs_xd - 1.0) * (abs_xd - 2.0) * poly;
+ if (xbits.is_neg()) {
+ double frac_x = xd - fputil::floor(xd);
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd);
+ }
+ } else if (x_abs < 0x4057e000u) {
+ // M3: t in [2.0, 3.373046875). lgamma(t) = (t-2) * P_M3(u). Degree 22.
+ constexpr double MID_M3 = 0x1.57e0000000000p+1;
+ constexpr double HW_M3 = 0x1.5f80000000000p-1;
+ constexpr double CHEB_M3[23] = {
+ 0x1.377590a2d969ep-1, 0x1.66c0d826233d6p-3, -0x1.3868a8220c21ap-7,
+ 0x1.89b1790a524e9p-11, -0x1.22a7465347af0p-14, 0x1.d450a3bd37b6ep-18,
+ -0x1.8e6073ea06924p-21, 0x1.5f7f615a0f66ap-24, -0x1.3e43e37b75ba3p-27,
+ 0x1.25b86af37bf32p-30, -0x1.13066b68132aep-33, 0x1.0472e866f0a58p-36,
+ -0x1.f1c3c21853665p-40, 0x1.df2b520b24b05p-43, -0x1.d016652966750p-46,
+ 0x1.c3ca97d4de8a2p-49, -0x1.b9c00d6288552p-52, 0x1.b1919741c4e8ep-55,
+ -0x1.aaf26c5eb3ff4p-58, 0x1.a5a74c8e09064p-61, -0x1.a1817398c1204p-64,
+ 0x1.9e43a9d4a06e1p-67, -0x1.95b4706f5cbbdp-70};
+ double u = (abs_xd - MID_M3) / HW_M3;
+ double u2 = 2.0 * u;
+ double b1 = CHEB_M3[22], b2 = 0.0;
+ for (int k = 21; k >= 1; --k) {
+ double b0 = fputil::multiply_add(u2, b1, CHEB_M3[k] - b2);
+ b2 = b1;
+ b1 = b0;
+ }
+ double poly = fputil::multiply_add(u, b1, CHEB_M3[0] - b2);
+ lgamma_val = (abs_xd - 2.0) * poly;
+ if (xbits.is_neg()) {
+ // Near the regular lgamma zero at x ~= -2.7475: subtractive cancellation
+ // in the reflection formula kills precision. Use a Taylor expansion
+ // centered at the zero. Range bits in (0x402f95c2, 0x40301b93).
+ // Coefficients adopted from CORE-MATH (Sibidanov, 2023).
+ if (LIBC_UNLIKELY(x_abs > 0x402f95c2u && x_abs < 0x40301b93u)) {
+ double h = (xd + 0x1.5fb410a1bd901p+1) - 0x1.a19a96d2e6f85p-54;
+ constexpr double C[8] = {-0x1.ea12da904b18cp+0, 0x1.3267f3c265a54p+3,
+ -0x1.4185ac30cadb3p+4, 0x1.f504accc3f2e4p+5,
+ -0x1.8588444c679b4p+7, 0x1.43740491dc22p+9,
+ -0x1.12400ea23f9e6p+11, 0x1.dac829f365795p+12};
+ double h2 = h * h, h4 = h2 * h2;
+ double p01 = fputil::multiply_add(h, C[1], C[0]);
+ double p23 = fputil::multiply_add(h, C[3], C[2]);
+ double p45 = fputil::multiply_add(h, C[5], C[4]);
+ double p67 = fputil::multiply_add(h, C[7], C[6]);
+ double p03 = fputil::multiply_add(h2, p23, p01);
+ double p47 = fputil::multiply_add(h2, p67, p45);
+ lgamma_val = h * fputil::multiply_add(h4, p47, p03);
+ } else if (LIBC_UNLIKELY(x_abs > 0x401ceccbu && x_abs < 0x401d95cau)) {
+ // Near the regular lgamma zero at x ~= -2.3614: same issue
+ double h = (xd + 0x1.3a7fc9600f86cp+1) + 0x1.55f64f98af8dp-55;
+ constexpr double C[7] = {0x1.83fe966af535fp+0, 0x1.36eebb002f61ap+2,
+ 0x1.694a60589a0b3p+0, 0x1.1718d7aedb0b5p+3,
+ 0x1.733a045eca0d3p+2, 0x1.8d4297421205bp+4,
+ 0x1.7feea5fb29965p+4};
+ double h2 = h * h, h4 = h2 * h2;
+ double p01 = fputil::multiply_add(h, C[1], C[0]);
+ double p23 = fputil::multiply_add(h, C[3], C[2]);
+ double p45 = fputil::multiply_add(h, C[5], C[4]);
+ double p46 = fputil::multiply_add(h2, C[6], p45);
+ double p03 = fputil::multiply_add(h2, p23, p01);
+ lgamma_val = h * fputil::multiply_add(h4, p46, p03);
+ } else if (LIBC_UNLIKELY(x_abs > 0x40492009u && x_abs < 0x404940efu)) {
+ // Near the regular lgamma zero at x ~= -3.1431: same issue
+ double h = (xd + 0x1.9260dbc9e59afp+1) + 0x1.f717cd335a7b3p-53;
+ constexpr double C[7] = {0x1.f20a65f2fac55p+2, 0x1.9d4d297715105p+4,
+ 0x1.c1137124d5b21p+6, 0x1.267203d24de38p+9,
+ 0x1.99a63399a0b44p+11, 0x1.2941214faaf0cp+14,
+ 0x1.bb912c0c9cdd1p+16};
+ double h2 = h * h, h4 = h2 * h2;
+ double p01 = fputil::multiply_add(h, C[1], C[0]);
+ double p23 = fputil::multiply_add(h, C[3], C[2]);
+ double p45 = fputil::multiply_add(h, C[5], C[4]);
+ double p46 = fputil::multiply_add(h2, C[6], p45);
+ double p03 = fputil::multiply_add(h2, p23, p01);
+ lgamma_val = h * fputil::multiply_add(h4, p46, p03);
+ } else {
+ double frac_x = xd - fputil::floor(xd);
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd);
+ }
+ }
+ } else {
+ // Large: |x| >= 3.373046875. Stirling + Bernoulli correction.
+ // lgamma(x) = (x-0.5)*log(x) - x + log(2*pi)/2 + (1/x)*P(1/x^2)
+ // = (x-0.5)*(log(x)-1) + STIR_CONST + (1/x)*P(1/x^2)
+ // STIR_CONST = log(2*pi)/2 - 0.5.
+ //
+ // log(x) = e*log(2) + log(m) where x = 2^e * m, m in [1, 2)
+ // = e*LOG2_HI + (e*LOG2_LO + log(m))
+ // LOG2_HI has 44 sig bits, so e*LOG2_HI is exact for |e| < 512.
+ // For huge positive x, lgamma(x) overflows float. Use a linear
+ // approximation in double that maps to the correct Inf/max_normal.
+ if (LIBC_UNLIKELY(!xbits.is_neg() && x >= 0x1.895f1cp+121f)) {
+ double r = fputil::multiply_add(xd, 0x1.4d3398p+6, 0x1.10f35ep+103);
+ float result = fputil::cast<float>(r);
+ if (FPBits(result).is_inf()) {
+ fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+ fputil::set_errno_if_required(ERANGE);
+ }
+ return result;
+ }
+
+ constexpr double LOG2_HI = 0x1.62e42fefa3800p-1;
+ constexpr double LOG2_LO = 0x1.ef00000000000p-45;
+ using DPBits = fputil::FPBits<double>;
+ DPBits abs_bits(abs_xd);
+ int e = static_cast<int>(abs_bits.get_biased_exponent()) - 0x3ff;
+ // m in [1, 2)
+ DPBits m_bits = abs_bits;
+ m_bits.set_biased_exponent(0x3ff);
+ double m = m_bits.get_val();
+ double log_m = math::log(m);
+ double e_d = static_cast<double>(e);
+ double log_hi = e_d * LOG2_HI; // exact
+ double log_lo = fputil::multiply_add(e_d, LOG2_LO, log_m);
+ fputil::DoubleDouble log_x = fputil::exact_add(log_hi, log_lo);
+
+ double xm = abs_xd - 0.5;
+ // (xm) * log_x
+ fputil::DoubleDouble prod = fputil::exact_mult(xm, log_x.hi);
+ prod.lo = fputil::multiply_add(xm, log_x.lo, prod.lo);
+
+ // result_main = (xm * log_x) - xm + STIR_CONST
+ fputil::DoubleDouble diff = fputil::exact_add(prod.hi, -xm);
+ diff.lo += prod.lo;
+ lgamma_val = diff.hi + (diff.lo + 0x1.acfe390c97d69p-2);
+
+ double inv_x = 1.0 / abs_xd;
+ double inv_x2 = inv_x * inv_x;
+
+ if (x_abs > 0x4b989680u) {
+ // |x| > 2e7 -> just 1/(12x). Next Bernoulli term is < 2^-80.
+ lgamma_val += inv_x * 0x1.5555555555555p-4;
+ } else if (x_abs > 0x44fa0000u) {
+ // |x| > 2000 -> 2-term BERN2.
+ constexpr double BERN2[2] = {0x1.5555555555555p-4, -0x1.6c16bfb7c65a8p-9};
+ lgamma_val += inv_x * fputil::multiply_add(inv_x2, BERN2[1], BERN2[0]);
+ } else if (x_abs > 0x42920000u) {
+ // |x| > 73 -> 4-term BERN4.
+ constexpr double BERN4[4] = {0x1.5555555555555p-4, -0x1.6c16c16c15f75p-9,
+ 0x1.a01a00593b36fp-11,
+ -0x1.37e91273668efp-11};
+ double inv_x4 = inv_x2 * inv_x2;
+ double p01 = fputil::multiply_add(inv_x2, BERN4[1], BERN4[0]);
+ double p23 = fputil::multiply_add(inv_x2, BERN4[3], BERN4[2]);
+ lgamma_val += inv_x * fputil::multiply_add(inv_x4, p23, p01);
+ } else {
+ // |x| in (3.373, 73]: degree-9 Chebyshev poly in s = 1/x^2 for the
+ // function h(s) = stir_resid(1/sqrt(s)) * (1/sqrt(s)), evaluated by
+ // Clenshaw. Then correction = h(s) * (1/x). Max abs error ~5.5e-18.
+ constexpr double CHEB_B10_MID = 0x1.68c7762f8b34ep-5;
+ constexpr double CHEB_B10_HW = 0x1.673ded08c6ba2p-5;
+ constexpr double CHEB_B10[10] = {
+ 0x1.54d759afaa7e1p-4, -0x1.f2c7cca18bb6bp-14,
+ 0x1.7601881c38683p-21, -0x1.5a6f789402bbcp-27,
+ 0x1.198036ed27e9bp-32, -0x1.5459d5559def8p-37,
+ 0x1.15c4fb498d074p-41, -0x1.1f1b0949ee034p-45,
+ 0x1.6759cacff87f6p-49, -0x1.046d1e7698673p-52};
+ double u = (inv_x2 - CHEB_B10_MID) / CHEB_B10_HW;
+ double u2 = 2.0 * u;
+ double b1 = CHEB_B10[9], b2 = 0.0;
+ for (int k = 8; k >= 1; --k) {
+ double b0 = fputil::multiply_add(u2, b1, CHEB_B10[k] - b2);
+ b2 = b1;
+ b1 = b0;
+ }
+ double poly = fputil::multiply_add(u, b1, CHEB_B10[0] - b2);
+ lgamma_val += inv_x * poly;
+ }
+
+ if (xbits.is_neg()) {
+ // Reflection: lgamma(x) = log(pi) - lgamma(|x|) -
+ // log(|x|*|sin(pi*frac_x)|)
+ double frac_x = xd - fputil::floor(xd);
+ lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd);
+ }
+ }
+
+ float result = fputil::cast<float>(lgamma_val);
+ if (LIBC_UNLIKELY(FPBits(result).is_inf())) {
+ fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+ fputil::set_errno_if_required(ERANGE);
+ }
+ return result;
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LGAMMAF_H
diff --git a/libc/src/__support/math/lgammaf16.h b/libc/src/__support/math/lgammaf16.h
index d93b801a4a38c..4abb8bc9aa361 100644
--- a/libc/src/__support/math/lgammaf16.h
+++ b/libc/src/__support/math/lgammaf16.h
@@ -13,7 +13,6 @@
#ifdef LIBC_TYPES_HAS_FLOAT16
-#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/NearestIntegerOperations.h"
@@ -23,100 +22,14 @@
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/macros/properties/types.h"
+#include "src/__support/math/gamma_util.h"
namespace LIBC_NAMESPACE_DECL {
namespace math {
-namespace lgammaf16_internal {
-
-LIBC_INLINE constexpr bool is_integer(float16 x) {
- using FPBits = fputil::FPBits<float16>;
- FPBits xbits(x);
- uint16_t x_u = xbits.uintval();
- unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
- unsigned lsb = static_cast<unsigned>(
- cpp::countr_zero(static_cast<uint16_t>(x_u | FPBits::EXP_MASK)));
- constexpr unsigned UNIT_EXPONENT =
- static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
- return x_e + lsb >= UNIT_EXPONENT;
-}
-
-// sin(pi*x) for x in (0, 1), written as (0.25 - u^2) * P(u^2) where u = x-0.5.
-// Coefficients for P are a degree-7 polynomial in u^2
-// approximating sin(pi*(u+0.5)) / (pi*(0.25-u^2))
-// with max error ~2^{-55}, generated by Sollya with:
-// > P = fpminimax(sin(pi*(x+0.5))/(pi*(0.25-x^2)),
-// [|0,2,4,6,8,10,12,14|], [|D...|], [-0.5, 0.5]);
-LIBC_INLINE constexpr double lg_sinpi(double x) {
- constexpr double COEFFS[8] = {0x000000000000001p+2, -0x1.de9e64df22ea4p+1,
- 0x1.472be122401f8p+0, -0x1.d4fcd82df91bp-3,
- 0x1.9f05c97e0aab2p-6, -0x1.f3091c427b611p-10,
- 0x1.b22c9bfdca547p-14, -0x1.15484325ef569p-18};
- double u = x - 0.5;
- double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
- double p01 = fputil::multiply_add(u2, COEFFS[1], COEFFS[0]);
- double p23 = fputil::multiply_add(u2, COEFFS[3], COEFFS[2]);
- double p45 = fputil::multiply_add(u2, COEFFS[5], COEFFS[4]);
- double p67 = fputil::multiply_add(u2, COEFFS[7], COEFFS[6]);
- double p03 = fputil::multiply_add(u4, p23, p01);
- double p47 = fputil::multiply_add(u4, p67, p45);
- return (0.25 - u2) * fputil::multiply_add(u8, p47, p03);
-}
-
-// Natural logarithm of x (x > 0), using 16-entry table + degree-7 polynomial.
-// Range reduction: x = 2^e * m where m in [1, 2), decomposed as
-// m = (1+i/16)*(1+z) so log(x) = e*log(2) + log(1+i/16) + log(1+z)
-// = e*log(2) + IL[i] + z*P(z).
-// P approximates log(1+z)/z on z in [-1/16, 1/16] with max
-// error ~2^{-54}, generated by Sollya with:
-// > P = fpminimax(log(1+x)/x, [|0,1,2,3,4,5,6,7|], [|D...|], [-1/16, 1/16]);
-LIBC_INLINE double lg_ln(double x) {
- using FPBits = fputil::FPBits<double>;
- uint64_t u = FPBits(x).uintval();
- int e = static_cast<int>(FPBits(x).get_biased_exponent()) - 0x3ff;
-
- // Coefficients for log(1 + z)/z on z in [-1/16, 1/16]
- constexpr double COEFFS[8] = {0x1.fffffffffff24p-1, -0x1.ffffffffd1d67p-2,
- 0x1.55555537802dep-2, -0x1.ffffeca81b866p-3,
- 0x1.999611761d772p-3, -0x1.54f3e581b61bfp-3,
- 0x1.1e642b4cb5143p-3, -0x1.9115a5af1e1edp-4};
- // IL[i] = log(1 + i/16) for i = 0..15
- constexpr double IL[16] = {
- 0x1.59caeec280116p-57, 0x1.f0a30c01162aap-5, 0x1.e27076e2af2ebp-4,
- 0x1.5ff3070a793d6p-3, 0x1.c8ff7c79a9a2p-3, 0x1.1675cababa60fp-2,
- 0x1.4618bc21c5ec2p-2, 0x1.739d7f6bbd007p-2, 0x1.9f323ecbf984dp-2,
- 0x1.c8ff7c79a9a21p-2, 0x1.f128f5faf06ecp-2, 0x1.0be72e4252a83p-1,
- 0x1.1e85f5e7040d1p-1, 0x1.307d7334f10bep-1, 0x1.41d8fe84672afp-1,
- 0x1.52a2d265bc5abp-1};
- // IX[i] = 1 / (1 + i/16) for i = 0..15
- constexpr double IX[16] = {
- 0x000000000000001p+0, 0x1.e1e1e1e1e1e1ep-1, 0x1.c71c71c71c71cp-1,
- 0x1.af286bca1af28p-1, 0x1.999999999999ap-1, 0x1.8618618618618p-1,
- 0x1.745d1745d1746p-1, 0x1.642c8590b2164p-1, 0x1.5555555555555p-1,
- 0x1.47ae147ae147bp-1, 0x1.3b13b13b13b14p-1, 0x1.2f684bda12f68p-1,
- 0x1.2492492492492p-1, 0x1.1a7b9611a7b96p-1, 0x1.1111111111111p-1,
- 0x1.0842108421084p-1};
-
- int i = static_cast<int>((u >> 48) & 0xf);
- // Reduce to mantissa in [1, 2)
- uint64_t mant_u = (u & (~uint64_t(0) >> 12)) | (uint64_t(0x3ff) << 52);
- double mant = FPBits(mant_u).get_val();
- double z = IX[i] * mant - 1.0, z2 = z * z, z4 = z2 * z2;
- double q01 = fputil::multiply_add(z, COEFFS[1], COEFFS[0]);
- double q23 = fputil::multiply_add(z, COEFFS[3], COEFFS[2]);
- double q45 = fputil::multiply_add(z, COEFFS[5], COEFFS[4]);
- double q67 = fputil::multiply_add(z, COEFFS[7], COEFFS[6]);
- double q03 = fputil::multiply_add(z2, q23, q01);
- double q47 = fputil::multiply_add(z2, q67, q45);
- return e * 0x1.62e42fefa39efp-1 + IL[i] +
- z * fputil::multiply_add(z4, q47, q03);
-}
-
-} // namespace lgammaf16_internal
-
LIBC_INLINE float16 lgammaf16(float16 x) {
- using namespace lgammaf16_internal;
+ using namespace gamma_internal;
using FPBits = fputil::FPBits<float16>;
FPBits xbits(x);
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 8d2a981679980..965c08990386f 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -376,6 +376,7 @@ add_math_entrypoint_object(issignalingf16)
add_math_entrypoint_object(issignalingf128)
add_math_entrypoint_object(issignalingbf16)
+add_math_entrypoint_object(lgammaf)
add_math_entrypoint_object(lgammaf16)
add_math_entrypoint_object(llogb)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 49ca911ef3005..acd6cd8abd01d 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1567,6 +1567,19 @@ add_entrypoint_object(
libc.src.__support.math.ilogbbf16
)
+add_entrypoint_object(
+ lgammaf
+ SRCS
+ lgammaf.cpp
+ HDRS
+ ../lgammaf.h
+ DEPENDS
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.math.lgammaf
+ libc.src.errno.errno
+)
+
add_entrypoint_object(
lgammaf16
SRCS
diff --git a/libc/src/math/generic/lgammaf.cpp b/libc/src/math/generic/lgammaf.cpp
new file mode 100644
index 0000000000000..e823d9d7cdb46
--- /dev/null
+++ b/libc/src/math/generic/lgammaf.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of lgammaf 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/lgammaf.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/math/lgammaf.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(float, lgammaf, (float x)) { return math::lgammaf(x); }
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/lgammaf.h b/libc/src/math/lgammaf.h
new file mode 100644
index 0000000000000..8a8fdd90dfee6
--- /dev/null
+++ b/libc/src/math/lgammaf.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for lgammaf -----------------------*- 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_LGAMMAF_H
+#define LLVM_LIBC_SRC_MATH_LGAMMAF_H
+
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float lgammaf(float x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_LGAMMAF_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 3278cf545c501..0d86a8fe0d2fa 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -243,6 +243,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ lgammaf_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ lgammaf_test.cpp
+ DEPENDS
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.math.lgammaf
+)
+
add_fp_unittest(
lgammaf16_test
NEED_MPFR
diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index 2d1301d3a1e66..f471e880752b1 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -800,3 +800,19 @@ add_fp_unittest(
LINK_LIBRARIES
-lpthread
)
+
+add_fp_unittest(
+ lgammaf_test
+ NO_RUN_POSTBUILD
+ NEED_MPFR
+ SUITE
+ libc_math_exhaustive_tests
+ SRCS
+ lgammaf_test.cpp
+ DEPENDS
+ .exhaustive_test
+ libc.src.math.lgammaf
+ libc.src.__support.FPUtil.fp_bits
+ LINK_LIBRARIES
+ -lpthread
+)
diff --git a/libc/test/src/math/exhaustive/lgammaf_test.cpp b/libc/test/src/math/exhaustive/lgammaf_test.cpp
new file mode 100644
index 0000000000000..743d43a70abb1
--- /dev/null
+++ b/libc/test/src/math/exhaustive/lgammaf_test.cpp
@@ -0,0 +1,35 @@
+//===-- Exhaustive test for lgammaf ---------------------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/math/lgammaf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLgammafTest = LIBC_NAMESPACE::testing::FPTest<float>;
+
+using LlvmLibcLgammafExhaustiveTest =
+ LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Lgamma,
+ LIBC_NAMESPACE::lgammaf>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Positive range: [0, +Inf)
+static constexpr uint32_t POS_START = 0x0000'0000U;
+static constexpr uint32_t POS_STOP = 0x7f80'0000U;
+
+// Negative range: [-0, -Inf)
+static constexpr uint32_t NEG_START = 0x8000'0000U;
+static constexpr uint32_t NEG_STOP = 0xff80'0000U;
+
+TEST_F(LlvmLibcLgammafExhaustiveTest, PositiveRange) {
+ test_full_range_all_roundings(POS_START, POS_STOP);
+}
+
+TEST_F(LlvmLibcLgammafExhaustiveTest, NegativeRange) {
+ test_full_range_all_roundings(NEG_START, NEG_STOP);
+}
diff --git a/libc/test/src/math/lgammaf_test.cpp b/libc/test/src/math/lgammaf_test.cpp
new file mode 100644
index 0000000000000..a49e730ce77ec
--- /dev/null
+++ b/libc/test/src/math/lgammaf_test.cpp
@@ -0,0 +1,100 @@
+//===-- Unittests for lgammaf ---------------------------------------------===//
+//
+// 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/stdint_proxy.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/libc_errno.h"
+#include "src/math/lgammaf.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcLgammafTest = LIBC_NAMESPACE::testing::FPTest<float>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Hard-to-round cases collected from prior exhaustive tests.
+// Each is a midpoint-case where small precision drift in the implementation
+// tips the double to float rounding direction.
+TEST_F(LlvmLibcLgammafTest, HardCases) {
+ static constexpr uint32_t INPUTS[] = {
+ 0x77ac5740u, // 0x1.58ace8p+112 (huge positive, Stirling path)
+ 0x87acf970u, // -0x1.f59f2ep-16 (small negative, reflection path)
+ 0xc0afda0bu, // -0x1.5fb416p+1 (negative middle, reflection through M3)
+ 0x3fc0737cu, // 0x1.80e6f8p+0 (medium near lgamma minimum)
+ 0x3fc07be9u, // 0x1.80f7d2p+0 (same)
+ 0x3fd0b2bfu, // 0x1.a1657ep+0 (same)
+ 0xb3000824u, // -0x1.001048p-25 (tiny negative)
+ 0x30f00a14u, // 0x1.e01428p-30 (tiny positive)
+ 0x00800000u, // 0x1p-126 (min normal positive)
+ 0x80800000u, // -0x1p-126 (min normal negative)
+ 0x3f800000u, // 0x1p+0 (lgamma(1) = 0 exactly)
+ 0x40000000u, // 0x1p+1 (lgamma(2) = 0 exactly)
+ };
+ for (uint32_t v : INPUTS) {
+ float x = FPBits(v).get_val();
+ libc_errno = 0;
+ float result = LIBC_NAMESPACE::lgammaf(x);
+ if (FPBits(result).is_nan() || FPBits(result).is_inf() || libc_errno != 0)
+ continue;
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammaf(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcLgammafTest, PositiveRange) {
+ constexpr uint32_t COUNT = 100'000;
+ // [min normal, max_normal].
+ constexpr uint32_t POS_START = 0x0080'0000U;
+ constexpr uint32_t POS_STOP = 0x7f7f'ffffU;
+ constexpr uint32_t STEP = (POS_STOP - POS_START) / COUNT;
+ for (uint32_t i = 0, v = POS_START; i <= COUNT; ++i, v += STEP) {
+ float x = FPBits(v).get_val();
+ libc_errno = 0;
+ float result = LIBC_NAMESPACE::lgammaf(x);
+ if (FPBits(result).is_nan() || FPBits(result).is_inf() || libc_errno != 0)
+ continue;
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammaf(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcLgammafTest, NegativeRange) {
+ constexpr uint32_t COUNT = 100'000;
+ //-max_normal, -min normal].
+ constexpr uint32_t NEG_START = 0x8080'0000U;
+ constexpr uint32_t NEG_STOP = 0xff7f'ffffU;
+ constexpr uint32_t STEP = (NEG_STOP - NEG_START) / COUNT;
+ for (uint32_t i = 0, v = NEG_START; i <= COUNT; ++i, v += STEP) {
+ float x = FPBits(v).get_val();
+ if (FPBits(v).is_nan() || FPBits(v).is_inf())
+ continue;
+ libc_errno = 0;
+ float result = LIBC_NAMESPACE::lgammaf(x);
+ if (FPBits(result).is_nan() || FPBits(result).is_inf() || libc_errno != 0)
+ continue;
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammaf(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcLgammafTest, SmallRange) {
+ constexpr float LO = -4.0f;
+ constexpr float HI = 4.0f;
+ constexpr uint32_t COUNT = 10'000;
+ for (uint32_t i = 0; i <= COUNT; ++i) {
+ float x =
+ LO + (HI - LO) * static_cast<float>(i) / static_cast<float>(COUNT);
+ libc_errno = 0;
+ float result = LIBC_NAMESPACE::lgammaf(x);
+ if (FPBits(result).is_nan() || FPBits(result).is_inf() || libc_errno != 0)
+ continue;
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Lgamma, x,
+ LIBC_NAMESPACE::lgammaf(x), 0.5);
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 3871798d00b7f..56d785d8c6139 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -145,6 +145,17 @@ add_fp_unittest(
libc.src.math.tanpif
)
+add_fp_unittest(
+ lgammaf_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ lgammaf_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.src.math.lgammaf
+)
+
add_fp_unittest(
lgammaf16_test
SUITE
diff --git a/libc/test/src/math/smoke/lgammaf_test.cpp b/libc/test/src/math/smoke/lgammaf_test.cpp
new file mode 100644
index 0000000000000..e02e2cb045223
--- /dev/null
+++ b/libc/test/src/math/smoke/lgammaf_test.cpp
@@ -0,0 +1,62 @@
+//===-- Unittests for lgammaf ---------------------------------------------===//
+//
+// 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 "src/math/lgammaf.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcLgammafTest = LIBC_NAMESPACE::testing::FPTest<float>;
+
+TEST_F(LlvmLibcLgammafTest, SpecialNumbers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::lgammaf(sNaN), FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::lgammaf(aNaN));
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::lgammaf(inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::lgammaf(neg_inf));
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(zero), FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(neg_zero),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST_F(LlvmLibcLgammafTest, NegativeIntegers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(-1.0f),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(-2.0f),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(-100.0f),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ // Large negative integer (still representable as float).
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(-0x1p23f),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST_F(LlvmLibcLgammafTest, ExactValues) {
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::lgammaf(1.0f));
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::lgammaf(2.0f));
+}
+
+TEST_F(LlvmLibcLgammafTest, Overflow) {
+ // lgamma(x) overflows float around x ~ 2^121. Pick a comfortable margin.
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::lgammaf(0x1p126f),
+ FE_OVERFLOW | FE_INEXACT);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 2719832f7b268..3a6b20b0d7800 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7169,6 +7169,17 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_gamma_util",
+ hdrs = ["src/__support/math/gamma_util.h"],
+ deps = [
+ ":__support_cpp_bit",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_macros_config",
+ ],
+)
+
libc_support_library(
name = "__support_math_getpayload",
hdrs = ["src/__support/math/getpayload.h"],
@@ -7594,11 +7605,29 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_lgammaf",
+ hdrs = ["src/__support/math/lgammaf.h"],
+ deps = [
+ ":__support_fputil_cast",
+ ":__support_fputil_except_value_utils",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer_operations",
+ ":__support_fputil_polyeval",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_gamma_util",
+ ":errno_macros",
+ ":fenv_macros",
+ ],
+)
+
libc_support_library(
name = "__support_math_lgammaf16",
hdrs = ["src/__support/math/lgammaf16.h"],
deps = [
- ":__support_cpp_bit",
":__support_fputil_cast",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
@@ -7608,6 +7637,7 @@ libc_support_library(
":__support_macros_config",
":__support_macros_optimization",
":__support_macros_properties_types",
+ ":__support_math_gamma_util",
":llvm_libc_macros_float16_macros",
],
)
@@ -11297,6 +11327,11 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "lgammaf",
+ additional_deps = [":__support_math_lgammaf"],
+)
+
libc_math_function(
name = "lgammaf16",
additional_deps = [":__support_math_lgammaf16"],
More information about the libc-commits
mailing list