remove internal function in lccdf and add two tests for expanded range

spinkney · spinkney · commit df0960030bbc · 2025-12-26T08:55:09.000-05:00
diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp
@@ -29,6 +29,66 @@ struct log_gamma_q_result {
   T dlog_q_da;  ///< d/da log(Q(a,z))
 };
 
+namespace internal {
+
+/**
+ * Compute log(Q(a,z)) using continued fraction expansion for upper incomplete
+ * gamma function.
+ *
+ * @tparam T_a Type of shape parameter a (double or fvar types)
+ * @tparam T_z Type of value parameter z (double or fvar types)
+ * @param a Shape parameter
+ * @param z Value at which to evaluate
+ * @param max_steps Maximum number of continued fraction iterations
+ * @param precision Convergence threshold
+ * @return log(Q(a,z)) with same type as T_a and T_z
+ */
+template <typename T_a, typename T_z>
+inline auto log_q_gamma_cf(const T_a& a, const T_z& z, int max_steps = 250,
+                           double precision = 1e-16) {
+  using stan::math::lgamma;
+  using stan::math::log;
+  using stan::math::value_of;
+  using std::fabs;
+  using T_return = return_type_t<T_a, T_z>;
+
+  const T_return a_ret = a;
+  const T_return z_ret = z;
+  const auto log_prefactor = a_ret * log(z_ret) - z_ret - lgamma(a_ret);
+
+  auto b = z_ret + 1.0 - a_ret;
+  auto C = (fabs(value_of(b)) >= EPSILON) ? b : T_return(EPSILON);
+  auto D = T_return(0.0);
+  auto f = C;
+
+  for (int i = 1; i <= max_steps; ++i) {
+    auto an = -i * (i - a_ret);
+    b += 2.0;
+
+    D = b + an * D;
+    if (fabs(value_of(D)) < EPSILON) {
+      D = T_return(EPSILON);
+    }
+    C = b + an / C;
+    if (fabs(value_of(C)) < EPSILON) {
+      C = T_return(EPSILON);
+    }
+
+    D = 1.0 / D;
+    auto delta = C * D;
+    f *= delta;
+
+    const double delta_m1 = value_of(fabs(value_of(delta) - 1.0));
+    if (delta_m1 < precision) {
+      break;
+    }
+  }
+
+  return log_prefactor - log(f);
+}
+
+}  // namespace internal
+
 /**
  * Compute log(Q(a,z)) and its gradient with respect to a using continued
  * fraction expansion, where Q(a,z) = Gamma(a,z) / Gamma(a) is the regularized
@@ -50,7 +110,6 @@ template <typename T_a, typename T_z>
 inline log_gamma_q_result<return_type_t<T_a, T_z>> log_gamma_q_dgamma(
     const T_a& a, const T_z& z, int max_steps = 250, double precision = 1e-16) {
   using std::exp;
-  using std::fabs;
   using std::log;
   using T_return = return_type_t<T_a, T_z>;
 
@@ -61,39 +120,8 @@ inline log_gamma_q_result<return_type_t<T_a, T_z>> log_gamma_q_dgamma(
 
   // For z > a + 1, use continued fraction for better numerical stability
   if (z_dbl > a_dbl + 1.0) {
-    // Continued fraction for Q(a,z) in log space
-    // log(Q(a,z)) = log_prefactor - log(continued_fraction)
-    const double log_prefactor = a_dbl * log(z_dbl) - z_dbl - lgamma(a_dbl);
-
-    double b = z_dbl + 1.0 - a_dbl;
-    double C = (fabs(b) >= EPSILON) ? b : EPSILON;
-    double D = 0.0;
-    double f = C;
-
-    for (int i = 1; i <= max_steps; ++i) {
-      const double an = -i * (i - a_dbl);
-      b += 2.0;
-
-      D = b + an * D;
-      if (fabs(D) < EPSILON) {
-        D = EPSILON;
-      }
-      C = b + an / C;
-      if (fabs(C) < EPSILON) {
-        C = EPSILON;
-      }
-
-      D = 1.0 / D;
-      const double delta = C * D;
-      f *= delta;
-
-      const double delta_m1 = fabs(delta - 1.0);
-      if (delta_m1 < precision) {
-        break;
-      }
-    }
-
-    result.log_q = log_prefactor - log(f);
+    result.log_q
+        = internal::log_q_gamma_cf(a_dbl, z_dbl, max_steps, precision);
 
     // For gradient, use: d/da log(Q) = (1/Q) * dQ/da
     // grad_reg_inc_gamma computes dQ/da
diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp
@@ -26,66 +26,6 @@
 namespace stan {
 namespace math {
 
-namespace internal {
-
-/**
- * Compute log(Q(a,x)) using continued fraction expansion for upper incomplete
- * gamma function. When used with fvar types, automatically computes
- * derivatives.
- *
- * @tparam T_a Type of shape parameter a (double or fvar types)
- * @param a Shape parameter
- * @param x Value at which to evaluate
- * @param max_steps Maximum number of continued fraction iterations
- * @param precision Convergence threshold
- * @return log(Q(a,x)) with same type as T_a
- */
-template <typename T_a, typename T_x>
-inline auto log_q_gamma_cf(const T_a& a, const T_x& x, int max_steps = 250,
-                           double precision = 1e-16) {
-  using stan::math::lgamma;
-  using stan::math::log;
-  using stan::math::value_of;
-  using std::fabs;
-  using T_return = return_type_t<T_a, T_x>;
-
-  const T_return a_ret = a;
-  const T_return x_ret = x;
-  const auto log_prefactor = a_ret * log(x_ret) - x_ret - lgamma(a_ret);
-
-  auto b = x_ret + 1.0 - a_ret;
-  auto C = (fabs(value_of(b)) >= EPSILON) ? b : T_return(EPSILON);
-  auto D = T_return(0.0);
-  auto f = C;
-
-  for (int i = 1; i <= max_steps; ++i) {
-    auto an = -i * (i - a_ret);
-    b += 2.0;
-
-    D = b + an * D;
-    if (fabs(value_of(D)) < EPSILON) {
-      D = T_a(EPSILON);
-    }
-    C = b + an / C;
-    if (fabs(value_of(C)) < EPSILON) {
-      C = T_a(EPSILON);
-    }
-
-    D = 1.0 / D;
-    auto delta = C * D;
-    f *= delta;
-
-    const double delta_m1 = value_of(fabs(value_of(delta) - 1.0));
-    if (delta_m1 < precision) {
-      break;
-    }
-  }
-
-  return log_prefactor - log(f);
-}
-
-}  // namespace internal
-
 template <typename T_y, typename T_shape, typename T_inv_scale>
 return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
                                                      const T_shape& alpha,
