From 4dcddfb8017818fa393a152cfe9d8aa53296e6fa Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 18 Apr 2022 09:22:31 -0700
Subject: [PATCH] Replace uses of log(tgamma(x)) with lgamma(x)

---
 include/boost/math/distributions/gamma.hpp         | 14 +++-----------
 include/boost/math/distributions/inverse_gamma.hpp |  4 ++--
 include/boost/math/distributions/poisson.hpp       |  4 ++--
 3 files changed, 7 insertions(+), 15 deletions(-)

diff --git a/include/boost/math/distributions/gamma.hpp b/include/boost/math/distributions/gamma.hpp
index 83c286595..28b7c55b0 100644
--- a/include/boost/math/distributions/gamma.hpp
+++ b/include/boost/math/distributions/gamma.hpp
@@ -152,7 +152,7 @@ template <class RealType, class Policy>
 inline RealType logpdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
 {
    BOOST_MATH_STD_USING  // for ADL of std functions
-   using boost::math::tgamma;
+   using boost::math::lgamma;
 
    static const char* function = "boost::math::logpdf(const gamma_distribution<%1%>&, %1%)";
 
@@ -169,16 +169,8 @@ inline RealType logpdf(const gamma_distribution<RealType, Policy>& dist, const R
    {
       return std::numeric_limits<RealType>::quiet_NaN();
    }
-   
-   // The following calculation does not always work with float so take the naive road out
-   BOOST_IF_CONSTEXPR(std::is_same<RealType, float>::value)
-   {
-      result = log(pdf(dist, x));
-   }
-   else
-   {
-      result = -k*log(theta) + (k-1)*log(x) - log(tgamma(k)) - (x/theta);
-   }
+
+   result = -k*log(theta) + (k-1)*log(x) - lgamma(k) - (x/theta);
    
    return result;
 } // logpdf
diff --git a/include/boost/math/distributions/inverse_gamma.hpp b/include/boost/math/distributions/inverse_gamma.hpp
index eb348180d..8c9e4763d 100644
--- a/include/boost/math/distributions/inverse_gamma.hpp
+++ b/include/boost/math/distributions/inverse_gamma.hpp
@@ -198,7 +198,7 @@ template <class RealType, class Policy>
 inline RealType logpdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
 {
    BOOST_MATH_STD_USING  // for ADL of std functions
-   using boost::math::tgamma;
+   using boost::math::lgamma;
 
    static const char* function = "boost::math::logpdf(const inverse_gamma_distribution<%1%>&, %1%)";
 
@@ -228,7 +228,7 @@ inline RealType logpdf(const inverse_gamma_distribution<RealType, Policy>& dist,
       return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
    }
 
-   return shape * log(scale) + (-shape-1)*log(x) - log(tgamma(shape)) - (scale/x);
+   return shape * log(scale) + (-shape-1)*log(x) - lgamma(shape) - (scale/x);
 } // pdf
 
 template <class RealType, class Policy>
diff --git a/include/boost/math/distributions/poisson.hpp b/include/boost/math/distributions/poisson.hpp
index dd5e8062f..5507360e8 100644
--- a/include/boost/math/distributions/poisson.hpp
+++ b/include/boost/math/distributions/poisson.hpp
@@ -279,7 +279,7 @@ namespace boost
       BOOST_FPU_EXCEPTION_GUARD
 
       BOOST_MATH_STD_USING // for ADL of std functions.
-      using boost::math::tgamma;
+      using boost::math::lgamma;
 
       RealType mean = dist.mean();
       // Error check:
@@ -302,7 +302,7 @@ namespace boost
       // Special case where k and lambda are both positive
       if(k > 0 && mean > 0)
       {
-        return -log(tgamma(k+1)) + k*log(mean) - mean;
+        return -lgamma(k+1) + k*log(mean) - mean;
       }
 
       result = log(pdf(dist, k));