[libc-commits] [libc] [llvm] [libc][math] Implement single precision lgamma function (PR #205490)
via libc-commits
libc-commits at lists.llvm.org
Tue Jul 28 07:02:37 PDT 2026
https://github.com/AnonMiraj updated https://github.com/llvm/llvm-project/pull/205490
>From 71a9f4d2a1e03a0be3e5ccdbeb55e4539109bce6 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/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 | 92 +----
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(+), 92 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 93808e4c471fc..1113dd2c50ffc 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -483,6 +483,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 71ecbe438bf7b..2a40e4dee9efa 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -495,6 +495,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 b38f60ab9d5f5..0998df01aecf1 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -491,6 +491,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 2e7b9276ab068..0ed35e71938a3 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -622,6 +622,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 8a29e77ffe0fd..e3f4b54e25eea 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -697,6 +697,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 0a4fb747c2940..63121f0dbab6b 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -702,6 +702,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 b86fcbc0abcea..7b52cedf99c24 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/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 9f4624682eafe..543190db38e42 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3784,6 +3784,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
@@ -3943,13 +3954,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
@@ -3959,6 +3990,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..007786368b7d5
--- /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 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 34f2c2717f493..28660fe92d5ad 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);
@@ -156,7 +69,6 @@ LIBC_INLINE float16 lgammaf16(float16 x) {
double abs_xd = xd < 0.0 ? -xd : xd;
double lgamma_val;
- // |x| = 0.66015625
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
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 8c9a184ca517f..8ac6c9106ad40 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -381,6 +381,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 d5661f6f55e91..4d6af7219a067 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1587,6 +1587,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 240689c0de02f..97c1a03a358a4 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -256,6 +256,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 81f76fa7681fd..71af8729193ab 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -158,6 +158,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 85cd4feddb958..5ae45b8deff3d 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -7909,6 +7909,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"],
@@ -8323,11 +8334,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",
@@ -8337,6 +8366,7 @@ libc_support_library(
":__support_macros_config",
":__support_macros_optimization",
":__support_macros_properties_types",
+ ":__support_math_gamma_util",
":llvm_libc_macros_float16_macros",
],
)
@@ -12080,6 +12110,11 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "lgammaf",
+ additional_deps = [":__support_math_lgammaf"],
+)
+
libc_math_function(
name = "lgammaf16",
additional_deps = [":__support_math_lgammaf16"],
>From 16b974ba60bc950b9c8d49483d1bfa7a03eba149 Mon Sep 17 00:00:00 2001
From: Anonmiraj <ezzibrahimx at gmail.com>
Date: Mon, 20 Jul 2026 04:06:33 +0300
Subject: [PATCH 2/3] optimize all paths
---
libc/src/__support/math/lgammaf.h | 513 +++++++++++++++++++-----------
1 file changed, 323 insertions(+), 190 deletions(-)
diff --git a/libc/src/__support/math/lgammaf.h b/libc/src/__support/math/lgammaf.h
index 9e4f4825d15c0..e2f0a16aa4b2b 100644
--- a/libc/src/__support/math/lgammaf.h
+++ b/libc/src/__support/math/lgammaf.h
@@ -14,7 +14,6 @@
#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"
@@ -26,8 +25,46 @@ namespace LIBC_NAMESPACE_DECL {
namespace math {
+namespace lgammaf_internal {
+
+// P_M2(d), d = t - 1.5, approximates lgamma(t)/((t-1)(t-2)) on [1, 2].
+// Degree 18 centered-monomial coefficients from a 250-bit Chebyshev fit,
+// max relative error 2^-51.6
+LIBC_INLINE double lgamma_m2_poly(double d) {
+ constexpr double POLY_M2[19] = {
+ 0x1.eeb95b094c191p-2, -0x1.2aed059bd613dp-3, 0x1.01af62a292eb9p-4,
+ -0x1.007aa83cc458bp-5, 0x1.13342351db068p-6, -0x1.34dbda42244dep-7,
+ 0x1.64f066cf96646p-8, -0x1.a5033bab00972p-9, 0x1.f814855fdfeb8p-10,
+ -0x1.31418be9f2e14p-10, 0x1.7510b1d1fa536p-11, -0x1.caccf4395927bp-12,
+ 0x1.1c0d19e4fef05p-12, -0x1.67d43c4be68f7p-13, 0x1.c3d6281f9ac35p-14,
+ -0x1.dd9e09255e712p-15, 0x1.2550c8f33e895p-15, -0x1.6cd0f980284cbp-15,
+ 0x1.dbd1c21b8dfa7p-16};
+ double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
+ double p01 = fputil::multiply_add(d, POLY_M2[1], POLY_M2[0]);
+ double p23 = fputil::multiply_add(d, POLY_M2[3], POLY_M2[2]);
+ double p45 = fputil::multiply_add(d, POLY_M2[5], POLY_M2[4]);
+ double p67 = fputil::multiply_add(d, POLY_M2[7], POLY_M2[6]);
+ double p89 = fputil::multiply_add(d, POLY_M2[9], POLY_M2[8]);
+ double p1011 = fputil::multiply_add(d, POLY_M2[11], POLY_M2[10]);
+ double p1213 = fputil::multiply_add(d, POLY_M2[13], POLY_M2[12]);
+ double p1415 = fputil::multiply_add(d, POLY_M2[15], POLY_M2[14]);
+ double p1617 = fputil::multiply_add(d, POLY_M2[17], POLY_M2[16]);
+ double q03 = fputil::multiply_add(d2, p23, p01);
+ double q47 = fputil::multiply_add(d2, p67, p45);
+ double q811 = fputil::multiply_add(d2, p1011, p89);
+ double q1215 = fputil::multiply_add(d2, p1415, p1213);
+ double q1618 = fputil::multiply_add(d2, POLY_M2[18], p1617);
+ double r07 = fputil::multiply_add(d4, q47, q03);
+ double r815 = fputil::multiply_add(d4, q1215, q811);
+ double s015 = fputil::multiply_add(d8, r815, r07);
+ return fputil::multiply_add(d16, q1618, s015);
+}
+
+} // namespace lgammaf_internal
+
LIBC_INLINE float lgammaf(float x) {
using namespace gamma_internal;
+ using namespace lgammaf_internal;
using FPBits = fputil::FPBits<float>;
FPBits xbits(x);
@@ -62,24 +99,6 @@ LIBC_INLINE float lgammaf(float x) {
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;
@@ -90,6 +109,16 @@ LIBC_INLINE float lgammaf(float x) {
// 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 fputil::ExceptValues<float, 4> LGAMMAF_EXCEPTS_TINY{{
+ // input, toward-zero result, RU, RD, RN
+ {0x9b7679ffu, 0x4247c72cu, 1, 0, 1},
+ {0x9e88452du, 0x4236bd8bu, 1, 0, 1},
+ {0xa77a8e47u, 0x42052b94u, 1, 0, 1},
+ {0xb0e17820u, 0x41a1d37bu, 1, 0, 0},
+ }};
+ if (auto r = LGAMMAF_EXCEPTS_TINY.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
constexpr double EULER_GAMMA = 0x1.2788cfc6fb619p-1;
double sign_corr = xbits.is_neg() ? EULER_GAMMA : -EULER_GAMMA;
return fputil::cast<float>(
@@ -97,116 +126,164 @@ LIBC_INLINE float lgammaf(float x) {
}
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);
+ // Small negative: x in (-0.66015625, -2^-23). Degree-25 monomial fit
+ // of g(x) = (lgamma(x) + log(-x)) / x (smooth on [-0.66, 0]), then
+ // lgamma(x) = x*g(x) - log(-x). d = x - MID is exact.
+ constexpr double MID_SN = -0x1.5200000000000p-2;
+ constexpr double POLY_SN[26] = {
+ -0x1.cf99908bbb1d7p-1, 0x1.386ced0346de0p+0, -0x1.d0b4c2d337e0fp-1,
+ 0x1.db2dcaaa022a9p-1, -0x1.1188790a072bap+0, 0x1.4e50420966212p+0,
+ -0x1.a787465391482p+0, 0x1.12d6375bce19ep+1, -0x1.6b05729f76a9cp+1,
+ 0x1.e602a94d73479p+1, -0x1.48d5f84ddbdefp+2, 0x1.c0e7f7d05fe6dp+2,
+ -0x1.34d52338a0d93p+3, 0x1.ab7e04eff787dp+3, -0x1.2724d2ace1c81p+4,
+ 0x1.9b17bee5ae20ap+4, -0x1.34beabd4deff5p+5, 0x1.bc510ff62dcf0p+5,
+ -0x1.7406b7594b6ddp+5, 0x1.b97b7265194a1p+5, -0x1.5e2ba44fee700p+8,
+ 0x1.18310d7d43c46p+9, 0x1.59472c8ca4419p+9, -0x1.2a61147e20ef1p+10,
+ -0x1.7cfcf784946edp+11, 0x1.2436ec1277b00p+12};
+ double d = xd - MID_SN;
+ double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
+ double p01 = fputil::multiply_add(d, POLY_SN[1], POLY_SN[0]);
+ double p23 = fputil::multiply_add(d, POLY_SN[3], POLY_SN[2]);
+ double p45 = fputil::multiply_add(d, POLY_SN[5], POLY_SN[4]);
+ double p67 = fputil::multiply_add(d, POLY_SN[7], POLY_SN[6]);
+ double p89 = fputil::multiply_add(d, POLY_SN[9], POLY_SN[8]);
+ double p1011 = fputil::multiply_add(d, POLY_SN[11], POLY_SN[10]);
+ double p1213 = fputil::multiply_add(d, POLY_SN[13], POLY_SN[12]);
+ double p1415 = fputil::multiply_add(d, POLY_SN[15], POLY_SN[14]);
+ double p1617 = fputil::multiply_add(d, POLY_SN[17], POLY_SN[16]);
+ double p1819 = fputil::multiply_add(d, POLY_SN[19], POLY_SN[18]);
+ double p2021 = fputil::multiply_add(d, POLY_SN[21], POLY_SN[20]);
+ double p2223 = fputil::multiply_add(d, POLY_SN[23], POLY_SN[22]);
+ double p2425 = fputil::multiply_add(d, POLY_SN[25], POLY_SN[24]);
+ double q03 = fputil::multiply_add(d2, p23, p01);
+ double q47 = fputil::multiply_add(d2, p67, p45);
+ double q811 = fputil::multiply_add(d2, p1011, p89);
+ double q1215 = fputil::multiply_add(d2, p1415, p1213);
+ double q1619 = fputil::multiply_add(d2, p1819, p1617);
+ double q2023 = fputil::multiply_add(d2, p2223, p2021);
+ double r07 = fputil::multiply_add(d4, q47, q03);
+ double r815 = fputil::multiply_add(d4, q1215, q811);
+ double r1623 = fputil::multiply_add(d4, q2023, q1619);
+ double s015 = fputil::multiply_add(d8, r815, r07);
+ double s1625 = fputil::multiply_add(d8, p2425, r1623);
+ double poly_g = fputil::multiply_add(d16, s1625, s015);
+ // lgamma(x) = x*g(x) - log(-x) with single rounding via FMA.
+ lgamma_val = fputil::multiply_add(xd, poly_g, -math::log(abs_xd));
} else {
+ // x = 0x1.f8a754p-9f
+ if (LIBC_UNLIKELY(xbits.uintval() == 0x3b7c53aau))
+ return fputil::round_result_slightly_up(0x1.63acc2p+2f);
+ // Small: t = x < 0.66015625. Degree-17 monomial fit of
+ // h(t) = (lgamma(t) + log(t)) / t (smooth on [0, 0.66]).
+ // d = t - MID is exact.
+ constexpr double MID_S = 0x1.5200000000000p-2;
+ constexpr double POLY_S[18] = {
+ -0x1.5dcd7586bfd88p-2, 0x1.3f88851c787c3p-1, -0x1.cdfaca3081737p-3,
+ 0x1.cfac7b321198bp-4, -0x1.0939a239f89b2p-4, 0x1.44cddc5f52a43p-5,
+ -0x1.9e4cc54acdcebp-6, 0x1.0f61691962355p-6, -0x1.6a429a6c89071p-7,
+ 0x1.ea573d148429cp-8, -0x1.4f6ad0c206ae6p-8, 0x1.cedaf9fc19cb5p-9,
+ -0x1.42333c4e88f7dp-9, 0x1.c2a0e3ffc858fp-10, -0x1.32153d4cb9316p-10,
+ 0x1.add42f56d3c7ap-11, -0x1.972ff2b1a91edp-11, 0x1.25f56c66e3049p-11};
+ double d = abs_xd - MID_S;
+ double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
+ double p01 = fputil::multiply_add(d, POLY_S[1], POLY_S[0]);
+ double p23 = fputil::multiply_add(d, POLY_S[3], POLY_S[2]);
+ double p45 = fputil::multiply_add(d, POLY_S[5], POLY_S[4]);
+ double p67 = fputil::multiply_add(d, POLY_S[7], POLY_S[6]);
+ double p89 = fputil::multiply_add(d, POLY_S[9], POLY_S[8]);
+ double p1011 = fputil::multiply_add(d, POLY_S[11], POLY_S[10]);
+ double p1213 = fputil::multiply_add(d, POLY_S[13], POLY_S[12]);
+ double p1415 = fputil::multiply_add(d, POLY_S[15], POLY_S[14]);
+ double p1617 = fputil::multiply_add(d, POLY_S[17], POLY_S[16]);
+ double q03 = fputil::multiply_add(d2, p23, p01);
+ double q47 = fputil::multiply_add(d2, p67, p45);
+ double q811 = fputil::multiply_add(d2, p1011, p89);
+ double q1215 = fputil::multiply_add(d2, p1415, p1213);
+ double r07 = fputil::multiply_add(d4, q47, q03);
+ double r815 = fputil::multiply_add(d4, q1215, q811);
+ double s015 = fputil::multiply_add(d8, r815, r07);
+ double poly_h = fputil::multiply_add(d16, p1617, s015);
// 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);
+ // x in (-1, -0.66015625]: Gamma(x) = Gamma(x+2)/(x(x+1)), so
+ // lgamma(x) = x(x+1)*P_M2(x+0.5) - log(-x(x+1)). x+0.5 and x+1 are
+ // exact (Sterbenz); x(x+1) = (t-1)(t-2) is the M2 prefactor.
+ double w = xd * (xd + 1.0);
+ double poly = lgamma_m2_poly(xd + 0.5);
+ lgamma_val = fputil::multiply_add(w, poly, -lg_ln(-w));
+ } else {
+ // M1: t in [0.66015625, 1.0). lgamma(t) = (t-1) * P_M1(d), d = t - MID.
+ // Degree-14 monomial fit of lgamma(t)/(t-1), max error 2^-51.3.
+ constexpr double MID_M1 = 0x1.a900000000000p-1;
+ constexpr double POLY_M1[15] = {
+ -0x1.75cb89aad8dc7p-1, 0x1.f958114ec790dp-1, -0x1.2b968b172a05cp-1,
+ 0x1.eaf0e2c8dc49fp-2, -0x1.c6efb493bfadfp-2, 0x1.c0a621f0e593ep-2,
+ -0x1.cb2416349da49p-2, 0x1.e19328033256dp-2, -0x1.010d0c64743f5p-1,
+ 0x1.162e14c80bdc9p-1, -0x1.3034dec6eb9a3p-1, 0x1.4ca601da8405ep-1,
+ -0x1.70f73638718d4p-1, 0x1.ded7569ab5a36p-1, -0x1.0fdd39040adfbp+0};
+ double d = abs_xd - MID_M1;
+ double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4;
+ double p01 = fputil::multiply_add(d, POLY_M1[1], POLY_M1[0]);
+ double p23 = fputil::multiply_add(d, POLY_M1[3], POLY_M1[2]);
+ double p45 = fputil::multiply_add(d, POLY_M1[5], POLY_M1[4]);
+ double p67 = fputil::multiply_add(d, POLY_M1[7], POLY_M1[6]);
+ double p89 = fputil::multiply_add(d, POLY_M1[9], POLY_M1[8]);
+ double p1011 = fputil::multiply_add(d, POLY_M1[11], POLY_M1[10]);
+ double p1213 = fputil::multiply_add(d, POLY_M1[13], POLY_M1[12]);
+ double q03 = fputil::multiply_add(d2, p23, p01);
+ double q47 = fputil::multiply_add(d2, p67, p45);
+ double q811 = fputil::multiply_add(d2, p1011, p89);
+ double q1214 = fputil::multiply_add(d2, POLY_M1[14], p1213);
+ double r07 = fputil::multiply_add(d4, q47, q03);
+ double r814 = fputil::multiply_add(d4, q1214, q811);
+ double poly = fputil::multiply_add(d8, r814, r07);
+ lgamma_val = (abs_xd - 1.0) * poly;
}
} 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);
+ // x in (-2, -1): Gamma(x) = Gamma(x+3)/(x(x+1)(x+2)), so
+ // lgamma(x) = (x+1)(x+2)*P_M2(x+1.5) - log(x(x+1)(x+2)). The shifts
+ // are exact (Sterbenz); (x+1)(x+2) = (t-1)(t-2) is the M2 prefactor.
+ double u = (xd + 1.0) * (xd + 2.0);
+ double poly = lgamma_m2_poly(xd + 1.5);
+ lgamma_val = fputil::multiply_add(u, poly, -lg_ln(xd * u));
+ } else {
+ // M2: t in [1.0, 2.0). lgamma(t) = (t-1)*(t-2) * P_M2(t - 1.5).
+ double d = abs_xd - 0x1.8p+0;
+ lgamma_val = (abs_xd - 1.0) * (abs_xd - 2.0) * lgamma_m2_poly(d);
}
} else if (x_abs < 0x4057e000u) {
- // M3: t in [2.0, 3.373046875). lgamma(t) = (t-2) * P_M3(u). Degree 22.
+ // M3: t in [2.0, 3.373046875). lgamma(t) = (t-2) * P_M3(d), d = t - MID.
+ // Degree-15 monomial fit of lgamma(t)/(t-2), max error 2^-49.7.
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);
+ constexpr double POLY_M3[16] = {
+ 0x1.3c4e36a0b4775p-1, 0x1.01f945be1325fp-2, -0x1.4203a95730c77p-5,
+ 0x1.227f82db7c7a1p-7, -0x1.32f092b6ec5a3p-9, 0x1.61df821ae7829p-11,
+ -0x1.aed17eeca55e6p-13, 0x1.100a8e7dcd62cp-14, -0x1.609aa16f960a0p-16,
+ 0x1.d1e293fa801bbp-18, -0x1.38b436ce1b3b9p-19, 0x1.a838eec563338p-21,
+ -0x1.1a6a9cf4aee4bp-22, 0x1.8387c2a068b06p-24, -0x1.5ee6ae8c133f7p-25,
+ 0x1.f195cf3f4b24ep-27};
+ double d = abs_xd - MID_M3;
+ double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4;
+ double p01 = fputil::multiply_add(d, POLY_M3[1], POLY_M3[0]);
+ double p23 = fputil::multiply_add(d, POLY_M3[3], POLY_M3[2]);
+ double p45 = fputil::multiply_add(d, POLY_M3[5], POLY_M3[4]);
+ double p67 = fputil::multiply_add(d, POLY_M3[7], POLY_M3[6]);
+ double p89 = fputil::multiply_add(d, POLY_M3[9], POLY_M3[8]);
+ double p1011 = fputil::multiply_add(d, POLY_M3[11], POLY_M3[10]);
+ double p1213 = fputil::multiply_add(d, POLY_M3[13], POLY_M3[12]);
+ double p1415 = fputil::multiply_add(d, POLY_M3[15], POLY_M3[14]);
+ double q03 = fputil::multiply_add(d2, p23, p01);
+ double q47 = fputil::multiply_add(d2, p67, p45);
+ double q811 = fputil::multiply_add(d2, p1011, p89);
+ double q1215 = fputil::multiply_add(d2, p1415, p1213);
+ double r07 = fputil::multiply_add(d4, q47, q03);
+ double r815 = fputil::multiply_add(d4, q1215, q811);
+ double poly = fputil::multiply_add(d8, r815, r07);
lgamma_val = (abs_xd - 2.0) * poly;
if (xbits.is_neg()) {
// Near the regular lgamma zero at x ~= -2.7475: subtractive cancellation
@@ -256,9 +333,45 @@ LIBC_INLINE float lgammaf(float x) {
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);
+ // x in (-3.373, -2): Gamma(x) = Gamma(x+5)/(x(x+1)(x+2)(x+3)(x+4)),
+ // t = x+5 in (1.627, 3), so lgamma(x) = (x+3)*Q_M3N(d) - log|prod|
+ // with d = x + 2.6865234375 (exact). Q_M3N is a degree-17 monomial
+ // fit of lgamma(t)/(t-2) on [1.627, 3]; the (t-2) = x+3 prefactor
+ // removes the lgamma zero at t = 2.
+ constexpr double POLY_M3N[18] = {
+ 0x1.091ff92b41f07p-1, 0x1.245eed13c42b0p-2,
+ -0x1.a6c3cf8a165cfp-5, 0x1.bc8d19a3ade02p-7,
+ -0x1.1231651010470p-8, 0x1.710f57a2abe2bp-10,
+ -0x1.061ddab9ab12cp-11, 0x1.81edb077ed799p-13,
+ -0x1.235fbc13dbf04p-14, 0x1.c02e6f8fb5dbfp-16,
+ -0x1.5d75e6a94d352p-17, 0x1.1384596ce083dp-18,
+ -0x1.b8da81a716039p-20, 0x1.61ac57fe0e288p-21,
+ -0x1.077101d536ad8p-22, 0x1.a59e2356ba870p-24,
+ -0x1.0cfcf8f102a43p-24, 0x1.c141a51f827f1p-26};
+ double d = xd + 0x1.57ep+1;
+ double d2 = d * d, d4 = d2 * d2, d8 = d4 * d4, d16 = d8 * d8;
+ double p01 = fputil::multiply_add(d, POLY_M3N[1], POLY_M3N[0]);
+ double p23 = fputil::multiply_add(d, POLY_M3N[3], POLY_M3N[2]);
+ double p45 = fputil::multiply_add(d, POLY_M3N[5], POLY_M3N[4]);
+ double p67 = fputil::multiply_add(d, POLY_M3N[7], POLY_M3N[6]);
+ double p89 = fputil::multiply_add(d, POLY_M3N[9], POLY_M3N[8]);
+ double p1011 = fputil::multiply_add(d, POLY_M3N[11], POLY_M3N[10]);
+ double p1213 = fputil::multiply_add(d, POLY_M3N[13], POLY_M3N[12]);
+ double p1415 = fputil::multiply_add(d, POLY_M3N[15], POLY_M3N[14]);
+ double p1617 = fputil::multiply_add(d, POLY_M3N[17], POLY_M3N[16]);
+ double q03 = fputil::multiply_add(d2, p23, p01);
+ double q47 = fputil::multiply_add(d2, p67, p45);
+ double q811 = fputil::multiply_add(d2, p1011, p89);
+ double q1215 = fputil::multiply_add(d2, p1415, p1213);
+ double r07 = fputil::multiply_add(d4, q47, q03);
+ double r815 = fputil::multiply_add(d4, q1215, q811);
+ double s015 = fputil::multiply_add(d8, r815, r07);
+ double poly = fputil::multiply_add(d16, p1617, s015);
+ double pa = xd * (xd + 1.0);
+ double pb = (xd + 2.0) * (xd + 3.0);
+ double prod = (pa * pb) * (xd + 4.0);
+ double aprod = prod < 0.0 ? -prod : prod;
+ lgamma_val = fputil::multiply_add(xd + 3.0, poly, -lg_ln(aprod));
}
}
} else {
@@ -266,10 +379,6 @@ LIBC_INLINE float lgammaf(float x) {
// 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)) {
@@ -282,80 +391,104 @@ LIBC_INLINE float lgammaf(float x) {
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);
-
+ // No cancellation in (x-0.5)*(log(x)-1): log(x)-1 > 0.2 on this range.
+ // Relative error ~2^-48; anything closer to a rounding boundary is in
+ // the exceptional cases tables below.
+ double lz = lg_ln(abs_xd);
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);
+ lgamma_val = fputil::multiply_add(xm, lz - 1.0, 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;
+ // For |x| >= 2^20 the 1/(12x) correction is below ~2^-47 of the result;
+ // skip it and the 1/x divide.
+ if (x_abs < 0x49800000u) {
+ double inv_x = 1.0 / abs_xd;
+ double inv_x2 = inv_x * inv_x;
+ if (x_abs > 0x44fa0000u) {
+ constexpr fputil::ExceptValues<float, 3> LGAMMAF_EXCEPTS_BERN2{{
+ // input, toward-zero result, RU, RD, RN
+ {0x46541516u, 0x47e1c01bu, 1, 0, 0},
+ {0x46b16323u, 0x48483adeu, 1, 0, 1},
+ {0xc6f7e151u, 0xc89116deu, 0, 1, 0},
+ }};
+ if (auto r = LGAMMAF_EXCEPTS_BERN2.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ // |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) {
+ // Exceptional cases of this range.
+ constexpr fputil::ExceptValues<float, 2> LGAMMAF_EXCEPTS_BERN4{{
+ // input, toward-zero result, RU, RD, RN
+ {0x449acf07u, 0x45ecd680u, 1, 0, 1},
+ {0xc33139a3u, 0xc43991afu, 0, 1, 0},
+ }};
+ if (auto r = LGAMMAF_EXCEPTS_BERN4.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ // |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 {
+ // Exceptional cases of this range.
+ constexpr fputil::ExceptValues<float, 2> LGAMMAF_EXCEPTS_B10{{
+ // input, toward-zero result, RU, RD, RN
+ {0x42468b59u, 0x430f25a7u, 1, 0, 0},
+ {0xc134eb14u, 0xc1875615u, 0, 1, 0},
+ }};
+ if (auto r = LGAMMAF_EXCEPTS_B10.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+ // |x| in (3.373, 73]: degree-10 monomial fit of
+ // h(s) = stir_resid(1/sqrt(s)) / sqrt(s), s = 1/x^2, on [0, 0.088];
+ // correction = h(s) * (1/x). Max error 2^-53.2.
+ constexpr double MID_B10 = 0x1.6880000000000p-5;
+ constexpr double POLY_B10[11] = {
+ 0x1.54d6b78cee955p-4, -0x1.635a5fb0cdf9fp-9,
+ 0x1.7b5253f44b255p-11, -0x1.f32907e7a7adap-12,
+ 0x1.1f269b6438739p-11, -0x1.e95dc64c5042cp-11,
+ 0x1.17ebce09c49a7p-9, -0x1.91599e7728747p-8,
+ 0x1.5a01e1f07b127p-6, -0x1.90e956754bc6ap-4,
+ 0x1.dfbed80c6f035p-2};
+ double u = inv_x2 - MID_B10;
+ double u2 = u * u, u4 = u2 * u2, u8 = u4 * u4;
+ double p01 = fputil::multiply_add(u, POLY_B10[1], POLY_B10[0]);
+ double p23 = fputil::multiply_add(u, POLY_B10[3], POLY_B10[2]);
+ double p45 = fputil::multiply_add(u, POLY_B10[5], POLY_B10[4]);
+ double p67 = fputil::multiply_add(u, POLY_B10[7], POLY_B10[6]);
+ double p89 = fputil::multiply_add(u, POLY_B10[9], POLY_B10[8]);
+ double q03 = fputil::multiply_add(u2, p23, p01);
+ double q47 = fputil::multiply_add(u2, p67, p45);
+ double q810 = fputil::multiply_add(u2, POLY_B10[10], p89);
+ double r07 = fputil::multiply_add(u4, q47, q03);
+ double poly = fputil::multiply_add(u8, q810, r07);
+ lgamma_val += inv_x * poly;
}
- double poly = fputil::multiply_add(u, b1, CHEB_B10[0] - b2);
- lgamma_val += inv_x * poly;
+ } else {
+ constexpr fputil::ExceptValues<float, 3> LGAMMAF_EXCEPTS_HUGE{{
+ // input, toward-zero result, RU, RD, RN
+ {0x65fca09fu, 0x68cead59u, 1, 0, 1},
+ {0x716e5dd5u, 0x747e2bb9u, 1, 0, 0},
+ {0x77ac5674u, 0x7acf27b2u, 1, 0, 1},
+ }};
+ if (auto r = LGAMMAF_EXCEPTS_HUGE.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
}
if (xbits.is_neg()) {
- // Reflection: lgamma(x) = log(pi) - lgamma(|x|) -
- // log(|x|*|sin(pi*frac_x)|)
+ // Reflection: lgamma(x) = log(pi) - lgamma(|x|) - log(|x|) -
+ // log(|sin(pi*frac_x)|). Reusing the already-computed lz = log(|x|)
+ // keeps the multiply out of the serial sinpi -> log chain.
double frac_x = xd - fputil::floor(xd);
- lgamma_val = 0x1.250d048e7a1bdp+0 - lgamma_val;
- lgamma_val -= lg_ln(lg_sinpi(frac_x) * abs_xd);
+ lgamma_val = (0x1.250d048e7a1bdp+0 - lgamma_val) - lz;
+ lgamma_val -= lg_ln(lg_sinpi(frac_x));
}
}
>From 7ce4e93f479322ec3c0b16aabea263900d1eb488 Mon Sep 17 00:00:00 2001
From: Anonmiraj <ezzibrahimx at gmail.com>
Date: Tue, 28 Jul 2026 17:02:20 +0300
Subject: [PATCH 3/3] fix formatting
---
libc/src/__support/math/lgammaf.h | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/libc/src/__support/math/lgammaf.h b/libc/src/__support/math/lgammaf.h
index e2f0a16aa4b2b..e5c0f818b44bd 100644
--- a/libc/src/__support/math/lgammaf.h
+++ b/libc/src/__support/math/lgammaf.h
@@ -405,7 +405,7 @@ LIBC_INLINE float lgammaf(float x) {
double inv_x2 = inv_x * inv_x;
if (x_abs > 0x44fa0000u) {
constexpr fputil::ExceptValues<float, 3> LGAMMAF_EXCEPTS_BERN2{{
- // input, toward-zero result, RU, RD, RN
+ // input, toward-zero result, RU, RD, RN
{0x46541516u, 0x47e1c01bu, 1, 0, 0},
{0x46b16323u, 0x48483adeu, 1, 0, 1},
{0xc6f7e151u, 0xc89116deu, 0, 1, 0},
@@ -420,7 +420,7 @@ LIBC_INLINE float lgammaf(float x) {
} else if (x_abs > 0x42920000u) {
// Exceptional cases of this range.
constexpr fputil::ExceptValues<float, 2> LGAMMAF_EXCEPTS_BERN4{{
- // input, toward-zero result, RU, RD, RN
+ // input, toward-zero result, RU, RD, RN
{0x449acf07u, 0x45ecd680u, 1, 0, 1},
{0xc33139a3u, 0xc43991afu, 0, 1, 0},
}};
@@ -438,7 +438,7 @@ LIBC_INLINE float lgammaf(float x) {
} else {
// Exceptional cases of this range.
constexpr fputil::ExceptValues<float, 2> LGAMMAF_EXCEPTS_B10{{
- // input, toward-zero result, RU, RD, RN
+ // input, toward-zero result, RU, RD, RN
{0x42468b59u, 0x430f25a7u, 1, 0, 0},
{0xc134eb14u, 0xc1875615u, 0, 1, 0},
}};
@@ -472,7 +472,7 @@ LIBC_INLINE float lgammaf(float x) {
}
} else {
constexpr fputil::ExceptValues<float, 3> LGAMMAF_EXCEPTS_HUGE{{
- // input, toward-zero result, RU, RD, RN
+ // input, toward-zero result, RU, RD, RN
{0x65fca09fu, 0x68cead59u, 1, 0, 1},
{0x716e5dd5u, 0x747e2bb9u, 1, 0, 0},
{0x77ac5674u, 0x7acf27b2u, 1, 0, 1},
More information about the libc-commits
mailing list