update to fix value_of and value_of_rec in gamma functions

Author: SteveBronder

stan/math/prim/fun/log_gamma_q_dgamma.hpp

@@ -15,13 +15,16 @@
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
 #include <cmath>
 
 namespace stan {
 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 <typename T_a, typename T_z>
 inline return_type_t<T_a, T_z> 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<T_a, T_z>;
   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<decltype(b_init)>(EPSILON);
+  T_return C = (fabs(value_of_rec(b_init)) >= EPSILON)
+                   ? b_init
+                   : std::decay_t<decltype(b_init)>(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<decltype(D)>(EPSILON);
+    D = (fabs(value_of_rec(D)) >= EPSILON)
+            ? D
+            : std::decay_t<decltype(D)>(EPSILON);
     C = b + an / C;
-    C = (fabs(value_of(C)) >= EPSILON) ? C : std::decay_t<decltype(C)>(EPSILON);
+    C = (fabs(value_of_rec(C)) >= EPSILON)
+            ? C
+            : std::decay_t<decltype(C)>(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<T_a, T_z> log_q_gamma_cf(const T_a& a, const T_z& z,
  */
 template <typename T_a, typename T_z>
 inline std::pair<return_type_t<T_a, T_z>, return_type_t<T_a, T_z>> 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<T_a, T_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<T_return, T_return> result{internal::log_q_gamma_cf(a_val, z_val, precision, max_steps), 0.0};
+    std::pair<T_return, T_return> 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);

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 <stan/math/fwd/meta/is_fvar.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
@@ -38,8 +39,8 @@ eval_q_cf(const T1& alpha, const T2& beta_y) {
   using ret_t = std::pair<scalar_t, scalar_t>;
   if constexpr (!any_fvar && is_autodiff_v<T_shape>) {
     std::pair<double, double> 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<ret_t>{std::nullopt};
@@ -127,10 +128,10 @@ inline return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(
   scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
   const size_t N = max_size(y, alpha, beta);
 
-  constexpr bool any_fvar = is_fvar<scalar_type_t<T_y>>::value
-                            || is_fvar<scalar_type_t<T_shape>>::value
-                            || is_fvar<scalar_type_t<T_inv_scale>>::value;
-  constexpr bool partials_fvar = is_fvar<T_partials_return>::value;
+  constexpr bool any_fvar = is_fvar_v<scalar_type_t<T_y>>
+                            || is_fvar_v<scalar_type_t<T_shape>>
+                            || is_fvar_v<scalar_type_t<T_inv_scale>>;
+  constexpr bool partials_fvar = is_fvar_v<T_partials_return>;
 
   for (size_t n = 0; n < N; n++) {
     // Explicit results for extreme values