Add ligamma and ligamma_lower, the logarithms of tgamma(a, x) and tgamma_lower(a, x)

doc/sf/igamma.qbk

@@ -43,6 +43,18 @@
 
    template <class T1, class T2, class ``__Policy``>
    BOOST_MATH_GPU_ENABLED ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&);
+
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma(T1 a, T2 z);
+
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma(T1 a, T2 z, const ``__Policy``&);
+
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma_lower(T1 a, T2 z);
+
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma_lower(T1 a, T2 z, const ``__Policy``&);
 
    }} // namespaces
 
@@ -130,6 +142,22 @@ Returns the full (non-normalised) upper incomplete gamma function of a and z:
 
 [equation igamma1]
 
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma(T1 a, T2 z);
+
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma(T1 a, T2 z, const ``__Policy``&);
+
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma_lower(T1 a, T2 z);
+
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` ligamma_lower(T1 a, T2 z, const ``__Policy``&);
+
+Return the natural logarithms of `tgamma(a, z)` and `tgamma_lower(a, z)`. They stay finite and accurate where those
+overflow or underflow: for example `ligamma(1000, 1000)` is about 5905 in double precision, where `log(tgamma(1000, 1000))`
+raises an overflow error. Where the result is exactly zero, its logarithm, -infinity, is reported as an __overflow_error.
+
 [h4 Accuracy]
 
 The following tables give peak and mean relative errors in over various domains of
@@ -303,8 +331,16 @@ in which case an expansion, for large x, of the (non-normalised) upper
 incomplete gamma function is used. The return is then normalised by subtracting
 the log of the gamma function and adding /a log(x)-x-log(x)/.  
 
+`ligamma` and `ligamma_lower` take the logarithm of `tgamma(a, z)` or `tgamma_lower(a, z)` where that is representable.
+Otherwise they add `lgamma(a)` to the logarithm of the regularized function, or of `log1p` of minus its complement when that
+is the smaller, except where the regularized function underflows: then they use /a log(z) - z/ plus the logarithm of the
+continued fraction (upper) or series (lower) that it multiplies, the same normalization by /z[super a]e[super -z]/ as the
+G-function of Abergel and Moisan.
+
 [h4 References]
 
+* R. Abergel and L. Moisan, Algorithm 1006: Fast and Accurate Evaluation of a Generalized Incomplete Gamma Function,
+ACM TOMS, Vol 46, No 1, 2020.
 * N. M. Temme, A Set of Algorithms for the Incomplete Gamma Functions,
 Probability in the Engineering and Informational Sciences, 8, 1994.
 * N. M. Temme, The Asymptotic Expansion of the Incomplete Gamma Functions,

include/boost/math/special_functions/gamma.hpp

@@ -237,7 +237,7 @@ BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T gamma_imp(T z, const Policy& pol
             return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
          /*
          * Result can never be subnormal as we have a value > 1 in the numerator:
-         if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
+         if((BOOST_MATH_NAMESPACE::fpclassify)(result) == (int)FP_SUBNORMAL)
             return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
             */
          BOOST_MATH_INSTRUMENT_VARIABLE(result);
@@ -1862,6 +1862,67 @@ BOOST_MATH_GPU_ENABLED T lgamma_incomplete_lower_imp(T a, T x, const Policy& pol
    return log1p(-gamma_q(a, x, pol), pol);
 }
 
+//
+// Logarithm of the non-normalised upper (upper == true) or lower incomplete gamma function.
+//
+// Where the non-normalised value is representable it is the logarithm of that. Otherwise this
+// uses the regularised value, the target T = Q (upper) or P (lower), and its complement C:
+//
+//    T >= 1/2:                 lgamma(a) + log1p(-C)
+//    T representable:          lgamma(a) + log(T)
+//    T underflows:             a log(x) - x + log(continued fraction), upper, x > a
+//                              a log(x) - x + log(series / a),         lower, x < a
+//
+// The last two are the methods that gamma_incomplete_imp uses in logs for large a, here applied
+// only far out in the tail, where they converge quickly. As in Abergel and Moisan's G-function
+// (ACM TOMS 46(1), Algorithm 1006), the prefix x^a e^-x is never formed.
+//
+template <class T, class Policy>
+BOOST_MATH_GPU_ENABLED T ligamma_imp(T a, T x, const Policy& pol, bool upper)
+{
+   BOOST_MATH_STD_USING
+   const char* function = upper ? "boost::math::ligamma<%1%>(%1%, %1%)" : "boost::math::ligamma_lower<%1%>(%1%, %1%)";
+
+   // These checks also reject NaN arguments:
+   if(!(BOOST_MATH_NAMESPACE::isfinite)(a))
+      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be finite (got a=%1%).", a, pol);
+   if(a <= 0)
+      return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
+   if(!(x >= 0))
+      return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
+
+   // At the ends of the range, the target is either the complete gamma function or zero:
+   bool target_is_zero = upper ? (BOOST_MATH_NAMESPACE::isinf)(x) : (x == 0);
+   if(target_is_zero)
+      return -policies::raise_overflow_error<T>(function, nullptr, pol);
+   if(upper ? (x == 0) : (BOOST_MATH_NAMESPACE::isinf)(x))
+      return BOOST_MATH_NAMESPACE::lgamma(a, pol);
+
+   // Under- and overflow are ours to handle, so stop the evaluations below from reporting them:
+   typedef typename policies::normalise<Policy, policies::underflow_error<policies::ignore_error>, policies::overflow_error<policies::ignore_error> >::type quiet_policy;
+
+   if(a < max_factorial<T>::value)
+   {
+      // The non-normalised value cannot overflow, as it is at most tgamma(a):
+      T value = gamma_incomplete_imp(a, x, false, upper, quiet_policy(), static_cast<T*>(nullptr));
+      if(value >= tools::min_value<T>())
+         return log(value);
+   }
+   T target = gamma_incomplete_imp(a, x, true, upper, quiet_policy(), static_cast<T*>(nullptr));
+   if(target > 0.5f)
+      return BOOST_MATH_NAMESPACE::lgamma(a, pol) + BOOST_MATH_NAMESPACE::log1p(-gamma_incomplete_imp(a, x, true, !upper, quiet_policy(), static_cast<T*>(nullptr)), pol);
+   if(target >= tools::min_value<T>())
+      return BOOST_MATH_NAMESPACE::lgamma(a, pol) + log(target);
+
+   // The target underflows, so x is far out in its tail:
+   T result = a * log(x) - x;
+   if(upper)
+      result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
+   else
+      result += log(detail::lower_gamma_series(a, x, pol, T(0)) / a);
+   return result;
+}
+
 //
 // Ratios of two gamma functions:
 //
@@ -2526,6 +2587,52 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> lgamma_p(T1 a, T2 z)
 {
    return lgamma_p(a, z, policies::policy<>());
 }
+
+BOOST_MATH_EXPORT template <class T1, class T2, class Policy>
+BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> ligamma(T1 a, T2 z, const Policy&)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef tools::promote_args_t<T1, T2> result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+      detail::ligamma_imp(static_cast<value_type>(a), static_cast<value_type>(z), forwarding_policy(), true), "boost::math::ligamma<%1%>(%1%, %1%)");
+}
+
+BOOST_MATH_EXPORT template <class T1, class T2>
+BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> ligamma(T1 a, T2 z)
+{
+   return ligamma(a, z, policies::policy<>());
+}
+
+BOOST_MATH_EXPORT template <class T1, class T2, class Policy>
+BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> ligamma_lower(T1 a, T2 z, const Policy&)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef tools::promote_args_t<T1, T2> result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+      detail::ligamma_imp(static_cast<value_type>(a), static_cast<value_type>(z), forwarding_policy(), false), "boost::math::ligamma_lower<%1%>(%1%, %1%)");
+}
+
+BOOST_MATH_EXPORT template <class T1, class T2>
+BOOST_MATH_GPU_ENABLED inline tools::promote_args_t<T1, T2> ligamma_lower(T1 a, T2 z)
+{
+   return ligamma_lower(a, z, policies::policy<>());
+}
 //
 // Regularised lower incomplete gamma:
 //

include/boost/math/special_functions/math_fwd.hpp

@@ -594,6 +594,18 @@ BOOST_MATH_NAMESPACE_BEGIN // Math functions (in roughly alphabetic order).
    BOOST_MATH_EXPORT template <class RT1, class RT2, class Policy>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> lgamma_p(RT1 a, RT2 z, const Policy&);
 
+   BOOST_MATH_EXPORT template <class RT1, class RT2>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> ligamma(RT1 a, RT2 z);
+
+   BOOST_MATH_EXPORT template <class RT1, class RT2, class Policy>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> ligamma(RT1 a, RT2 z, const Policy&);
+
+   BOOST_MATH_EXPORT template <class RT1, class RT2>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> ligamma_lower(RT1 a, RT2 z);
+
+   BOOST_MATH_EXPORT template <class RT1, class RT2, class Policy>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> ligamma_lower(RT1 a, RT2 z, const Policy&);
+
    BOOST_MATH_EXPORT template <class RT1, class RT2>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2> gamma_p(RT1 a, RT2 z);
 
@@ -1566,6 +1578,12 @@ BOOST_MATH_NAMESPACE_BEGIN // Math functions (in roughly alphabetic order).
 \
    template <class RT1, class RT2>\
    BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2> lgamma_p(RT1 a, RT2 z){ return BOOST_MATH_NAMESPACE::lgamma_p(a, z, Policy()); }\
+\
+   template <class RT1, class RT2>\
+   BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2> ligamma(RT1 a, RT2 z){ return BOOST_MATH_NAMESPACE::ligamma(a, z, Policy()); }\
+\
+   template <class RT1, class RT2>\
+   BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2> ligamma_lower(RT1 a, RT2 z){ return BOOST_MATH_NAMESPACE::ligamma_lower(a, z, Policy()); }\
 \
    template <class RT1, class RT2>\
    BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2> gamma_p(RT1 a, RT2 z){ return BOOST_MATH_NAMESPACE::gamma_p(a, z, Policy()); }\

reporting/performance/ligamma_performance.cpp

@@ -0,0 +1,97 @@
+//  (C) Copyright Nick Thompson 2026.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+//
+// Compares ligamma(a, x) with log(tgamma(a, x)) where both are finite, even in float: a log-uniform in [0.5, 30],
+// x/a log-uniform in [0.1, 10]. Then, for a log-uniform in [1e3, 1e5], where tgamma(a, x) overflows,
+// compares ligamma with lgamma(a) + lgamma_q(a, x).
+//
+// Promotion is disabled so that each type is timed in its own arithmetic: by default,
+// float is evaluated in double, and double in long double where that is wider.
+//
+
+#include <cmath>
+#include <random>
+#include <utility>
+#include <vector>
+#include <benchmark/benchmark.h>
+#include <boost/math/special_functions/gamma.hpp>
+
+using no_promotion = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
+
+template <typename Real>
+std::vector<std::pair<Real, Real>> arguments(double log_a_min, double log_a_max)
+{
+    std::mt19937_64 mt {42};
+    std::uniform_real_distribution<double> log_a(log_a_min, log_a_max);
+    std::uniform_real_distribution<double> log_ratio(-1, 1);
+    std::vector<std::pair<Real, Real>> v(1024);
+    for (auto& p : v)
+    {
+        double a = std::pow(10.0, log_a(mt));
+        p = {static_cast<Real>(a), static_cast<Real>(a * std::pow(10.0, log_ratio(mt)))};
+    }
+    return v;
+}
+
+template <typename Real>
+void ligamma_performance(benchmark::State& state)
+{
+    auto v = arguments<Real>(std::log10(0.5), std::log10(30.0));
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(boost::math::ligamma(p.first, p.second, no_promotion()));
+    }
+}
+
+template <typename Real>
+void log_tgamma_performance(benchmark::State& state)
+{
+    using std::log;
+    auto v = arguments<Real>(std::log10(0.5), std::log10(30.0));
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(log(boost::math::tgamma(p.first, p.second, no_promotion())));
+    }
+}
+
+template <typename Real>
+void ligamma_large_a_performance(benchmark::State& state)
+{
+    auto v = arguments<Real>(3, 5);
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(boost::math::ligamma(p.first, p.second, no_promotion()));
+    }
+}
+
+template <typename Real>
+void lgamma_plus_lgamma_q_large_a_performance(benchmark::State& state)
+{
+    auto v = arguments<Real>(3, 5);
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(boost::math::lgamma(p.first, no_promotion()) + boost::math::lgamma_q(p.first, p.second, no_promotion()));
+    }
+}
+
+BENCHMARK_TEMPLATE(ligamma_performance, float);
+BENCHMARK_TEMPLATE(ligamma_performance, double);
+BENCHMARK_TEMPLATE(ligamma_performance, long double);
+BENCHMARK_TEMPLATE(log_tgamma_performance, float);
+BENCHMARK_TEMPLATE(log_tgamma_performance, double);
+BENCHMARK_TEMPLATE(log_tgamma_performance, long double);
+BENCHMARK_TEMPLATE(ligamma_large_a_performance, double);
+BENCHMARK_TEMPLATE(lgamma_plus_lgamma_q_large_a_performance, double);
+
+BENCHMARK_MAIN();

test/Jamfile.v2

@@ -489,6 +489,7 @@ test-suite special_fun :
           <toolset>intel:<pch>off
         : test_ibeta_inv_ab_real_concept3  ]
    [ run test_igamma.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
+   [ run test_ligamma.cpp ]
    [ run test_igamma_inv.cpp test_instances//test_instances pch_light  /boost/test//boost_unit_test_framework
         : # command line
         : # input files

test/compile_test/instantiate.hpp

@@ -265,6 +265,8 @@ void instantiate(RealType)
    boost::math::gamma_q(v1, v2);
    boost::math::lgamma_q(v1, v2);
    boost::math::lgamma_p(v1, v2);
+   boost::math::ligamma(v1, v2);
+   boost::math::ligamma_lower(v1, v2);
    boost::math::gamma_p_inv(v1, v2);
    boost::math::gamma_q_inv(v1, v2);
    boost::math::gamma_p_inva(v1, v2);
@@ -549,6 +551,8 @@ void instantiate(RealType)
    boost::math::gamma_q(v1 * 1, v2 + 0);
    boost::math::lgamma_q(v1 * 1, v2 + 0);
    boost::math::lgamma_p(v1 * 1, v2 + 0);
+   boost::math::ligamma(v1 * 1, v2 + 0);
+   boost::math::ligamma_lower(v1 * 1, v2 + 0);
    boost::math::gamma_p_inv(v1 * 1, v2 + 0);
    boost::math::gamma_q_inv(v1 * 1, v2 + 0);
    boost::math::gamma_p_inva(v1 * 1, v2 + 0);
@@ -805,6 +809,8 @@ void instantiate(RealType)
    boost::math::gamma_q(v1, v2, pol);
    boost::math::lgamma_q(v1, v2, pol);
    boost::math::lgamma_p(v1, v2, pol);
+   boost::math::ligamma(v1, v2, pol);
+   boost::math::ligamma_lower(v1, v2, pol);
    boost::math::gamma_p_inv(v1, v2, pol);
    boost::math::gamma_q_inv(v1, v2, pol);
    boost::math::gamma_p_inva(v1, v2, pol);
@@ -1087,6 +1093,8 @@ void instantiate(RealType)
    test::gamma_q(v1, v2);
    test::lgamma_q(v1, v2);
    test::lgamma_p(v1, v2);
+   test::ligamma(v1, v2);
+   test::ligamma_lower(v1, v2);
    test::gamma_p_inv(v1, v2);
    test::gamma_q_inv(v1, v2);
    test::gamma_p_inva(v1, v2);
@@ -1373,6 +1381,8 @@ void instantiate_mixed(RealType)
    boost::math::gamma_q(i, s);
    boost::math::lgamma_q(i, s);
    boost::math::lgamma_p(i, s);
+   boost::math::ligamma(i, s);
+   boost::math::ligamma_lower(i, s);
    boost::math::gamma_q(fr, lr);
    boost::math::gamma_p_inv(i, fr);
    boost::math::gamma_q_inv(s, fr);
@@ -1593,6 +1603,8 @@ void instantiate_mixed(RealType)
    boost::math::gamma_q(i, s, pol);
    boost::math::lgamma_q(i, s, pol);
    boost::math::lgamma_p(i, s, pol);
+   boost::math::ligamma(i, s, pol);
+   boost::math::ligamma_lower(i, s, pol);
    boost::math::gamma_q(fr, lr, pol);
    boost::math::gamma_p_inv(i, fr, pol);
    boost::math::gamma_q_inv(s, fr, pol);
@@ -1812,6 +1824,8 @@ void instantiate_mixed(RealType)
    test::gamma_q(fr, lr);
    test::lgamma_q(fr, lr);
    test::lgamma_p(fr, lr);
+   test::ligamma(fr, lr);
+   test::ligamma_lower(fr, lr);
    test::gamma_p_inv(i, fr);
    test::gamma_q_inv(s, fr);
    test::gamma_p_inva(i, lr);