diff --git a/test/unit/math/prim/prob/gamma_lccdf_test.cpp b/test/unit/math/prim/prob/gamma_lccdf_test.cpp
@@ -91,6 +91,26 @@ TEST(ProbGamma, lccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) {
   EXPECT_DOUBLE_EQ(new_val, expected);
 }
 
+TEST(ProbGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) {
+  using stan::math::gamma_lcdf;
+  using stan::math::gamma_lccdf;
+  using stan::math::log1m_exp;
+  using stan::math::negative_infinity;
+
+  double y = 20000.0;
+  double alpha = 800.0;
+  double beta = 0.1;
+
+  double lcdf = gamma_lcdf(y, alpha, beta);
+  double log1m_lcdf = log1m_exp(lcdf);
+  double lccdf = gamma_lccdf(y, alpha, beta);
+
+  EXPECT_EQ(lcdf, 0.0);
+  EXPECT_EQ(log1m_lcdf, negative_infinity());
+  EXPECT_TRUE(std::isfinite(lccdf));
+  EXPECT_LT(lccdf, 0.0);
+}
+
 TEST(ProbGamma, lccdf_large_alpha_large_y) {
   using stan::math::gamma_lccdf;
 
diff --git a/test/unit/math/rev/prob/gamma_lccdf_test.cpp b/test/unit/math/rev/prob/gamma_lccdf_test.cpp
@@ -272,6 +272,31 @@ TEST(ProbDistributionsGamma,
   EXPECT_LE(grads[0], 0.0);
 }
 
+TEST(ProbDistributionsGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) {
+  using stan::math::gamma_lcdf;
+  using stan::math::gamma_lccdf;
+  using stan::math::log1m_exp;
+  using stan::math::negative_infinity;
+  using stan::math::var;
+
+  double y_d = 20000.0;
+  double alpha_d = 800.0;
+  double beta_d = 0.1;
+
+  double lcdf = gamma_lcdf(y_d, alpha_d, beta_d);
+  double log1m_lcdf = log1m_exp(lcdf);
+
+  var y_v = y_d;
+  var alpha_v = alpha_d;
+  var beta_v = beta_d;
+  var lccdf_var = gamma_lccdf(y_v, alpha_v, beta_v);
+
+  EXPECT_EQ(lcdf, 0.0);
+  EXPECT_EQ(log1m_lcdf, negative_infinity());
+  EXPECT_TRUE(std::isfinite(lccdf_var.val()));
+  EXPECT_LT(lccdf_var.val(), 0.0);
+}
+
 TEST(ProbDistributionsGamma, lccdf_extreme_values_large) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
