[libc-commits] [libc] [llvm] [libc][math] Implement half precision lgamma function (PR #192834)
via libc-commits
libc-commits at lists.llvm.org
Sat Jul 18 14:59:47 PDT 2026
================
@@ -0,0 +1,307 @@
+//===-- 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] = {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 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) {
----------------
lntue wrote:
nit: add comment on its real value.
https://github.com/llvm/llvm-project/pull/192834
More information about the libc-commits
mailing list