From 93d47500967ec7847a3d3fc8576c6cd33aca3ecd Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Thu, 5 Feb 2026 10:28:11 -0500 Subject: [PATCH 01/11] Revert "Merge pull request #3272 from stan-dev/revert-3266-fix-gamma-lccdf-v2" This reverts commit 42063be9833129cdf1d4b2299616ec8aba237d80, reversing changes made to bf9fc6bc00539c5ddb31d4ed622b08b19ea36c44. --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 156 ++++++++++++++++++ stan/math/prim/prob/gamma_lccdf.hpp | 148 ++++++++++++----- test/unit/math/prim/prob/gamma_lccdf_test.cpp | 68 ++++++++ test/unit/math/rev/prob/gamma_lccdf_test.cpp | 97 +++++++++++ 4 files changed, 430 insertions(+), 39 deletions(-) create mode 100644 stan/math/prim/fun/log_gamma_q_dgamma.hpp diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp new file mode 100644 index 00000000000..e264c5d50a5 --- /dev/null +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -0,0 +1,156 @@ +#ifndef STAN_MATH_PRIM_FUN_LOG_GAMMA_Q_DGAMMA_HPP +#define STAN_MATH_PRIM_FUN_LOG_GAMMA_Q_DGAMMA_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Result structure containing log(Q(a,z)) and its gradient with respect to a. + * + * @tparam T return type + */ +template +struct log_gamma_q_result { + T log_q; ///< log(Q(a,z)) where Q is upper regularized incomplete gamma + 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 precision Convergence threshold + * @param max_steps Maximum number of continued fraction iterations + * @return log(Q(a,z)) with the return type of T_a and T_z + */ +template +inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, + double precision = 1e-16, + int max_steps = 250) { + using stan::math::lgamma; + using stan::math::log; + using stan::math::value_of; + using std::fabs; + using T_return = return_type_t; + + const T_return log_prefactor = a * log(z) - z - lgamma(a); + + T_return b = z + 1.0 - a; + T_return C = (fabs(value_of(b)) >= EPSILON) ? b : T_return(EPSILON); + T_return D = 0.0; + T_return f = C; + + for (int i = 1; i <= max_steps; ++i) { + T_a an = -i * (i - a); + b += 2.0; + + D = b + an * D; + if (fabs(D) < EPSILON) { + D = EPSILON; + } + C = b + an / C; + if (fabs(C) < EPSILON) { + C = EPSILON; + } + + D = inv(D); + T_return 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 + * upper incomplete gamma function. + * + * This uses a continued fraction representation for numerical stability when + * computing the upper incomplete gamma function in log space, along with + * analytical gradient computation. + * + * @tparam T_a type of the shape parameter + * @tparam T_z type of the value parameter + * @param a shape parameter (must be positive) + * @param z value parameter (must be non-negative) + * @param precision convergence threshold + * @param max_steps maximum iterations for continued fraction + * @return structure containing log(Q(a,z)) and d/da log(Q(a,z)) + */ +template +inline log_gamma_q_result> log_gamma_q_dgamma( + const T_a& a, const T_z& z, double precision = 1e-16, int max_steps = 250) { + using std::exp; + using std::log; + using T_return = return_type_t; + + const double a_dbl = value_of(a); + const double z_dbl = value_of(z); + + log_gamma_q_result result; + + // For z > a + 1, use continued fraction for better numerical stability + if (z_dbl > a_dbl + 1.0) { + result.log_q = internal::log_q_gamma_cf(a_dbl, z_dbl, precision, max_steps); + + // For gradient, use: d/da log(Q) = (1/Q) * dQ/da + // grad_reg_inc_gamma computes dQ/da + const double Q_val = exp(result.log_q); + const double dQ_da + = grad_reg_inc_gamma(a_dbl, z_dbl, tgamma(a_dbl), digamma(a_dbl)); + result.dlog_q_da = dQ_da / Q_val; + + } else { + // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy + const double P_val = gamma_p(a_dbl, z_dbl); + result.log_q = log1m(P_val); + + // Gradient: d/da log(Q) = (1/Q) * dQ/da + // grad_reg_inc_gamma computes dQ/da + const double Q_val = exp(result.log_q); + if (Q_val > 0) { + const double dQ_da + = grad_reg_inc_gamma(a_dbl, z_dbl, tgamma(a_dbl), digamma(a_dbl)); + result.dlog_q_da = dQ_da / Q_val; + } else { + // Fallback if Q rounds to zero - use asymptotic approximation + result.dlog_q_da = log(z_dbl) - digamma(a_dbl); + } + } + + return result; +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index a670cefcecf..eada09359e9 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -6,15 +6,20 @@ #include #include #include -#include +#include +#include #include +#include +#include #include +#include #include #include #include #include #include -#include +#include +#include #include #include @@ -24,10 +29,9 @@ namespace math { template inline return_type_t gamma_lccdf( const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { - using T_partials_return = partials_return_t; using std::exp; using std::log; - using std::pow; + using T_partials_return = partials_return_t; using T_y_ref = ref_type_t; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; @@ -51,61 +55,127 @@ inline return_type_t gamma_lccdf( scalar_seq_view y_vec(y_ref); scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); - size_t N = max_size(y, alpha, beta); - - // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(y); i++) { - if (y_vec.val(i) == 0) { - // LCCDF(0) = log(P(Y > 0)) = log(1) = 0 - return ops_partials.build(0.0); - } - } + const size_t N = max_size(y, alpha, beta); + + constexpr bool any_fvar = is_fvar>::value + || is_fvar>::value + || is_fvar>::value; + constexpr bool partials_fvar = is_fvar::value; for (size_t n = 0; n < N; n++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { - // LCCDF(∞) = log(P(Y > ∞)) = log(0) = -∞ + const T_partials_return y_dbl = y_vec.val(n); + if (y_dbl == 0.0) { + continue; + } + if (y_dbl == INFTY) { return ops_partials.build(negative_infinity()); } - const T_partials_return y_dbl = y_vec.val(n); const T_partials_return alpha_dbl = alpha_vec.val(n); const T_partials_return beta_dbl = beta_vec.val(n); - const T_partials_return beta_y_dbl = beta_dbl * y_dbl; - // Qn = 1 - Pn - const T_partials_return Qn = gamma_q(alpha_dbl, beta_y_dbl); - const T_partials_return log_Qn = log(Qn); + const T_partials_return beta_y = beta_dbl * y_dbl; + if (beta_y == INFTY) { + return ops_partials.build(negative_infinity()); + } + bool use_cf = beta_y > alpha_dbl + 1.0; + T_partials_return log_Qn; + [[maybe_unused]] T_partials_return dlogQ_dalpha = 0.0; + + // Branch by autodiff type first, then handle use_cf logic inside each path + if constexpr (!any_fvar && is_autodiff_v) { + // var-only path: use log_gamma_q_dgamma which computes both log_q + // and its gradient analytically with double inputs + const double beta_y_dbl = value_of_rec(beta_y); + const double alpha_dbl_val = value_of_rec(alpha_dbl); + + if (use_cf) { + auto log_q_result = log_gamma_q_dgamma(alpha_dbl_val, beta_y_dbl); + log_Qn = log_q_result.log_q; + dlogQ_dalpha = log_q_result.dlog_q_da; + } else { + const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); + log_Qn = log1m(Pn); + const T_partials_return Qn = exp(log_Qn); + + // Check if we need to fallback to continued fraction + bool need_cf_fallback + = !std::isfinite(value_of_rec(log_Qn)) || Qn <= 0.0; + if (need_cf_fallback && beta_y > 0.0) { + auto log_q_result = log_gamma_q_dgamma(alpha_dbl_val, beta_y_dbl); + log_Qn = log_q_result.log_q; + dlogQ_dalpha = log_q_result.dlog_q_da; + } else { + dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha_dbl, beta_y) / Qn; + } + } + } else if constexpr (partials_fvar && is_autodiff_v) { + // fvar path: use unit derivative trick to compute gradients + T_partials_return alpha_unit = alpha_dbl; + alpha_unit.d_ = 1; + T_partials_return beta_unit = beta_y; + beta_unit.d_ = 0; + + if (use_cf) { + log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); + T_partials_return log_Qn_fvar + = internal::log_q_gamma_cf(alpha_unit, beta_unit); + dlogQ_dalpha = log_Qn_fvar.d_; + } else { + const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); + log_Qn = log1m(Pn); + + if (!std::isfinite(value_of_rec(log_Qn)) && beta_y > 0.0) { + // Fallback to continued fraction + log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); + T_partials_return log_Qn_fvar + = internal::log_q_gamma_cf(alpha_unit, beta_unit); + dlogQ_dalpha = log_Qn_fvar.d_; + } else { + T_partials_return log_Qn_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + dlogQ_dalpha = log_Qn_fvar.d_; + } + } + } else { + // No alpha derivative needed (alpha is constant or double-only) + if (use_cf) { + log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); + } else { + const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); + log_Qn = log1m(Pn); + + if (!std::isfinite(value_of_rec(log_Qn)) && beta_y > 0.0) { + log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); + } + } + } + if (!std::isfinite(value_of_rec(log_Qn))) { + return ops_partials.build(negative_infinity()); + } P += log_Qn; - if constexpr (is_any_autodiff_v) { - const T_partials_return log_y_dbl = log(y_dbl); - const T_partials_return log_beta_dbl = log(beta_dbl); - const T_partials_return log_pdf - = alpha_dbl * log_beta_dbl - lgamma(alpha_dbl) - + (alpha_dbl - 1.0) * log_y_dbl - beta_y_dbl; - const T_partials_return common_term = exp(log_pdf - log_Qn); + if constexpr (is_autodiff_v || is_autodiff_v) { + const T_partials_return log_y = log(y_dbl); + const T_partials_return alpha_minus_one = fma(alpha_dbl, log_y, -log_y); + + const T_partials_return log_pdf = alpha_dbl * log(beta_dbl) + - lgamma(alpha_dbl) + alpha_minus_one + - beta_y; + + const T_partials_return hazard = exp(log_pdf - log_Qn); // f/Q if constexpr (is_autodiff_v) { - // d/dy log(1-F(y)) = -f(y)/(1-F(y)) - partials<0>(ops_partials)[n] -= common_term; + partials<0>(ops_partials)[n] -= hazard; } if constexpr (is_autodiff_v) { - // d/dbeta log(1-F(y)) = -y*f(y)/(beta*(1-F(y))) - partials<2>(ops_partials)[n] -= y_dbl / beta_dbl * common_term; + partials<2>(ops_partials)[n] -= (y_dbl / beta_dbl) * hazard; } } - if constexpr (is_autodiff_v) { - const T_partials_return digamma_val = digamma(alpha_dbl); - const T_partials_return gamma_val = tgamma(alpha_dbl); - // d/dalpha log(1-F(y)) = grad_upper_inc_gamma / (1-F(y)) - partials<1>(ops_partials)[n] - += grad_reg_inc_gamma(alpha_dbl, beta_y_dbl, gamma_val, digamma_val) - / Qn; + partials<1>(ops_partials)[n] += dlogQ_dalpha; } } return ops_partials.build(P); diff --git a/test/unit/math/prim/prob/gamma_lccdf_test.cpp b/test/unit/math/prim/prob/gamma_lccdf_test.cpp index 2893f2f0166..e0c84861e0c 100644 --- a/test/unit/math/prim/prob/gamma_lccdf_test.cpp +++ b/test/unit/math/prim/prob/gamma_lccdf_test.cpp @@ -66,6 +66,51 @@ TEST(ProbGamma, lccdf_small_alpha_small_y) { EXPECT_LT(result, 0.0); } +TEST(ProbGamma, lccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) { + using stan::math::gamma_lccdf; + using stan::math::gamma_p; + using stan::math::gamma_q; + using stan::math::log1m; + + // For large alpha and very small y, the CCDF is extremely close to 1. + // The old implementation computed `log(gamma_q(alpha, beta * y))`, which can + // round to `log(1) == 0`. The updated implementation uses `log1m(gamma_p)`, + // which preserves the tiny negative value. + double y = 1e-8; + double alpha = 31.25; + double beta = 1.0; + + double new_val = gamma_lccdf(y, alpha, beta); + double expected = log1m(gamma_p(alpha, beta * y)); + + // Old code: log(gamma_q(alpha, beta * y)) + double old_val = std::log(gamma_q(alpha, beta * y)); + + EXPECT_EQ(old_val, 0.0); + EXPECT_LT(new_val, 0.0); + EXPECT_DOUBLE_EQ(new_val, expected); +} + +TEST(ProbGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) { + using stan::math::gamma_lccdf; + using stan::math::gamma_lcdf; + 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; @@ -154,6 +199,29 @@ TEST(ProbGamma, lccdf_extreme_large_alpha) { EXPECT_TRUE(std::isfinite(result)); } +TEST(ProbGamma, lccdf_large_alpha_1000_beta_3) { + using stan::math::gamma_lccdf; + + // Large alpha = 1000, beta = 3 + double alpha = 1000.0; + double beta = 3.0; + + // Test various y values + std::vector y_values = {100.0, 300.0, 333.333, 400.0, 500.0}; + + for (double y : y_values) { + double result = gamma_lccdf(y, alpha, beta); + + // Result should be finite + EXPECT_TRUE(std::isfinite(result)) + << "Failed for y=" << y << ", alpha=" << alpha << ", beta=" << beta; + + // Result should be <= 0 (log of probability) + EXPECT_LE(result, 0.0) << "Positive value for y=" << y + << ", alpha=" << alpha << ", beta=" << beta; + } +} + TEST(ProbGamma, lccdf_monotonic_in_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 index 1066773e060..aa1f9a1a942 100644 --- a/test/unit/math/rev/prob/gamma_lccdf_test.cpp +++ b/test/unit/math/rev/prob/gamma_lccdf_test.cpp @@ -230,6 +230,73 @@ TEST(ProbDistributionsGamma, lccdf_extreme_values_small) { } } +TEST(ProbDistributionsGamma, + lccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) { + using stan::math::gamma_lccdf; + using stan::math::gamma_p; + using stan::math::gamma_q; + using stan::math::log1m; + using stan::math::var; + + // Same comparison as the prim test, but also exercises autodiff for + // alpha > 30. + double y_d = 1e-8; + double alpha_d = 31.25; + double beta_d = 1.0; + + 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); + + // Old code: log(gamma_q(alpha, beta * y)) + double old_val = std::log(gamma_q(alpha_d, beta_d * y_d)); + double expected = log1m(gamma_p(alpha_d, beta_d * y_d)); + + EXPECT_EQ(old_val, 0.0); + EXPECT_LT(lccdf_var.val(), 0.0); + EXPECT_DOUBLE_EQ(lccdf_var.val(), expected); + + std::vector vars = {y_v, alpha_v, beta_v}; + std::vector grads; + lccdf_var.grad(vars, grads); + + for (size_t i = 0; i < grads.size(); ++i) { + EXPECT_FALSE(std::isnan(grads[i])) << "Gradient " << i << " is NaN"; + EXPECT_TRUE(std::isfinite(grads[i])) + << "Gradient " << i << " is not finite"; + } + + // d/dy log(CCDF) should be <= 0 (can underflow to -0) + EXPECT_LE(grads[0], 0.0); +} + +TEST(ProbDistributionsGamma, lccdf_log1m_exp_lcdf_rounds_to_inf) { + using stan::math::gamma_lccdf; + using stan::math::gamma_lcdf; + 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; @@ -258,6 +325,36 @@ TEST(ProbDistributionsGamma, lccdf_extreme_values_large) { } } +TEST(ProbDistributionsGamma, lccdf_large_alpha_1000_beta_3) { + using stan::math::gamma_lccdf; + using stan::math::var; + + // Large alpha = 1000, beta = 3 + // Note: This test only checks values, not gradients, as large alpha values + // can cause numerical issues with gradient computation + double alpha_d = 1000.0; + double beta_d = 3.0; + + // Test various y values + std::vector y_values = {100.0, 300.0, 333.333, 400.0, 500.0}; + + for (double y_d : y_values) { + 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); + + // Value should be finite and <= 0 + EXPECT_TRUE(std::isfinite(lccdf_var.val())) + << "Failed for y=" << y_d << ", alpha=" << alpha_d + << ", beta=" << beta_d; + EXPECT_LE(lccdf_var.val(), 0.0) + << "Positive value for y=" << y_d << ", alpha=" << alpha_d + << ", beta=" << beta_d; + } +} + TEST(ProbDistributionsGamma, lccdf_alpha_one_derivatives) { using stan::math::gamma_lccdf; using stan::math::var; From e9a3c56b523683dfea32687f62e7adfeede5c63e Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Tue, 13 Jan 2026 15:55:27 -0500 Subject: [PATCH 02/11] refactor gamma_lccdf --- stan/math/prim/prob/gamma_lccdf.hpp | 171 +++++++++++++++------------- 1 file changed, 95 insertions(+), 76 deletions(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index eada09359e9..aab4590dd96 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -25,10 +25,85 @@ namespace stan { namespace math { +namespace internal { + template + struct Q_eval { + T log_Q{0.0}; + T dlogQ_dalpha{0.0}; + bool ok{false}; + }; + + /** + * Computes log q and d(log q) / d(alpha) using continued fraction. + */ + template + static inline Q_eval eval_q_cf(const T& alpha, + const T& beta_y) { + Q_eval out; + if constexpr (!any_fvar && is_autodiff_v) { + auto log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); + out.log_Q = log_q_result.log_q; + out.dlogQ_dalpha = log_q_result.dlog_q_da; + } else { + out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); + if constexpr (is_autodiff_v) { + if constexpr (!partials_fvar) { + out.dlogQ_dalpha + = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), + digamma(alpha)) / exp(out.log_Q); + } else { + T alpha_unit = alpha; + alpha_unit.d_ = 1; + T beta_y_unit = beta_y; + beta_y_unit.d_ = 0; + T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); + out.dlogQ_dalpha = log_Q_fvar.d_; + } + } + } + + out.ok = std::isfinite(value_of_rec(out.log_Q)); + return out; + } + + /** + * Computes log q and d(log q) / d(alpha) using log1m. + */ + template + static inline Q_eval eval_q_log1m(const T& alpha, + const T& beta_y) { + Q_eval out; + out.log_Q = log1m(gamma_p(alpha, beta_y)); + + if (!std::isfinite(value_of_rec(out.log_Q))) { + out.ok = false; + return out; + } + + if constexpr (is_autodiff_v) { + if constexpr (partials_fvar) { + T alpha_unit = alpha; + alpha_unit.d_ = 1; + T beta_unit = beta_y; + beta_unit.d_ = 0; + T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + out.dlogQ_dalpha = log_Q_fvar.d_; + } else { + out.dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); + } + } + + out.ok = true; + return out; + } +} template -inline return_type_t gamma_lccdf( - const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { +inline return_type_t gamma_lccdf(const T_y& y, + const T_shape& alpha, + const T_inv_scale& beta) { using std::exp; using std::log; using T_partials_return = partials_return_t; @@ -81,91 +156,35 @@ inline return_type_t gamma_lccdf( return ops_partials.build(negative_infinity()); } - bool use_cf = beta_y > alpha_dbl + 1.0; - T_partials_return log_Qn; - [[maybe_unused]] T_partials_return dlogQ_dalpha = 0.0; - - // Branch by autodiff type first, then handle use_cf logic inside each path - if constexpr (!any_fvar && is_autodiff_v) { - // var-only path: use log_gamma_q_dgamma which computes both log_q - // and its gradient analytically with double inputs - const double beta_y_dbl = value_of_rec(beta_y); - const double alpha_dbl_val = value_of_rec(alpha_dbl); - - if (use_cf) { - auto log_q_result = log_gamma_q_dgamma(alpha_dbl_val, beta_y_dbl); - log_Qn = log_q_result.log_q; - dlogQ_dalpha = log_q_result.dlog_q_da; - } else { - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); - log_Qn = log1m(Pn); - const T_partials_return Qn = exp(log_Qn); - - // Check if we need to fallback to continued fraction - bool need_cf_fallback - = !std::isfinite(value_of_rec(log_Qn)) || Qn <= 0.0; - if (need_cf_fallback && beta_y > 0.0) { - auto log_q_result = log_gamma_q_dgamma(alpha_dbl_val, beta_y_dbl); - log_Qn = log_q_result.log_q; - dlogQ_dalpha = log_q_result.dlog_q_da; - } else { - dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha_dbl, beta_y) / Qn; - } - } - } else if constexpr (partials_fvar && is_autodiff_v) { - // fvar path: use unit derivative trick to compute gradients - T_partials_return alpha_unit = alpha_dbl; - alpha_unit.d_ = 1; - T_partials_return beta_unit = beta_y; - beta_unit.d_ = 0; - - if (use_cf) { - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - T_partials_return log_Qn_fvar - = internal::log_q_gamma_cf(alpha_unit, beta_unit); - dlogQ_dalpha = log_Qn_fvar.d_; - } else { - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); - log_Qn = log1m(Pn); - - if (!std::isfinite(value_of_rec(log_Qn)) && beta_y > 0.0) { - // Fallback to continued fraction - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - T_partials_return log_Qn_fvar - = internal::log_q_gamma_cf(alpha_unit, beta_unit); - dlogQ_dalpha = log_Qn_fvar.d_; - } else { - T_partials_return log_Qn_fvar = log1m(gamma_p(alpha_unit, beta_unit)); - dlogQ_dalpha = log_Qn_fvar.d_; - } - } + const bool use_continued_fraction = beta_y > alpha_dbl + 1.0; + internal::Q_eval result; + if (use_continued_fraction) { + result = internal::eval_q_cf(alpha_dbl, beta_y); } else { - // No alpha derivative needed (alpha is constant or double-only) - if (use_cf) { - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - } else { - const T_partials_return Pn = gamma_p(alpha_dbl, beta_y); - log_Qn = log1m(Pn); + result = internal::eval_q_log1m(alpha_dbl, beta_y); - if (!std::isfinite(value_of_rec(log_Qn)) && beta_y > 0.0) { - log_Qn = internal::log_q_gamma_cf(alpha_dbl, beta_y); - } + if (!result.ok && beta_y > 0.0) { + // Fallback to continued fraction if log1m fails + result = internal::eval_q_cf(alpha_dbl, beta_y); } } - if (!std::isfinite(value_of_rec(log_Qn))) { + if (!result.ok) { return ops_partials.build(negative_infinity()); } - P += log_Qn; + + P += result.log_Q; if constexpr (is_autodiff_v || is_autodiff_v) { const T_partials_return log_y = log(y_dbl); const T_partials_return alpha_minus_one = fma(alpha_dbl, log_y, -log_y); - const T_partials_return log_pdf = alpha_dbl * log(beta_dbl) - - lgamma(alpha_dbl) + alpha_minus_one - - beta_y; + const T_partials_return log_pdf + = alpha_dbl * log(beta_dbl) - lgamma(alpha_dbl) + alpha_minus_one - beta_y; - const T_partials_return hazard = exp(log_pdf - log_Qn); // f/Q + const T_partials_return hazard = exp(log_pdf - result.log_Q); // f/Q if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] -= hazard; @@ -175,7 +194,7 @@ inline return_type_t gamma_lccdf( } } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += dlogQ_dalpha; + partials<1>(ops_partials)[n] += result.dlogQ_dalpha; } } return ops_partials.build(P); From 86ac5617162c98dcad1f1bb6fc240792db23185a Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 5 Feb 2026 10:36:44 -0500 Subject: [PATCH 03/11] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/gamma_lccdf.hpp | 147 ++++++++++++++-------------- 1 file changed, 73 insertions(+), 74 deletions(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index aab4590dd96..6c7a0ebd9d2 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -26,84 +26,81 @@ namespace stan { namespace math { namespace internal { - template - struct Q_eval { - T log_Q{0.0}; - T dlogQ_dalpha{0.0}; - bool ok{false}; - }; - - /** - * Computes log q and d(log q) / d(alpha) using continued fraction. - */ - template - static inline Q_eval eval_q_cf(const T& alpha, - const T& beta_y) { - Q_eval out; - if constexpr (!any_fvar && is_autodiff_v) { - auto log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); - out.log_Q = log_q_result.log_q; - out.dlogQ_dalpha = log_q_result.dlog_q_da; - } else { - out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); - if constexpr (is_autodiff_v) { - if constexpr (!partials_fvar) { - out.dlogQ_dalpha - = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), - digamma(alpha)) / exp(out.log_Q); - } else { - T alpha_unit = alpha; - alpha_unit.d_ = 1; - T beta_y_unit = beta_y; - beta_y_unit.d_ = 0; - T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); - out.dlogQ_dalpha = log_Q_fvar.d_; - } - } - } - - out.ok = std::isfinite(value_of_rec(out.log_Q)); - return out; - } - - /** - * Computes log q and d(log q) / d(alpha) using log1m. - */ - template - static inline Q_eval eval_q_log1m(const T& alpha, - const T& beta_y) { - Q_eval out; - out.log_Q = log1m(gamma_p(alpha, beta_y)); - - if (!std::isfinite(value_of_rec(out.log_Q))) { - out.ok = false; - return out; - } - +template +struct Q_eval { + T log_Q{0.0}; + T dlogQ_dalpha{0.0}; + bool ok{false}; +}; + +/** + * Computes log q and d(log q) / d(alpha) using continued fraction. + */ +template +static inline Q_eval eval_q_cf(const T& alpha, const T& beta_y) { + Q_eval out; + if constexpr (!any_fvar && is_autodiff_v) { + auto log_q_result + = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); + out.log_Q = log_q_result.log_q; + out.dlogQ_dalpha = log_q_result.dlog_q_da; + } else { + out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); if constexpr (is_autodiff_v) { - if constexpr (partials_fvar) { + if constexpr (!partials_fvar) { + out.dlogQ_dalpha + = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), digamma(alpha)) + / exp(out.log_Q); + } else { T alpha_unit = alpha; alpha_unit.d_ = 1; - T beta_unit = beta_y; - beta_unit.d_ = 0; - T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + T beta_y_unit = beta_y; + beta_y_unit.d_ = 0; + T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); out.dlogQ_dalpha = log_Q_fvar.d_; - } else { - out.dlogQ_dalpha = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); } } + } + + out.ok = std::isfinite(value_of_rec(out.log_Q)); + return out; +} - out.ok = true; +/** + * Computes log q and d(log q) / d(alpha) using log1m. + */ +template +static inline Q_eval eval_q_log1m(const T& alpha, const T& beta_y) { + Q_eval out; + out.log_Q = log1m(gamma_p(alpha, beta_y)); + + if (!std::isfinite(value_of_rec(out.log_Q))) { + out.ok = false; return out; } + + if constexpr (is_autodiff_v) { + if constexpr (partials_fvar) { + T alpha_unit = alpha; + alpha_unit.d_ = 1; + T beta_unit = beta_y; + beta_unit.d_ = 0; + T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + out.dlogQ_dalpha = log_Q_fvar.d_; + } else { + out.dlogQ_dalpha + = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); + } + } + + out.ok = true; + return out; } +} // namespace internal template -inline return_type_t gamma_lccdf(const T_y& y, - const T_shape& alpha, - const T_inv_scale& beta) { +inline return_type_t gamma_lccdf( + const T_y& y, const T_shape& alpha, const T_inv_scale& beta) { using std::exp; using std::log; using T_partials_return = partials_return_t; @@ -159,16 +156,17 @@ inline return_type_t gamma_lccdf(const T_y& y, const bool use_continued_fraction = beta_y > alpha_dbl + 1.0; internal::Q_eval result; if (use_continued_fraction) { - result = internal::eval_q_cf(alpha_dbl, beta_y); - } else { - result = internal::eval_q_log1m(alpha_dbl, beta_y); + } else { + result + = internal::eval_q_log1m( + alpha_dbl, beta_y); if (!result.ok && beta_y > 0.0) { // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf(alpha_dbl, beta_y); + result = internal::eval_q_cf(alpha_dbl, beta_y); } } if (!result.ok) { @@ -181,8 +179,9 @@ inline return_type_t gamma_lccdf(const T_y& y, const T_partials_return log_y = log(y_dbl); const T_partials_return alpha_minus_one = fma(alpha_dbl, log_y, -log_y); - const T_partials_return log_pdf - = alpha_dbl * log(beta_dbl) - lgamma(alpha_dbl) + alpha_minus_one - beta_y; + const T_partials_return log_pdf = alpha_dbl * log(beta_dbl) + - lgamma(alpha_dbl) + alpha_minus_one + - beta_y; const T_partials_return hazard = exp(log_pdf - result.log_Q); // f/Q From c150da9bd6bddb7806127403ae2772dc0548e2b7 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 4 Mar 2026 11:34:30 -0500 Subject: [PATCH 04/11] cleanup for gamma_lccdf --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 50 ++++++++--------------- stan/math/prim/prob/gamma_lccdf.hpp | 46 +++++++++------------ 2 files changed, 37 insertions(+), 59 deletions(-) diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index e264c5d50a5..b2fc239f908 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -48,18 +49,13 @@ template inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, double precision = 1e-16, int max_steps = 250) { - using stan::math::lgamma; - using stan::math::log; - using stan::math::value_of; - using std::fabs; using T_return = return_type_t; + const auto log_prefactor = a * log(z) - z - lgamma(a); - const T_return log_prefactor = a * log(z) - z - lgamma(a); - - T_return b = z + 1.0 - a; - T_return C = (fabs(value_of(b)) >= EPSILON) ? b : T_return(EPSILON); - T_return D = 0.0; - T_return f = C; + auto b = z + 1.0 - a; + auto C = (fabs(value_of(b)) >= EPSILON) ? b : std::decay_t(EPSILON); + auto D = 0.0; + auto f = C; for (int i = 1; i <= max_steps; ++i) { T_a an = -i * (i - a); @@ -75,7 +71,7 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, } D = inv(D); - T_return delta = C * D; + auto delta = C * D; f *= delta; const double delta_m1 = value_of(fabs(value_of(delta) - 1.0)); @@ -83,7 +79,6 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, break; } } - return log_prefactor - log(f); } @@ -109,45 +104,36 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, template inline log_gamma_q_result> log_gamma_q_dgamma( const T_a& a, const T_z& z, double precision = 1e-16, int max_steps = 250) { - using std::exp; - using std::log; using T_return = return_type_t; - - const double a_dbl = value_of(a); - const double z_dbl = value_of(z); - - log_gamma_q_result result; - + const auto a_dbl = value_of(a); + const auto z_dbl = value_of(z); // For z > a + 1, use continued fraction for better numerical stability if (z_dbl > a_dbl + 1.0) { - result.log_q = internal::log_q_gamma_cf(a_dbl, z_dbl, precision, max_steps); - + log_gamma_q_result result{internal::log_q_gamma_cf(a_dbl, z_dbl, precision, max_steps), 0.0}; // For gradient, use: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da - const double Q_val = exp(result.log_q); - const double dQ_da + const auto Q_val = exp(result.log_q); + const auto dQ_da = grad_reg_inc_gamma(a_dbl, z_dbl, tgamma(a_dbl), digamma(a_dbl)); result.dlog_q_da = dQ_da / Q_val; - + return result; } else { // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy - const double P_val = gamma_p(a_dbl, z_dbl); - result.log_q = log1m(P_val); - + const auto P_val = gamma_p(a_dbl, z_dbl); + log_gamma_q_result result{log1m(P_val), 0.0}; // Gradient: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da - const double Q_val = exp(result.log_q); + const auto Q_val = exp(result.log_q); if (Q_val > 0) { - const double dQ_da + const auto dQ_da = grad_reg_inc_gamma(a_dbl, z_dbl, tgamma(a_dbl), digamma(a_dbl)); result.dlog_q_da = dQ_da / Q_val; } else { // Fallback if Q rounds to zero - use asymptotic approximation result.dlog_q_da = log(z_dbl) - digamma(a_dbl); } + return result; } - - return result; } } // namespace math diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index 6c7a0ebd9d2..ad6b14ec15f 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -22,7 +22,7 @@ #include #include #include - +#include namespace stan { namespace math { namespace internal { @@ -30,15 +30,14 @@ template struct Q_eval { T log_Q{0.0}; T dlogQ_dalpha{0.0}; - bool ok{false}; }; /** * Computes log q and d(log q) / d(alpha) using continued fraction. */ -template -static inline Q_eval eval_q_cf(const T& alpha, const T& beta_y) { - Q_eval out; +template +static inline auto eval_q_cf(const T& alpha, const T& beta_y) { + Q_eval> out; if constexpr (!any_fvar && is_autodiff_v) { auto log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); @@ -62,23 +61,23 @@ static inline Q_eval eval_q_cf(const T& alpha, const T& beta_y) { } } - out.ok = std::isfinite(value_of_rec(out.log_Q)); - return out; + if (std::isfinite(value_of_rec(out.log_Q))) { + return std::optional>>{out}; + } else { + return std::optional>>{std::nullopt}; + } } /** * Computes log q and d(log q) / d(alpha) using log1m. */ -template +template static inline Q_eval eval_q_log1m(const T& alpha, const T& beta_y) { Q_eval out; out.log_Q = log1m(gamma_p(alpha, beta_y)); - if (!std::isfinite(value_of_rec(out.log_Q))) { - out.ok = false; return out; } - if constexpr (is_autodiff_v) { if constexpr (partials_fvar) { T alpha_unit = alpha; @@ -92,8 +91,6 @@ static inline Q_eval eval_q_log1m(const T& alpha, const T& beta_y) { = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); } } - - out.ok = true; return out; } } // namespace internal @@ -154,26 +151,21 @@ inline return_type_t gamma_lccdf( } const bool use_continued_fraction = beta_y > alpha_dbl + 1.0; - internal::Q_eval result; + std::optional> result; if (use_continued_fraction) { - result = internal::eval_q_cf(alpha_dbl, beta_y); + result = internal::eval_q_cf(alpha_dbl, beta_y); } else { - result - = internal::eval_q_log1m( - alpha_dbl, beta_y); - - if (!result.ok && beta_y > 0.0) { + result = internal::eval_q_log1m(alpha_dbl, beta_y); + if (!result && beta_y > 0.0) { // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf(alpha_dbl, beta_y); + result = internal::eval_q_cf(alpha_dbl, beta_y); } } - if (!result.ok) { + if (!result) { return ops_partials.build(negative_infinity()); } - P += result.log_Q; + P += result->log_Q; if constexpr (is_autodiff_v || is_autodiff_v) { const T_partials_return log_y = log(y_dbl); @@ -183,7 +175,7 @@ inline return_type_t gamma_lccdf( - lgamma(alpha_dbl) + alpha_minus_one - beta_y; - const T_partials_return hazard = exp(log_pdf - result.log_Q); // f/Q + const T_partials_return hazard = exp(log_pdf - result->log_Q); // f/Q if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] -= hazard; @@ -193,7 +185,7 @@ inline return_type_t gamma_lccdf( } } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += result.dlogQ_dalpha; + partials<1>(ops_partials)[n] += result->dlogQ_dalpha; } } return ops_partials.build(P); From 70c7601765b89a2c921c60b4f4a3465430367fc5 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 4 Mar 2026 12:38:16 -0500 Subject: [PATCH 05/11] more cleanup --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 67 +++++-------- stan/math/prim/prob/gamma_lccdf.hpp | 112 +++++++++++----------- 2 files changed, 78 insertions(+), 101 deletions(-) diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index b2fc239f908..c9c110ae56c 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -20,17 +20,6 @@ namespace stan { namespace math { -/** - * Result structure containing log(Q(a,z)) and its gradient with respect to a. - * - * @tparam T return type - */ -template -struct log_gamma_q_result { - T log_q; ///< log(Q(a,z)) where Q is upper regularized incomplete gamma - T dlog_q_da; ///< d/da log(Q(a,z)) -}; - namespace internal { /** @@ -41,40 +30,32 @@ namespace internal { * @tparam T_z Type of value parameter z (double or fvar types) * @param a Shape parameter * @param z Value at which to evaluate - * @param precision Convergence threshold + * @param precision Convergence threshold, default of sqrt(machine_epsilon) * @param max_steps Maximum number of continued fraction iterations * @return log(Q(a,z)) with the return type of T_a and T_z */ template inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, - double precision = 1e-16, + double precision = 1.49012e-08, int max_steps = 250) { using T_return = return_type_t; const auto log_prefactor = a * log(z) - z - lgamma(a); - auto b = z + 1.0 - a; - auto C = (fabs(value_of(b)) >= EPSILON) ? b : std::decay_t(EPSILON); + auto b_init = z + 1.0 - a; + auto C = (fabs(value_of(b_init)) >= EPSILON) ? b_init : std::decay_t(EPSILON); auto D = 0.0; auto f = C; - for (int i = 1; i <= max_steps; ++i) { T_a an = -i * (i - a); - b += 2.0; - + const auto b = b_init + 2.0 * i; D = b + an * D; - if (fabs(D) < EPSILON) { - D = EPSILON; - } + D = (fabs(value_of(D)) >= EPSILON) ? D : std::decay_t(EPSILON); C = b + an / C; - if (fabs(C) < EPSILON) { - C = EPSILON; - } - + C = (fabs(value_of(C)) >= EPSILON) ? C : std::decay_t(EPSILON); D = inv(D); auto delta = C * D; f *= delta; - - const double delta_m1 = value_of(fabs(value_of(delta) - 1.0)); + const auto delta_m1 = fabs(value_of(delta) - 1.0); if (delta_m1 < precision) { break; } @@ -97,40 +78,40 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, * @tparam T_z type of the value parameter * @param a shape parameter (must be positive) * @param z value parameter (must be non-negative) - * @param precision convergence threshold + * @param precision convergence threshold, default of sqrt(machine_epsilon) * @param max_steps maximum iterations for continued fraction * @return structure containing log(Q(a,z)) and d/da log(Q(a,z)) */ template -inline log_gamma_q_result> log_gamma_q_dgamma( - const T_a& a, const T_z& z, double precision = 1e-16, int max_steps = 250) { +inline auto log_gamma_q_dgamma( + const T_a& a, const T_z& z, double precision = 1.49012e-08, int max_steps = 250) { using T_return = return_type_t; - const auto a_dbl = value_of(a); - const auto z_dbl = value_of(z); + const auto a_val = value_of(a); + const auto z_val = value_of(z); // For z > a + 1, use continued fraction for better numerical stability - if (z_dbl > a_dbl + 1.0) { - log_gamma_q_result result{internal::log_q_gamma_cf(a_dbl, z_dbl, precision, max_steps), 0.0}; + if (z_val > a_val + 1.0) { + std::pair result{internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0}; // For gradient, use: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da - const auto Q_val = exp(result.log_q); + const auto Q_val = exp(result.first); const auto dQ_da - = grad_reg_inc_gamma(a_dbl, z_dbl, tgamma(a_dbl), digamma(a_dbl)); - result.dlog_q_da = dQ_da / Q_val; + = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val)); + result.second = dQ_da / Q_val; return result; } else { // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy - const auto P_val = gamma_p(a_dbl, z_dbl); - log_gamma_q_result result{log1m(P_val), 0.0}; + const auto P_val = gamma_p(a_val, z_val); + std::pair result{log1m(P_val), 0.0}; // Gradient: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da - const auto Q_val = exp(result.log_q); + const auto Q_val = exp(result.first); if (Q_val > 0) { const auto dQ_da - = grad_reg_inc_gamma(a_dbl, z_dbl, tgamma(a_dbl), digamma(a_dbl)); - result.dlog_q_da = dQ_da / Q_val; + = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val)); + result.second = dQ_da / Q_val; } else { // Fallback if Q rounds to zero - use asymptotic approximation - result.dlog_q_da = log(z_dbl) - digamma(a_dbl); + result.second = log(z_val) - digamma(a_val); } return result; } diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index ad6b14ec15f..97df0e29fbf 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -26,72 +26,70 @@ namespace stan { namespace math { namespace internal { -template -struct Q_eval { - T log_Q{0.0}; - T dlogQ_dalpha{0.0}; -}; /** * Computes log q and d(log q) / d(alpha) using continued fraction. */ -template -static inline auto eval_q_cf(const T& alpha, const T& beta_y) { - Q_eval> out; +template +inline auto eval_q_cf(const T1& alpha, const T2& beta_y) { + using scalar_t = return_type_t; + using ret_t = std::pair; if constexpr (!any_fvar && is_autodiff_v) { auto log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); - out.log_Q = log_q_result.log_q; - out.dlogQ_dalpha = log_q_result.dlog_q_da; + if (likely(std::isfinite(value_of_rec(log_q_result.first)))) { + return std::optional{log_q_result}; + } else { + return std::optional{std::nullopt}; + } } else { - out.log_Q = internal::log_q_gamma_cf(alpha, beta_y); + ret_t out{internal::log_q_gamma_cf(alpha, beta_y), 0.0}; + if (unlikely(!std::isfinite(value_of_rec(out.first)))) { + return std::optional{std::nullopt}; + } if constexpr (is_autodiff_v) { if constexpr (!partials_fvar) { - out.dlogQ_dalpha + out.second = grad_reg_inc_gamma(alpha, beta_y, tgamma(alpha), digamma(alpha)) - / exp(out.log_Q); + / exp(out.first); } else { - T alpha_unit = alpha; + auto alpha_unit = alpha; alpha_unit.d_ = 1; - T beta_y_unit = beta_y; + auto beta_y_unit = beta_y; beta_y_unit.d_ = 0; - T log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); - out.dlogQ_dalpha = log_Q_fvar.d_; + auto log_Q_fvar = internal::log_q_gamma_cf(alpha_unit, beta_y_unit); + out.second = log_Q_fvar.d_; } } - } - - if (std::isfinite(value_of_rec(out.log_Q))) { - return std::optional>>{out}; - } else { - return std::optional>>{std::nullopt}; + return std::optional{out}; } } /** * Computes log q and d(log q) / d(alpha) using log1m. */ -template -static inline Q_eval eval_q_log1m(const T& alpha, const T& beta_y) { - Q_eval out; - out.log_Q = log1m(gamma_p(alpha, beta_y)); - if (!std::isfinite(value_of_rec(out.log_Q))) { - return out; +template +inline auto eval_q_log1m(const T1& alpha, const T2& beta_y) { + using scalar_t = return_type_t; + using ret_t = std::pair; + ret_t out{log1m(gamma_p(alpha, beta_y)), 0.0}; + if (unlikely(!std::isfinite(value_of_rec(out.first)))) { + return std::optional{std::nullopt}; } if constexpr (is_autodiff_v) { if constexpr (partials_fvar) { - T alpha_unit = alpha; + auto alpha_unit = alpha; alpha_unit.d_ = 1; - T beta_unit = beta_y; + auto beta_unit = beta_y; beta_unit.d_ = 0; - T log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); - out.dlogQ_dalpha = log_Q_fvar.d_; + auto log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); + out.second = log_Q_fvar.d_; } else { - out.dlogQ_dalpha - = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.log_Q); + out.second + = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.first); } } - return out; + return std::optional{out}; } } // namespace internal @@ -134,58 +132,56 @@ inline return_type_t gamma_lccdf( for (size_t n = 0; n < N; n++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero - const T_partials_return y_dbl = y_vec.val(n); - if (y_dbl == 0.0) { + const auto y_val = y_vec.val(n); + if (y_val == 0.0) { continue; } - if (y_dbl == INFTY) { + if (y_val == INFTY) { return ops_partials.build(negative_infinity()); } - const T_partials_return alpha_dbl = alpha_vec.val(n); - const T_partials_return beta_dbl = beta_vec.val(n); + const auto alpha_val = alpha_vec.val(n); + const auto beta_val = beta_vec.val(n); - const T_partials_return beta_y = beta_dbl * y_dbl; + const auto beta_y = beta_val * y_val; if (beta_y == INFTY) { return ops_partials.build(negative_infinity()); } - - const bool use_continued_fraction = beta_y > alpha_dbl + 1.0; - std::optional> result; - if (use_continued_fraction) { - result = internal::eval_q_cf(alpha_dbl, beta_y); + std::optional> result; + if (beta_y > alpha_val + 1.0) { + result = internal::eval_q_cf(alpha_val, beta_y); } else { - result = internal::eval_q_log1m(alpha_dbl, beta_y); + result = internal::eval_q_log1m(alpha_val, beta_y); if (!result && beta_y > 0.0) { // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf(alpha_dbl, beta_y); + result = internal::eval_q_cf(alpha_val, beta_y); } } - if (!result) { + if (unlikely(!result)) { return ops_partials.build(negative_infinity()); } - P += result->log_Q; + P += result->first; if constexpr (is_autodiff_v || is_autodiff_v) { - const T_partials_return log_y = log(y_dbl); - const T_partials_return alpha_minus_one = fma(alpha_dbl, log_y, -log_y); + const auto log_y = log(y_val); + const auto alpha_minus_one = fma(alpha_val, log_y, -log_y); - const T_partials_return log_pdf = alpha_dbl * log(beta_dbl) - - lgamma(alpha_dbl) + alpha_minus_one + const auto log_pdf = alpha_val * log(beta_val) + - lgamma(alpha_val) + alpha_minus_one - beta_y; - const T_partials_return hazard = exp(log_pdf - result->log_Q); // f/Q + const auto hazard = exp(log_pdf - result->first); // f/Q if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] -= hazard; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] -= (y_dbl / beta_dbl) * hazard; + partials<2>(ops_partials)[n] -= (y_val / beta_val) * hazard; } } if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] += result->dlogQ_dalpha; + partials<1>(ops_partials)[n] += result->second; } } return ops_partials.build(P); From 6e1dc0af543749daaa7a84e813bf3ad51fd55b09 Mon Sep 17 00:00:00 2001 From: Daniel Lee Date: Wed, 11 Mar 2026 22:20:17 -0400 Subject: [PATCH 06/11] explicit types --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 32 +++++++++++------------ stan/math/prim/prob/gamma_lccdf.hpp | 25 ++++++++++-------- 2 files changed, 30 insertions(+), 27 deletions(-) diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index c9c110ae56c..32f040c3224 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -39,23 +39,23 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, double precision = 1.49012e-08, int max_steps = 250) { using T_return = return_type_t; - const auto log_prefactor = a * log(z) - z - lgamma(a); + const T_return log_prefactor = a * log(z) - z - lgamma(a); - auto b_init = z + 1.0 - a; - auto C = (fabs(value_of(b_init)) >= EPSILON) ? b_init : std::decay_t(EPSILON); - auto D = 0.0; - auto f = C; + T_return b_init = z + 1.0 - a; + T_return C = (fabs(value_of(b_init)) >= EPSILON) ? b_init : std::decay_t(EPSILON); + T_return D = 0.0; + T_return f = C; for (int i = 1; i <= max_steps; ++i) { T_a an = -i * (i - a); - const auto b = b_init + 2.0 * i; + const T_return b = b_init + 2.0 * i; D = b + an * D; D = (fabs(value_of(D)) >= EPSILON) ? D : std::decay_t(EPSILON); C = b + an / C; C = (fabs(value_of(C)) >= EPSILON) ? C : std::decay_t(EPSILON); D = inv(D); - auto delta = C * D; + const T_return delta = C * D; f *= delta; - const auto delta_m1 = fabs(value_of(delta) - 1.0); + const double delta_m1 = fabs(value_of(delta) - 1.0); if (delta_m1 < precision) { break; } @@ -83,30 +83,30 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, * @return structure containing log(Q(a,z)) and d/da log(Q(a,z)) */ template -inline auto log_gamma_q_dgamma( +inline std::pair, return_type_t> log_gamma_q_dgamma( const T_a& a, const T_z& z, double precision = 1.49012e-08, int max_steps = 250) { using T_return = return_type_t; - const auto a_val = value_of(a); - const auto z_val = value_of(z); + const double a_val = value_of(a); + const double z_val = value_of(z); // For z > a + 1, use continued fraction for better numerical stability if (z_val > a_val + 1.0) { std::pair result{internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0}; // For gradient, use: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da - const auto Q_val = exp(result.first); - const auto dQ_da + const T_return Q_val = exp(result.first); + const double dQ_da = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val)); result.second = dQ_da / Q_val; return result; } else { // For z <= a + 1, use log1m(P(a,z)) for better numerical accuracy - const auto P_val = gamma_p(a_val, z_val); + const double P_val = gamma_p(a_val, z_val); std::pair result{log1m(P_val), 0.0}; // Gradient: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da - const auto Q_val = exp(result.first); + const T_return Q_val = exp(result.first); if (Q_val > 0) { - const auto dQ_da + const double dQ_da = grad_reg_inc_gamma(a_val, z_val, tgamma(a_val), digamma(a_val)); result.second = dQ_da / Q_val; } else { diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index 97df0e29fbf..67c6295ebeb 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -23,6 +23,7 @@ #include #include #include + namespace stan { namespace math { namespace internal { @@ -31,11 +32,12 @@ namespace internal { * Computes log q and d(log q) / d(alpha) using continued fraction. */ template -inline auto eval_q_cf(const T1& alpha, const T2& beta_y) { +inline std::optional, return_type_t>> +eval_q_cf(const T1& alpha, const T2& beta_y) { using scalar_t = return_type_t; using ret_t = std::pair; if constexpr (!any_fvar && is_autodiff_v) { - auto log_q_result + std::pair log_q_result = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); if (likely(std::isfinite(value_of_rec(log_q_result.first)))) { return std::optional{log_q_result}; @@ -69,7 +71,8 @@ inline auto eval_q_cf(const T1& alpha, const T2& beta_y) { * Computes log q and d(log q) / d(alpha) using log1m. */ template -inline auto eval_q_log1m(const T1& alpha, const T2& beta_y) { +inline std::optional, return_type_t>> +eval_q_log1m(const T1& alpha, const T2& beta_y) { using scalar_t = return_type_t; using ret_t = std::pair; ret_t out{log1m(gamma_p(alpha, beta_y)), 0.0}; @@ -132,7 +135,7 @@ inline return_type_t gamma_lccdf( for (size_t n = 0; n < N; n++) { // Explicit results for extreme values // The gradients are technically ill-defined, but treated as zero - const auto y_val = y_vec.val(n); + const T_partials_return y_val = y_vec.val(n); if (y_val == 0.0) { continue; } @@ -140,10 +143,10 @@ inline return_type_t gamma_lccdf( return ops_partials.build(negative_infinity()); } - const auto alpha_val = alpha_vec.val(n); - const auto beta_val = beta_vec.val(n); + const T_partials_return alpha_val = alpha_vec.val(n); + const T_partials_return beta_val = beta_vec.val(n); - const auto beta_y = beta_val * y_val; + const T_partials_return beta_y = beta_val * y_val; if (beta_y == INFTY) { return ops_partials.build(negative_infinity()); } @@ -164,14 +167,14 @@ inline return_type_t gamma_lccdf( P += result->first; if constexpr (is_autodiff_v || is_autodiff_v) { - const auto log_y = log(y_val); - const auto alpha_minus_one = fma(alpha_val, log_y, -log_y); + const T_partials_return log_y = log(y_val); + const T_partials_return alpha_minus_one = fma(alpha_val, log_y, -log_y); - const auto log_pdf = alpha_val * log(beta_val) + const T_partials_return log_pdf = alpha_val * log(beta_val) - lgamma(alpha_val) + alpha_minus_one - beta_y; - const auto hazard = exp(log_pdf - result->first); // f/Q + const T_partials_return hazard = exp(log_pdf - result->first); // f/Q if constexpr (is_autodiff_v) { partials<0>(ops_partials)[n] -= hazard; From 77240117f10bc7e920259e4de935108d74c6b866 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Mon, 16 Mar 2026 17:13:20 -0400 Subject: [PATCH 07/11] update to fix value_of and value_of_rec in gamma functions --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 27 +++++++++++++++++------ stan/math/prim/prob/gamma_lccdf.hpp | 13 ++++++----- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index 32f040c3224..2b68c9e2e40 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -15,6 +15,7 @@ #include #include #include +#include #include namespace stan { @@ -22,6 +23,8 @@ namespace math { namespace internal { +constexpr double LOG_Q_GAMMA_CF_PRECISION = 1.49012e-12; + /** * Compute log(Q(a,z)) using continued fraction expansion for upper incomplete * gamma function. @@ -36,26 +39,33 @@ namespace internal { */ template inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, - double precision = 1.49012e-08, + double precision + = LOG_Q_GAMMA_CF_PRECISION, int max_steps = 250) { using T_return = return_type_t; const T_return log_prefactor = a * log(z) - z - lgamma(a); T_return b_init = z + 1.0 - a; - T_return C = (fabs(value_of(b_init)) >= EPSILON) ? b_init : std::decay_t(EPSILON); + T_return C = (fabs(value_of_rec(b_init)) >= EPSILON) + ? b_init + : std::decay_t(EPSILON); T_return D = 0.0; T_return f = C; for (int i = 1; i <= max_steps; ++i) { T_a an = -i * (i - a); const T_return b = b_init + 2.0 * i; D = b + an * D; - D = (fabs(value_of(D)) >= EPSILON) ? D : std::decay_t(EPSILON); + D = (fabs(value_of_rec(D)) >= EPSILON) + ? D + : std::decay_t(EPSILON); C = b + an / C; - C = (fabs(value_of(C)) >= EPSILON) ? C : std::decay_t(EPSILON); + C = (fabs(value_of_rec(C)) >= EPSILON) + ? C + : std::decay_t(EPSILON); D = inv(D); const T_return delta = C * D; f *= delta; - const double delta_m1 = fabs(value_of(delta) - 1.0); + const double delta_m1 = fabs(value_of_rec(delta) - 1.0); if (delta_m1 < precision) { break; } @@ -84,13 +94,16 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, */ template inline std::pair, return_type_t> log_gamma_q_dgamma( - const T_a& a, const T_z& z, double precision = 1.49012e-08, int max_steps = 250) { + const T_a& a, const T_z& z, + double precision = internal::LOG_Q_GAMMA_CF_PRECISION, + int max_steps = 250) { using T_return = return_type_t; const double a_val = value_of(a); const double z_val = value_of(z); // For z > a + 1, use continued fraction for better numerical stability if (z_val > a_val + 1.0) { - std::pair result{internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0}; + std::pair result{ + internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0}; // For gradient, use: d/da log(Q) = (1/Q) * dQ/da // grad_reg_inc_gamma computes dQ/da const T_return Q_val = exp(result.first); diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index 67c6295ebeb..ffd3ad318b9 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -1,6 +1,7 @@ #ifndef STAN_MATH_PRIM_PROB_GAMMA_LCCDF_HPP #define STAN_MATH_PRIM_PROB_GAMMA_LCCDF_HPP +#include #include #include #include @@ -38,8 +39,8 @@ eval_q_cf(const T1& alpha, const T2& beta_y) { using ret_t = std::pair; if constexpr (!any_fvar && is_autodiff_v) { std::pair log_q_result - = log_gamma_q_dgamma(value_of_rec(alpha), value_of_rec(beta_y)); - if (likely(std::isfinite(value_of_rec(log_q_result.first)))) { + = log_gamma_q_dgamma(value_of(alpha), value_of(beta_y)); + if (likely(std::isfinite(log_q_result.first))) { return std::optional{log_q_result}; } else { return std::optional{std::nullopt}; @@ -127,10 +128,10 @@ inline return_type_t gamma_lccdf( scalar_seq_view beta_vec(beta_ref); const size_t N = max_size(y, alpha, beta); - constexpr bool any_fvar = is_fvar>::value - || is_fvar>::value - || is_fvar>::value; - constexpr bool partials_fvar = is_fvar::value; + constexpr bool any_fvar = is_fvar_v> + || is_fvar_v> + || is_fvar_v>; + constexpr bool partials_fvar = is_fvar_v; for (size_t n = 0; n < N; n++) { // Explicit results for extreme values From 8e29962778f26e338e0a2691742fe9de3a3fba11 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 16 Mar 2026 20:49:27 -0400 Subject: [PATCH 08/11] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/log_gamma_q_dgamma.hpp | 18 ++++++++---------- stan/math/prim/prob/gamma_lccdf.hpp | 21 ++++++++++++--------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/fun/log_gamma_q_dgamma.hpp b/stan/math/prim/fun/log_gamma_q_dgamma.hpp index 2b68c9e2e40..742b23784ba 100644 --- a/stan/math/prim/fun/log_gamma_q_dgamma.hpp +++ b/stan/math/prim/fun/log_gamma_q_dgamma.hpp @@ -55,13 +55,11 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, T_a an = -i * (i - a); const T_return b = b_init + 2.0 * i; D = b + an * D; - D = (fabs(value_of_rec(D)) >= EPSILON) - ? D - : std::decay_t(EPSILON); + D = (fabs(value_of_rec(D)) >= EPSILON) ? D + : std::decay_t(EPSILON); C = b + an / C; - C = (fabs(value_of_rec(C)) >= EPSILON) - ? C - : std::decay_t(EPSILON); + C = (fabs(value_of_rec(C)) >= EPSILON) ? C + : std::decay_t(EPSILON); D = inv(D); const T_return delta = C * D; f *= delta; @@ -93,10 +91,10 @@ inline return_type_t log_q_gamma_cf(const T_a& a, const T_z& z, * @return structure containing log(Q(a,z)) and d/da log(Q(a,z)) */ template -inline std::pair, return_type_t> log_gamma_q_dgamma( - const T_a& a, const T_z& z, - double precision = internal::LOG_Q_GAMMA_CF_PRECISION, - int max_steps = 250) { +inline std::pair, return_type_t> +log_gamma_q_dgamma(const T_a& a, const T_z& z, + double precision = internal::LOG_Q_GAMMA_CF_PRECISION, + int max_steps = 250) { using T_return = return_type_t; const double a_val = value_of(a); const double z_val = value_of(z); diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index ffd3ad318b9..d625c61908f 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -32,7 +32,8 @@ namespace internal { /** * Computes log q and d(log q) / d(alpha) using continued fraction. */ -template +template inline std::optional, return_type_t>> eval_q_cf(const T1& alpha, const T2& beta_y) { using scalar_t = return_type_t; @@ -89,8 +90,7 @@ eval_q_log1m(const T1& alpha, const T2& beta_y) { auto log_Q_fvar = log1m(gamma_p(alpha_unit, beta_unit)); out.second = log_Q_fvar.d_; } else { - out.second - = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.first); + out.second = -grad_reg_lower_inc_gamma(alpha, beta_y) / exp(out.first); } } return std::optional{out}; @@ -128,9 +128,9 @@ inline return_type_t gamma_lccdf( scalar_seq_view beta_vec(beta_ref); const size_t N = max_size(y, alpha, beta); - constexpr bool any_fvar = is_fvar_v> - || is_fvar_v> - || is_fvar_v>; + constexpr bool any_fvar + = is_fvar_v> || is_fvar_v> || is_fvar_v>; constexpr bool partials_fvar = is_fvar_v; for (size_t n = 0; n < N; n++) { @@ -153,12 +153,15 @@ inline return_type_t gamma_lccdf( } std::optional> result; if (beta_y > alpha_val + 1.0) { - result = internal::eval_q_cf(alpha_val, beta_y); + result = internal::eval_q_cf(alpha_val, + beta_y); } else { - result = internal::eval_q_log1m(alpha_val, beta_y); + result + = internal::eval_q_log1m(alpha_val, beta_y); if (!result && beta_y > 0.0) { // Fallback to continued fraction if log1m fails - result = internal::eval_q_cf(alpha_val, beta_y); + result = internal::eval_q_cf( + alpha_val, beta_y); } } if (unlikely(!result)) { From 352ef47853d6d9bd4e5e06eaa483cb3df444b7e8 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 17 Mar 2026 11:57:54 -0400 Subject: [PATCH 09/11] breakdown large compile time conditional statement into smaller lines --- stan/math/prim/prob/gamma_lccdf.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index d625c61908f..5545f2b7a7c 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -128,9 +128,10 @@ inline return_type_t gamma_lccdf( scalar_seq_view beta_vec(beta_ref); const size_t N = max_size(y, alpha, beta); - constexpr bool any_fvar - = is_fvar_v> || is_fvar_v> || is_fvar_v>; + constexpr bool is_y_fvar = is_fvar_v>; + constexpr bool is_shape_fvar = is_fvar_v>; + constexpr bool is_beta_fvar = is_fvar_v>; + constexpr bool any_fvar = is_y_fvar || is_shape_fvar || is_beta_fvar; constexpr bool partials_fvar = is_fvar_v; for (size_t n = 0; n < N; n++) { From f16a0c0636795bec5520935bc1c192c247386f4e Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 17 Mar 2026 14:32:59 -0400 Subject: [PATCH 10/11] remove fwd in prim include --- stan/math/prim/prob/gamma_lccdf.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp index 5545f2b7a7c..c5538fa5966 100644 --- a/stan/math/prim/prob/gamma_lccdf.hpp +++ b/stan/math/prim/prob/gamma_lccdf.hpp @@ -1,7 +1,6 @@ #ifndef STAN_MATH_PRIM_PROB_GAMMA_LCCDF_HPP #define STAN_MATH_PRIM_PROB_GAMMA_LCCDF_HPP -#include #include #include #include From b0feb956656ce2c1fe9d04f12c178fd582f8d017 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 17 Mar 2026 14:52:45 -0400 Subject: [PATCH 11/11] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/meta/is_fvar.hpp | 1 - test/unit/math/rev/prob/gamma_lccdf_test.cpp | 4 +++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/fwd/meta/is_fvar.hpp b/stan/math/fwd/meta/is_fvar.hpp index 8102aeb0fb3..e208d08bc71 100644 --- a/stan/math/fwd/meta/is_fvar.hpp +++ b/stan/math/fwd/meta/is_fvar.hpp @@ -21,6 +21,5 @@ struct is_fvar>::value>> : std::true_type {}; - } // namespace stan #endif diff --git a/test/unit/math/rev/prob/gamma_lccdf_test.cpp b/test/unit/math/rev/prob/gamma_lccdf_test.cpp index 72bb85f2342..d069a365578 100644 --- a/test/unit/math/rev/prob/gamma_lccdf_test.cpp +++ b/test/unit/math/rev/prob/gamma_lccdf_test.cpp @@ -231,7 +231,9 @@ TEST_F(AgradRev, ProbDistributionsGamma_lccdf_extreme_values_small) { } } -TEST_F(AgradRev, ProbDistributionsGammalccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) { +TEST_F( + AgradRev, + ProbDistributionsGammalccdf_alpha_gt_30_small_y_old_code_rounds_to_zero) { using stan::math::gamma_lccdf; using stan::math::gamma_p; using stan::math::gamma_q;