Fix catastrophic cancellation error in poisson logpdf when k and mu are large

include/boost/math/distributions/poisson.hpp

@@ -50,6 +50,7 @@
 #include <boost/math/special_functions/factorials.hpp> // factorials.
 #include <boost/math/tools/roots.hpp> // for root finding.
 #include <boost/math/distributions/detail/inv_discrete_quantile.hpp>
+#include <boost/math/constants/constants.hpp>
 
 namespace boost
 {
@@ -141,6 +142,55 @@ namespace boost
         return true;
       } // bool check_dist_and_prob
 
+      template <class RealType>
+      BOOST_MATH_GPU_ENABLED inline RealType stirlerr(const RealType& n) {
+        BOOST_MATH_STD_USING // for ADL of std functions.
+        using boost::math::lgamma;
+
+        // Stirling's series coefficients
+        const RealType S0 = RealType(1)/12;
+        const RealType S1 = RealType(1)/360;
+        const RealType S2 = RealType(1)/1260;
+        const RealType S3 = RealType(1)/1680;
+        const RealType S4 = RealType(1)/1188;
+
+        // Use Stirling's series if n is small; use the direct formula otherwise
+        bool is_small = n < 15;
+        if (is_small) {
+          return lgamma(n + 1) - (n * log(n) - n + 0.5 * log(2 * boost::math::constants::pi<RealType>() * n));
+        } else {
+          RealType n2 = n * n;
+          return (S0 - (S1 - (S2 - (S3 - S4/n2)/n2)/n2)/n2)/n;
+        }
+
+      }
+
+      template <class RealType>
+      BOOST_MATH_GPU_ENABLED inline RealType bd0(const RealType& mean, const RealType& k) {
+        BOOST_MATH_STD_USING // for ADL of std functions.
+
+        // Calculate v = (k - mean) / (k + mean) from Loader (2000) approximation
+        bool is_close = abs(k - mean) < RealType(0.1) * (k + mean);
+
+        if (is_close) { // Use the series approximation if |v| < 0.1
+          RealType v = (k - mean) / (k + mean);
+          RealType v2 = v * v;
+          RealType series_term = ((k - mean) * (k - mean)) / (k + mean);
+
+          RealType term = 2 * k * v;
+          for (int i = 1; i < 11; ++i) {
+            term *= v2;
+            series_term += term / (2 * i + 1);
+          }
+          return series_term;
+        } else {
+          // Use the direct formula if |v| >= 0.1
+          RealType direct = (k == 0) ? RealType(0) : k * log(k / mean) + mean - k;
+          return direct;
+        }
+
+      }
+
     } // namespace poisson_detail
 
     BOOST_MATH_EXPORT template <class RealType = double, class Policy = policies::policy<> >
@@ -304,7 +354,8 @@ namespace boost
       // Special case where k and lambda are both positive
       if(k > 0 && mean > 0)
       {
-        return -lgamma(k+1) + k*log(mean) - mean;
+        // Use the Loader (2000) saddle-point approximation for logpdf calculation
+        return -poisson_detail::stirlerr(k) - poisson_detail::bd0(mean, k) - RealType(0.5) * log(2 * boost::math::constants::pi<RealType>() * k);
       }
 
       result = log(pdf(dist, k));

test/test_poisson.cpp

@@ -244,7 +244,68 @@ void test_spots(RealType)
       static_cast<RealType>(20)),   //  K>> mean 
       log(static_cast<RealType>(8.277463646553730E-009)), // probability.
          tolerance);
-
+
+   // New test cases for Loader (2000) saddle-point approximation. Probs already
+   // in log space. Values calculated using mpmath (1000-digit precision).
+   BOOST_CHECK_CLOSE(
+     logpdf(poisson_distribution<RealType>(static_cast<RealType>(14)), // mean 14.
+      static_cast<RealType>(14)),
+      static_cast<RealType>(-2.244418568125061), // probability (already in log space).
+         tolerance);
+
+   BOOST_CHECK_CLOSE(
+      logpdf(poisson_distribution<RealType>(static_cast<RealType>(20)), // mean 20.
+         static_cast<RealType>(18)),
+         static_cast<RealType>(-2.472264284061216), // probability (already in log space).
+         tolerance);
+
+   // Cases below require around 15+ significant decimal digits to represent
+  // k / mean meaningfully, so skip for float.
+  if (std::numeric_limits<RealType>::digits10 > 15)
+  {
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(1000000)), // mean 1000000.
+        static_cast<RealType>(1300000)),
+        static_cast<RealType>(-41081.501683746894),
+           tolerance);
+
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(1e8)), // mean 1e8.
+        static_cast<RealType>(8e7)),
+        static_cast<RealType>(-2148525.91257035),
+           tolerance);
+
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(1e10)), // mean 1e10.
+        static_cast<RealType>(105e9)),
+        static_cast<RealType>(-151894402015.7727),
+           tolerance);
+
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(1e15)), // mean 1e15.
+        static_cast<RealType>(8e14)), // |v| > 0.1 boundary
+        static_cast<RealType>(-21485158948650.273),
+           tolerance);
+
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(1e15)), // mean 1e15.
+        static_cast<RealType>(12e14)), // |v| < 0.1 boundary
+        static_cast<RealType>(-18785868152763.832),
+           tolerance);
+
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(1e16)), // mean 1e16.
+        static_cast<RealType>(1e16)), // old formula returns 0.0 here
+        static_cast<RealType>(-19.339619277157038),
+           tolerance);
+
+    BOOST_CHECK_CLOSE(
+       logpdf(poisson_distribution<RealType>(static_cast<RealType>(5e15)), // mean 5e15.
+        static_cast<RealType>(5e15)),
+        static_cast<RealType>(-18.993045686877064),
+           tolerance);
+  }
+
   // CDF
   BOOST_CHECK_CLOSE(
      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.