explicit types

Author: syclik

stan/math/prim/fun/log_gamma_q_dgamma.hpp

@@ -39,23 +39,23 @@ 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,
                                               int max_steps = 250) {
   using T_return = return_type_t<T_a, T_z>;
-  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<decltype(b_init)>(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<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 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<decltype(D)>(EPSILON);
     C = b + an / C;
     C = (fabs(value_of(C)) >= EPSILON) ? C : std::decay_t<decltype(C)>(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<T_a, T_z> 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 <typename T_a, typename T_z>
-inline auto log_gamma_q_dgamma(
+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) {
   using T_return = return_type_t<T_a, T_z>;
-  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<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 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<T_return, T_return> 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 {

stan/math/prim/prob/gamma_lccdf.hpp

@@ -23,6 +23,7 @@
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 #include <optional>
+
 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 <bool any_fvar, bool partials_fvar, typename T_shape, typename T1, typename T2>
-inline auto eval_q_cf(const T1& alpha, const T2& beta_y) {
+inline std::optional<std::pair<return_type_t<T1, T2>, return_type_t<T1, T2>>>
+eval_q_cf(const T1& alpha, const T2& beta_y) {
   using scalar_t = return_type_t<T1, T2>;
   using ret_t = std::pair<scalar_t, scalar_t>;
   if constexpr (!any_fvar && is_autodiff_v<T_shape>) {
-    auto log_q_result
+    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)))) {
       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 <bool partials_fvar, typename T_shape, typename T1, typename T2>
-inline auto eval_q_log1m(const T1& alpha, const T2& beta_y) {
+inline std::optional<std::pair<return_type_t<T1, T2>, return_type_t<T1, T2>>>
+eval_q_log1m(const T1& alpha, const T2& beta_y) {
   using scalar_t = return_type_t<T1, T2>;
   using ret_t = std::pair<scalar_t, scalar_t>;
   ret_t out{log1m(gamma_p(alpha, beta_y)), 0.0};
@@ -132,18 +135,18 @@ inline return_type_t<T_y, T_shape, T_inv_scale> 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;
     }
     if (y_val == INFTY) {
       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<T_y, T_shape, T_inv_scale> gamma_lccdf(
     P += result->first;
 
     if constexpr (is_autodiff_v<T_y> || is_autodiff_v<T_inv_scale>) {
-      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<T_y>) {
         partials<0>(ops_partials)[n] -= hazard;