From 76ba3037aa2113c5555e208a00601ec62a3d29e4 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 19 Jun 2025 13:50:04 +0800
Subject: [PATCH 1/9] Initial implementation

---
 stan/math/fwd/prob.hpp               |  1 +
 stan/math/fwd/prob/student_t_qf.hpp  | 82 ++++++++++++++++++++++++++++
 stan/math/prim/prob.hpp              |  1 +
 stan/math/prim/prob/student_t_qf.hpp | 37 +++++++++++++
 stan/math/rev/prob.hpp               |  1 +
 stan/math/rev/prob/student_t_qf.hpp  | 80 +++++++++++++++++++++++++++
 6 files changed, 202 insertions(+)
 create mode 100644 stan/math/fwd/prob/student_t_qf.hpp
 create mode 100644 stan/math/prim/prob/student_t_qf.hpp
 create mode 100644 stan/math/rev/prob/student_t_qf.hpp

diff --git a/stan/math/fwd/prob.hpp b/stan/math/fwd/prob.hpp
index 5d280c06525..c390e191db5 100644
--- a/stan/math/fwd/prob.hpp
+++ b/stan/math/fwd/prob.hpp
@@ -5,5 +5,6 @@
 #include <stan/math/fwd/fun/Eigen_NumTraits.hpp>
 
 #include <stan/math/fwd/prob/std_normal_log_qf.hpp>
+#include <stan/math/fwd/prob/student_t_qf.hpp>
 
 #endif
diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
new file mode 100644
index 00000000000..f2ba00f27f3
--- /dev/null
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -0,0 +1,82 @@
+#ifndef STAN_MATH_FWD_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_FWD_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/fun/sqrt.hpp>
+#include <stan/math/fwd/fun/inv_inc_beta.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T_p, typename T_mu, typename T_sigma, typename T_nu,
+          typename FvarT = return_type_t<T_p, T_mu, T_sigma, T_nu>,
+          require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
+          require_fvar_t<FvarT>* = nullptr>
+FvarT student_t_qf(const T_p& p, const T_mu& mu,
+                    const T_sigma& sigma, const T_nu& nu) {
+  static constexpr const char* function = "student_t_qf";
+  using T_partials = partials_type_t<FvarT>;
+
+  const auto& p_val = value_of(p);
+  const auto& mu_val = value_of(mu);
+  const auto& sigma_val = value_of(sigma);
+  const auto& nu_val = value_of(nu);
+
+  check_nonnegative(function, "Degrees of freedom parameter", nu_val);
+  check_positive(function, "Scale parameter", sigma_val);
+  check_bounded(function, "Probability parameter", p_val, 0.0, 1.0);
+
+  if (p_val == 0.0) {
+    return {NEGATIVE_INFTY, 0.0};
+  } else if (p_val == 1.0) {
+    return {INFTY, 0.0};
+  } else if (p_val == 0.5) {
+    return {mu_val, 0.0};
+  }
+
+
+  const T_partials p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+  const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+
+  T_partials ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
+  FvarT rtn = mu_val + p_sign * sigma_val * sqrt(nu_val)
+                * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+
+  if constexpr (!is_constant<T_p>::value) {
+    rtn.d_ += p.d_ * exp(log(p_val) - student_t_lpdf(rtn, mu_val, sigma_val, nu_val));
+  }
+
+  if constexpr (!is_constant<T_mu>::value) {
+    if (p_val > 0.0 && p_val < 1.0) {
+      mu.d_ += rtn.d_;
+    }
+  }
+
+  if constexpr (!is_constant<T_sigma>::value) {
+    if (p_val > 0.0 && p_val < 1.0) {
+      sigma.d_ += rtn.d_;
+    }
+  }
+
+  if constexpr (!is_constant<T_nu>::value) {
+    const T_partials half_nu = nu_val / 2.0;
+    const T_partials hyper_arg = hypergeometric_3F2({0.5, half_nu, half_nu}, {1 + half_nu, 1 + half_nu}, ibeta_arg);
+    const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
+    const T_partials digamma_a1 = digamma(half_nu);
+    const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
+
+    const T_partials arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
+    const T_partials arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg) * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+
+    const T_partials num1 = p_sign * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
+    const T_partials den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
+    nu.d_ += rtn.d_ * p_sign * num1 / den1;
+  }
+
+  return rtn;
+
+}
+}
+}
+
+#endif
diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index e04c7e87172..4e3207f0c68 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -287,6 +287,7 @@
 #include <stan/math/prim/prob/student_t_lccdf.hpp>
 #include <stan/math/prim/prob/student_t_lcdf.hpp>
 #include <stan/math/prim/prob/student_t_lpdf.hpp>
+#include <stan/math/prim/prob/student_t_qf.hpp>
 #include <stan/math/prim/prob/student_t_rng.hpp>
 #include <stan/math/prim/prob/uniform_ccdf_log.hpp>
 #include <stan/math/prim/prob/uniform_cdf.hpp>
diff --git a/stan/math/prim/prob/student_t_qf.hpp b/stan/math/prim/prob/student_t_qf.hpp
new file mode 100644
index 00000000000..4e6078ecf10
--- /dev/null
+++ b/stan/math/prim/prob/student_t_qf.hpp
@@ -0,0 +1,37 @@
+#ifndef STAN_MATH_PRIM_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_PRIM_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
+#include <stan/math/prim/fun/inv_inc_beta.hpp>
+
+namespace stan {
+namespace math {
+
+double student_t_qf(const double p, const double mu,
+                    const double sigma, const double nu) {
+  static constexpr const char* function = "student_t_qf";
+  check_nonnegative(function, "Degrees of freedom parameter", nu);
+  check_positive(function, "Scale parameter", sigma);
+  check_bounded(function, "Probability parameter", p, 0.0, 1.0);
+
+  if (p == 0.0) {
+    return NEGATIVE_INFTY;
+  } else if (p == 1.0) {
+    return INFTY;
+  } else if (p == 0.5) {
+    return mu;
+  }
+
+  const double p_val_flip = p < 0.5 ? p : 1.0 - p;
+  const double p_sign = p < 0.5 ? -1.0 : 1.0;
+  const auto& ibeta_arg = inv_inc_beta(0.5 * nu, 0.5, 2 * p_val_flip);
+
+  return mu + p_sign * sigma * math::sqrt(nu)
+          * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+}
+}
+}
+
+#endif
diff --git a/stan/math/rev/prob.hpp b/stan/math/rev/prob.hpp
index 98a6fee613f..f1f2427fdd3 100644
--- a/stan/math/rev/prob.hpp
+++ b/stan/math/rev/prob.hpp
@@ -4,5 +4,6 @@
 #include <stan/math/prim/fun/Eigen.hpp>
 
 #include <stan/math/rev/prob/std_normal_log_qf.hpp>
+#include <stan/math/rev/prob/student_t_qf.hpp>
 
 #endif
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
new file mode 100644
index 00000000000..6fab70b4fe0
--- /dev/null
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -0,0 +1,80 @@
+#ifndef STAN_MATH_REV_PROB_STUDENT_T_QF_HPP
+#define STAN_MATH_REV_PROB_STUDENT_T_QF_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/sqrt.hpp>
+#include <stan/math/rev/fun/inv_inc_beta.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename T_p, typename T_mu, typename T_sigma, typename T_nu,
+          require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
+          require_any_var_t<T_p, T_mu, T_sigma, T_nu>* = nullptr>
+var student_t_qf(const T_p& p, const T_mu& mu,
+                    const T_sigma& sigma, const T_nu& nu) {
+  static constexpr const char* function = "student_t_qf";
+  const double p_val = value_of(p);
+  const double mu_val = value_of(mu);
+  const double sigma_val = value_of(sigma);
+  const double nu_val = value_of(nu);
+
+  check_nonnegative(function, "Degrees of freedom parameter", nu_val);
+  check_positive(function, "Scale parameter", sigma_val);
+  check_bounded(function, "Probability parameter", p_val, 0.0, 1.0);
+
+  double ret_val(0);
+  double ibeta_arg(0);
+
+  if (p == 0.0) {
+    ret_val = NEGATIVE_INFTY
+  } else if (p == 1.0) {
+    ret_val = INFTY;
+  } else if (p == 0.5) {
+    ret_val = mu_val;
+  } else {
+    const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+    const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+    ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
+
+    ret_val = mu_val + p_sign * sigma_val * math::sqrt(nu_val)
+                  * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+  }
+
+  return make_callback_var(ret_val, [p, mu, sigma, nu, ibeta_arg, nu_val, p_sign](auto& vi) mutable {
+    if constexpr (!is_constant<T_p>::value) {
+      if (p.val() > 0.0 && p.val() < 1.0) {
+        p.adj() += vi.adj() * exp(log(p_val) - student_t_lpdf(ret_val, mu.val(), sigma.val(), nu_val));
+      }
+    }
+    if constexpr (!is_constant<T_mu>::value) {
+      if (p_val > 0.0 && p_val < 1.0) {
+        mu.adj() += vi.adj();
+      }
+    }
+    if constexpr (!is_constant<T_sigma>::value) {
+      if (p_val > 0.0 && p_val < 1.0) {
+        sigma.adj() += vi.adj();
+      }
+    }
+    if constexpr (!is_constant<T_nu>::value) {
+      const double half_nu = nu_val / 2.0;
+      const double hyper_arg = hypergeometric_3F2({0.5, half_nu, half_nu}, {1 + half_nu, 1 + half_nu}, ibeta_arg);
+      const double hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
+      const double digamma_a1 = digamma(half_nu);
+      const double digamma_a2 = digamma((1 + nu_val) / 2);
+
+      const double arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
+      const double arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg) * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+
+      const double num1 = s * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
+      const double den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
+      nu.adj() += vi.adj() * p_sign * num1 / den1;
+    }
+  });
+
+}
+}
+}
+
+#endif

From f8c67018d4d85a36c03e278ee4dd5c7c666e2d15 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 19 Jun 2025 13:53:13 +0800
Subject: [PATCH 2/9] Compile fixs

---
 stan/math/rev/prob/student_t_qf.hpp | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index 6fab70b4fe0..39e31565dd5 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -27,7 +27,7 @@ var student_t_qf(const T_p& p, const T_mu& mu,
   double ibeta_arg(0);
 
   if (p == 0.0) {
-    ret_val = NEGATIVE_INFTY
+    ret_val = NEGATIVE_INFTY;
   } else if (p == 1.0) {
     ret_val = INFTY;
   } else if (p == 0.5) {
@@ -41,7 +41,8 @@ var student_t_qf(const T_p& p, const T_mu& mu,
                   * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
-  return make_callback_var(ret_val, [p, mu, sigma, nu, ibeta_arg, nu_val, p_sign](auto& vi) mutable {
+  return make_callback_var(ret_val, [p, mu, sigma, nu, ibeta_arg, nu_val, ret_val](auto& vi) mutable {
+    const double p_val = value_of(p);
     if constexpr (!is_constant<T_p>::value) {
       if (p.val() > 0.0 && p.val() < 1.0) {
         p.adj() += vi.adj() * exp(log(p_val) - student_t_lpdf(ret_val, mu.val(), sigma.val(), nu_val));
@@ -58,6 +59,8 @@ var student_t_qf(const T_p& p, const T_mu& mu,
       }
     }
     if constexpr (!is_constant<T_nu>::value) {
+      const double sigma_val = value_of(sigma);
+      const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
       const double half_nu = nu_val / 2.0;
       const double hyper_arg = hypergeometric_3F2({0.5, half_nu, half_nu}, {1 + half_nu, 1 + half_nu}, ibeta_arg);
       const double hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
@@ -67,7 +70,7 @@ var student_t_qf(const T_p& p, const T_mu& mu,
       const double arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
       const double arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg) * (log(ibeta_arg) - digamma_a1 + digamma_a2);
 
-      const double num1 = s * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
+      const double num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
       const double den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
       nu.adj() += vi.adj() * p_sign * num1 / den1;
     }

From d8ba50cc8185a34922e46989dcc8e8cca10b4893 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 19 Jun 2025 15:48:38 +0800
Subject: [PATCH 3/9] Initial tests

---
 stan/math/fwd/prob/student_t_qf.hpp           | 47 ++++++++++---------
 stan/math/prim/prob/student_t_qf.hpp          |  4 +-
 stan/math/rev/prob/student_t_qf.hpp           | 47 ++++++++++++-------
 test/unit/math/mix/prob/student_t_qf_test.cpp | 22 +++++++++
 .../unit/math/prim/prob/student_t_qf_test.cpp | 28 +++++++++++
 5 files changed, 106 insertions(+), 42 deletions(-)
 create mode 100644 test/unit/math/mix/prob/student_t_qf_test.cpp
 create mode 100644 test/unit/math/prim/prob/student_t_qf_test.cpp

diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
index f2ba00f27f3..5ed093dd515 100644
--- a/stan/math/fwd/prob/student_t_qf.hpp
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -4,23 +4,24 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun/sqrt.hpp>
 #include <stan/math/fwd/fun/inv_inc_beta.hpp>
+#include <stan/math/prim/prob/student_t_lpdf.hpp>
 
 namespace stan {
 namespace math {
 
-template <typename T_p, typename T_mu, typename T_sigma, typename T_nu,
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
           typename FvarT = return_type_t<T_p, T_mu, T_sigma, T_nu>,
           require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
           require_fvar_t<FvarT>* = nullptr>
-FvarT student_t_qf(const T_p& p, const T_mu& mu,
-                    const T_sigma& sigma, const T_nu& nu) {
+FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                    const T_sigma& sigma) {
   static constexpr const char* function = "student_t_qf";
   using T_partials = partials_type_t<FvarT>;
 
   const auto& p_val = value_of(p);
+  const auto& nu_val = value_of(nu);
   const auto& mu_val = value_of(mu);
   const auto& sigma_val = value_of(sigma);
-  const auto& nu_val = value_of(nu);
 
   check_nonnegative(function, "Degrees of freedom parameter", nu_val);
   check_positive(function, "Scale parameter", sigma_val);
@@ -36,31 +37,25 @@ FvarT student_t_qf(const T_p& p, const T_mu& mu,
 
 
   const T_partials p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
-  const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+  const double p_sign = value_of_rec(p_val) < 0.5 ? -1.0 : 1.0;
 
   T_partials ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
-  FvarT rtn = mu_val + p_sign * sigma_val * sqrt(nu_val)
+  T_partials rtn_val = mu_val + p_sign * sigma_val * sqrt(nu_val)
                 * math::sqrt(-1.0 + 1.0 / ibeta_arg);
 
-  if constexpr (!is_constant<T_p>::value) {
-    rtn.d_ += p.d_ * exp(log(p_val) - student_t_lpdf(rtn, mu_val, sigma_val, nu_val));
-  }
-
-  if constexpr (!is_constant<T_mu>::value) {
-    if (p_val > 0.0 && p_val < 1.0) {
-      mu.d_ += rtn.d_;
-    }
-  }
+  FvarT rtn(rtn_val, 0.0);
 
-  if constexpr (!is_constant<T_sigma>::value) {
-    if (p_val > 0.0 && p_val < 1.0) {
-      sigma.d_ += rtn.d_;
-    }
+  if constexpr (!is_constant<T_p>::value) {
+    rtn.d_ += p.d_ * exp(-student_t_lpdf(rtn_val, nu_val, mu_val, sigma_val));
   }
 
   if constexpr (!is_constant<T_nu>::value) {
     const T_partials half_nu = nu_val / 2.0;
-    const T_partials hyper_arg = hypergeometric_3F2({0.5, half_nu, half_nu}, {1 + half_nu, 1 + half_nu}, ibeta_arg);
+    Eigen::Matrix<T_partials, Eigen::Dynamic, 1> hyper_arg_a(3);
+    hyper_arg_a << 0.5, half_nu, half_nu;
+    Eigen::Matrix<T_partials, Eigen::Dynamic, 1> hyper_arg_b(2);
+    hyper_arg_b << 1 + half_nu, 1 + half_nu;
+    const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
     const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
     const T_partials digamma_a1 = digamma(half_nu);
     const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
@@ -68,9 +63,17 @@ FvarT student_t_qf(const T_p& p, const T_mu& mu,
     const T_partials arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
     const T_partials arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg) * (log(ibeta_arg) - digamma_a1 + digamma_a2);
 
-    const T_partials num1 = p_sign * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
+    const T_partials num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
     const T_partials den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
-    nu.d_ += rtn.d_ * p_sign * num1 / den1;
+    rtn.d_ += nu.d_ * p_sign * num1 / den1;
+  }
+
+  if constexpr (!is_constant<T_mu>::value) {
+    rtn.d_ += mu.d_;
+  }
+
+  if constexpr (!is_constant<T_sigma>::value) {
+    rtn.d_ += sigma.d_ * p_sign * sqrt(nu_val) * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
   return rtn;
diff --git a/stan/math/prim/prob/student_t_qf.hpp b/stan/math/prim/prob/student_t_qf.hpp
index 4e6078ecf10..b84ccb37357 100644
--- a/stan/math/prim/prob/student_t_qf.hpp
+++ b/stan/math/prim/prob/student_t_qf.hpp
@@ -9,8 +9,8 @@
 namespace stan {
 namespace math {
 
-double student_t_qf(const double p, const double mu,
-                    const double sigma, const double nu) {
+double student_t_qf(const double p, const double nu, const double mu,
+                    const double sigma) {
   static constexpr const char* function = "student_t_qf";
   check_nonnegative(function, "Degrees of freedom parameter", nu);
   check_positive(function, "Scale parameter", sigma);
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index 39e31565dd5..53aadb1e3c3 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -4,20 +4,21 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/fun/sqrt.hpp>
 #include <stan/math/rev/fun/inv_inc_beta.hpp>
+#include <stan/math/prim/prob/student_t_lpdf.hpp>
 
 namespace stan {
 namespace math {
 
-template <typename T_p, typename T_mu, typename T_sigma, typename T_nu,
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
           require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
           require_any_var_t<T_p, T_mu, T_sigma, T_nu>* = nullptr>
-var student_t_qf(const T_p& p, const T_mu& mu,
-                    const T_sigma& sigma, const T_nu& nu) {
+var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                    const T_sigma& sigma) {
   static constexpr const char* function = "student_t_qf";
   const double p_val = value_of(p);
+  const double nu_val = value_of(nu);
   const double mu_val = value_of(mu);
   const double sigma_val = value_of(sigma);
-  const double nu_val = value_of(nu);
 
   check_nonnegative(function, "Degrees of freedom parameter", nu_val);
   check_positive(function, "Scale parameter", sigma_val);
@@ -41,28 +42,28 @@ var student_t_qf(const T_p& p, const T_mu& mu,
                   * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
-  return make_callback_var(ret_val, [p, mu, sigma, nu, ibeta_arg, nu_val, ret_val](auto& vi) mutable {
+  return make_callback_var(ret_val, [p, mu, sigma, nu, ibeta_arg](auto& vi) mutable {
     const double p_val = value_of(p);
+    const double mu_val = value_of(mu);
+    const double sigma_val = value_of(sigma);
+    const double nu_val = value_of(nu);
+    const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
+    const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+
     if constexpr (!is_constant<T_p>::value) {
       if (p.val() > 0.0 && p.val() < 1.0) {
-        p.adj() += vi.adj() * exp(log(p_val) - student_t_lpdf(ret_val, mu.val(), sigma.val(), nu_val));
-      }
-    }
-    if constexpr (!is_constant<T_mu>::value) {
-      if (p_val > 0.0 && p_val < 1.0) {
-        mu.adj() += vi.adj();
-      }
-    }
-    if constexpr (!is_constant<T_sigma>::value) {
-      if (p_val > 0.0 && p_val < 1.0) {
-        sigma.adj() += vi.adj();
+        p.adj() += vi.adj() * exp(-student_t_lpdf(vi.val(), nu_val, mu_val, sigma_val));
       }
     }
+
     if constexpr (!is_constant<T_nu>::value) {
       const double sigma_val = value_of(sigma);
-      const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
       const double half_nu = nu_val / 2.0;
-      const double hyper_arg = hypergeometric_3F2({0.5, half_nu, half_nu}, {1 + half_nu, 1 + half_nu}, ibeta_arg);
+      Eigen::VectorXd hyper_arg_a(3);
+      hyper_arg_a << 0.5, half_nu, half_nu;
+      Eigen::VectorXd hyper_arg_b(2);
+      hyper_arg_b << 1 + half_nu, 1 + half_nu;
+      const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
       const double hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
       const double digamma_a1 = digamma(half_nu);
       const double digamma_a2 = digamma((1 + nu_val) / 2);
@@ -74,6 +75,16 @@ var student_t_qf(const T_p& p, const T_mu& mu,
       const double den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
       nu.adj() += vi.adj() * p_sign * num1 / den1;
     }
+    if constexpr (!is_constant<T_mu>::value) {
+      if (p_val > 0.0 && p_val < 1.0) {
+        mu.adj() += vi.adj();
+      }
+    }
+    if constexpr (!is_constant<T_sigma>::value) {
+      if (p_val > 0.0 && p_val < 1.0) {
+        sigma.adj() += vi.adj() * p_sign * sqrt(value_of(nu)) * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+      }
+    }
   });
 
 }
diff --git a/test/unit/math/mix/prob/student_t_qf_test.cpp b/test/unit/math/mix/prob/student_t_qf_test.cpp
new file mode 100644
index 00000000000..07c400a729d
--- /dev/null
+++ b/test/unit/math/mix/prob/student_t_qf_test.cpp
@@ -0,0 +1,22 @@
+#include <stan/math/rev.hpp>
+#include <stan/math/prim.hpp>
+#include <test/unit/math/rev/util.hpp>
+#include <test/unit/math/test_ad.hpp>
+
+TEST_F(AgradRev, mathMixProb_student_t_qf) {
+  auto f = [](const auto& p, const auto& nu, const auto& sigma) {
+    const auto& p_scaled = stan::math::inv_logit(p);
+    const auto& nu_scaled = stan::math::exp(nu);
+    const auto& sigma_scaled = stan::math::exp(sigma);
+    return stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled);
+  };
+  auto f2 = [](const auto& nu, const auto& sigma) {
+    const auto& p_scaled = stan::math::inv_logit(0.8);
+    const auto& nu_scaled = stan::math::exp(nu);
+    const auto& sigma_scaled = stan::math::exp(sigma);
+    return stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled);
+  };
+  stan::test::expect_ad(f, 0.3, 0.5, 3);
+  stan::test::expect_ad(f, 0.8, 0.5, 0.1);
+  stan::test::expect_ad(f, 0.1, 3, 3);
+}
diff --git a/test/unit/math/prim/prob/student_t_qf_test.cpp b/test/unit/math/prim/prob/student_t_qf_test.cpp
new file mode 100644
index 00000000000..3c7005867ac
--- /dev/null
+++ b/test/unit/math/prim/prob/student_t_qf_test.cpp
@@ -0,0 +1,28 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <limits>
+
+TEST(MathFunctions, student_t_qf_vals) {
+  using stan::math::student_t_qf;
+
+  EXPECT_FLOAT_EQ(student_t_qf(0.2, 4, 2, 3), -0.8228937317);
+  EXPECT_FLOAT_EQ(student_t_qf(0.8, 4, 2, 3), 4.822893732);
+  EXPECT_FLOAT_EQ(student_t_qf(0.8, 3, 0, 1), 0.9784723124);
+
+  EXPECT_FLOAT_EQ(student_t_qf(0.5, 3, 0, 1), 0.0);
+}
+
+TEST(MathFunctions, student_t_qf_inf) {
+  using stan::math::student_t_qf;
+  long double log_p = std::numeric_limits<long double>::min();
+  const double inf = std::numeric_limits<double>::infinity();
+  EXPECT_EQ(student_t_qf(0, 2, 0, 1), -inf);
+  EXPECT_EQ(student_t_qf(1, 2, 0, 1), inf);
+}
+
+TEST(MathFunctions, student_t_qf_nan) {
+  using stan::math::student_t_qf;
+  double nan = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(student_t_qf(nan, 5, 0, 1), std::domain_error);
+  EXPECT_THROW(student_t_qf(2.1, 5, 0, 1), std::domain_error);
+}

From 44aed3ebeb97feda614d0e0d1fb7f26c50f7cd52 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 19 Jun 2025 18:00:23 +0800
Subject: [PATCH 4/9] Working scalar implementations

---
 stan/math/fwd/prob/student_t_qf.hpp           |  6 ++----
 stan/math/rev/prob/student_t_qf.hpp           |  6 ++----
 test/unit/math/mix/prob/student_t_qf_test.cpp | 12 ++++++++----
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
index 5ed093dd515..2b27a08259e 100644
--- a/stan/math/fwd/prob/student_t_qf.hpp
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -51,10 +51,8 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
 
   if constexpr (!is_constant<T_nu>::value) {
     const T_partials half_nu = nu_val / 2.0;
-    Eigen::Matrix<T_partials, Eigen::Dynamic, 1> hyper_arg_a(3);
-    hyper_arg_a << 0.5, half_nu, half_nu;
-    Eigen::Matrix<T_partials, Eigen::Dynamic, 1> hyper_arg_b(2);
-    hyper_arg_b << 1 + half_nu, 1 + half_nu;
+    Eigen::Matrix<T_partials, 3, 1> hyper_arg_a{0.5, half_nu, half_nu};
+    Eigen::Matrix<T_partials, 2, 1> hyper_arg_b{1 + half_nu, 1 + half_nu};
     const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
     const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
     const T_partials digamma_a1 = digamma(half_nu);
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index 53aadb1e3c3..96917ef0bd1 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -59,10 +59,8 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     if constexpr (!is_constant<T_nu>::value) {
       const double sigma_val = value_of(sigma);
       const double half_nu = nu_val / 2.0;
-      Eigen::VectorXd hyper_arg_a(3);
-      hyper_arg_a << 0.5, half_nu, half_nu;
-      Eigen::VectorXd hyper_arg_b(2);
-      hyper_arg_b << 1 + half_nu, 1 + half_nu;
+      Eigen::Vector3d hyper_arg_a{0.5, half_nu, half_nu};
+      Eigen::Vector2d hyper_arg_b{1 + half_nu, 1 + half_nu};
       const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
       const double hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
       const double digamma_a1 = digamma(half_nu);
diff --git a/test/unit/math/mix/prob/student_t_qf_test.cpp b/test/unit/math/mix/prob/student_t_qf_test.cpp
index 07c400a729d..51b33fc284c 100644
--- a/test/unit/math/mix/prob/student_t_qf_test.cpp
+++ b/test/unit/math/mix/prob/student_t_qf_test.cpp
@@ -10,13 +10,17 @@ TEST_F(AgradRev, mathMixProb_student_t_qf) {
     const auto& sigma_scaled = stan::math::exp(sigma);
     return stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled);
   };
-  auto f2 = [](const auto& nu, const auto& sigma) {
-    const auto& p_scaled = stan::math::inv_logit(0.8);
+  auto f2 = [](const auto& p, const auto& nu, const auto& mu) {
+    const auto& p_scaled = stan::math::inv_logit(p);
     const auto& nu_scaled = stan::math::exp(nu);
-    const auto& sigma_scaled = stan::math::exp(sigma);
-    return stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled);
+    return stan::math::student_t_qf(p_scaled, nu_scaled, mu, 2.5);
   };
+
   stan::test::expect_ad(f, 0.3, 0.5, 3);
   stan::test::expect_ad(f, 0.8, 0.5, 0.1);
   stan::test::expect_ad(f, 0.1, 3, 3);
+
+  stan::test::expect_ad(f2, 0.3, 0.5, 3);
+  stan::test::expect_ad(f2, 0.8, 0.5, 0.1);
+  stan::test::expect_ad(f2, 0.1, 3, 3);
 }

From bd9e001dcb969833f4e96a638a2008dbc69c9a19 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 19 Jun 2025 21:52:27 +0800
Subject: [PATCH 5/9] Proper container

---
 stan/math/prim/meta.hpp                       |  1 +
 stan/math/prim/meta/common_container_type.hpp | 72 +++++++++++++++++++
 stan/math/prim/prob/student_t_qf.hpp          | 39 ++++++++--
 test/unit/math/mix/prob/student_t_qf_test.cpp | 20 +++---
 .../unit/math/prim/prob/student_t_qf_test.cpp | 28 ++++++++
 5 files changed, 147 insertions(+), 13 deletions(-)
 create mode 100644 stan/math/prim/meta/common_container_type.hpp

diff --git a/stan/math/prim/meta.hpp b/stan/math/prim/meta.hpp
index dd6e3053eb6..5a36e92ba69 100644
--- a/stan/math/prim/meta.hpp
+++ b/stan/math/prim/meta.hpp
@@ -71,6 +71,7 @@
 #include <stan/math/prim/meta/ad_promotable.hpp>
 #include <stan/math/prim/meta/append_return_type.hpp>
 #include <stan/math/prim/meta/base_type.hpp>
+#include <stan/math/prim/meta/common_container_type.hpp>
 #include <stan/math/prim/meta/contains_std_vector.hpp>
 #include <stan/math/prim/meta/contains_tuple.hpp>
 #include <stan/math/prim/meta/error_index.hpp>
diff --git a/stan/math/prim/meta/common_container_type.hpp b/stan/math/prim/meta/common_container_type.hpp
new file mode 100644
index 00000000000..453c1ca0941
--- /dev/null
+++ b/stan/math/prim/meta/common_container_type.hpp
@@ -0,0 +1,72 @@
+#ifndef STAN_MATH_PRIM_META_COMMON_CONTAINER_TYPE_HPP
+#define STAN_MATH_PRIM_META_COMMON_CONTAINER_TYPE_HPP
+
+#include <stan/math/prim/meta/is_container.hpp>
+#include <stan/math/prim/meta/is_tuple.hpp>
+#include <stan/math/prim/meta/is_detected.hpp>
+#include <stan/math/prim/meta/is_stan_scalar.hpp>
+#include <stan/math/prim/meta/is_var_matrix.hpp>
+#include <stan/math/prim/meta/plain_type.hpp>
+#include <stan/math/prim/meta/return_type.hpp>
+#include <stan/math/prim/meta/promote_scalar_type.hpp>
+#include <type_traits>
+
+namespace stan {
+namespace internal {
+template <typename T1, typename T2, typename = void, typename = void>
+struct common_container_type_impl;
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_stan_scalar_t<T1>,
+                                  require_stan_scalar_t<T2>> {
+  using type = return_type_t<T1, T2>;
+};
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_container_t<T1>,
+                                  require_container_t<T2>> {
+  using return_t = return_type_t<T1, T2>;
+  using container_type_1 = math::promote_scalar_t<return_t, plain_type_t<T1>>;
+  using container_type_2 = math::promote_scalar_t<return_t, plain_type_t<T2>>;
+  using type =
+    std::conditional_t<
+      std::is_same<container_type_1, container_type_2>::value,
+      container_type_1,
+      void
+    >;
+};
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_stan_scalar_t<T1>,
+                                  require_container_t<T2>> {
+  using type = math::promote_scalar_t<return_type_t<T1, T2>, plain_type_t<T2>>;
+};
+
+template <typename T1, typename T2>
+struct common_container_type_impl<T1, T2, require_container_t<T1>,
+                                  require_stan_scalar_t<T2>> {
+  using type = math::promote_scalar_t<return_type_t<T1, T2>, plain_type_t<T1>>;
+};
+}
+
+template <typename... Ts>
+struct common_container_type;
+
+template <typename T>
+struct common_container_type<T> {
+  using type = typename internal::common_container_type_impl<
+      T, double>::type;  // Use double as a fallback for scalar types
+};
+
+template <typename T1, typename... Ts>
+struct common_container_type<T1, Ts...> {
+  using type =  typename internal::common_container_type_impl<
+      T1, typename common_container_type<Ts...>::type>::type;
+};
+
+template <typename... Ts>
+using common_container_t = typename common_container_type<Ts...>::type;
+
+}  // namespace stan
+
+#endif  // STAN_MATH_PRIM_META_PLAIN_TYPE_HPP
diff --git a/stan/math/prim/prob/student_t_qf.hpp b/stan/math/prim/prob/student_t_qf.hpp
index b84ccb37357..abbeb11f1ab 100644
--- a/stan/math/prim/prob/student_t_qf.hpp
+++ b/stan/math/prim/prob/student_t_qf.hpp
@@ -5,12 +5,16 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/sqrt.hpp>
 #include <stan/math/prim/fun/inv_inc_beta.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
 
 namespace stan {
 namespace math {
 
-double student_t_qf(const double p, const double nu, const double mu,
-                    const double sigma) {
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          require_all_stan_scalar_t<T_p, T_nu, T_mu, T_sigma>* = nullptr,
+          require_all_arithmetic_t<T_p, T_nu, T_mu, T_sigma>* = nullptr>
+double student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                    const T_sigma& sigma) {
   static constexpr const char* function = "student_t_qf";
   check_nonnegative(function, "Degrees of freedom parameter", nu);
   check_positive(function, "Scale parameter", sigma);
@@ -28,9 +32,36 @@ double student_t_qf(const double p, const double nu, const double mu,
   const double p_sign = p < 0.5 ? -1.0 : 1.0;
   const auto& ibeta_arg = inv_inc_beta(0.5 * nu, 0.5, 2 * p_val_flip);
 
-  return mu + p_sign * sigma * math::sqrt(nu)
-          * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+  return mu + p_sign * sigma * sqrt(nu) * sqrt(-1.0 + 1.0 / ibeta_arg);
 }
+
+template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
+          typename T_container = common_container_t<T_p, T_nu, T_mu, T_sigma>,
+          require_any_vector_t<T_p, T_nu, T_mu, T_sigma>* = nullptr,
+          require_not_t<std::is_void<T_container>>* = nullptr>
+auto student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
+                  const T_sigma& sigma) {
+  static constexpr const char* function = "student_t_qf";
+  const size_t max_size_all = max_size(p, nu, mu, sigma);
+  T_container result(max_size_all);
+
+  ref_type_t<T_p> p_ref = p;
+  ref_type_t<T_nu> nu_ref = nu;
+  ref_type_t<T_mu> mu_ref = mu;
+  ref_type_t<T_sigma> sigma_ref = sigma;
+
+  scalar_seq_view<ref_type_t<T_p>> p_vec(p_ref);
+  scalar_seq_view<ref_type_t<T_nu>> nu_vec(nu_ref);
+  scalar_seq_view<ref_type_t<T_mu>> mu_vec(mu_ref);
+  scalar_seq_view<ref_type_t<T_sigma>> sigma_vec(sigma_ref);
+
+  for (size_t i = 0; i < max_size_all; ++i) {
+    result[i] = student_t_qf(p_vec[i], nu_vec[i], mu_vec[i], sigma_vec[i]);
+  }
+
+  return result;
+}
+
 }
 }
 
diff --git a/test/unit/math/mix/prob/student_t_qf_test.cpp b/test/unit/math/mix/prob/student_t_qf_test.cpp
index 51b33fc284c..b82d81504c8 100644
--- a/test/unit/math/mix/prob/student_t_qf_test.cpp
+++ b/test/unit/math/mix/prob/student_t_qf_test.cpp
@@ -8,19 +8,21 @@ TEST_F(AgradRev, mathMixProb_student_t_qf) {
     const auto& p_scaled = stan::math::inv_logit(p);
     const auto& nu_scaled = stan::math::exp(nu);
     const auto& sigma_scaled = stan::math::exp(sigma);
-    return stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled);
-  };
-  auto f2 = [](const auto& p, const auto& nu, const auto& mu) {
-    const auto& p_scaled = stan::math::inv_logit(p);
-    const auto& nu_scaled = stan::math::exp(nu);
-    return stan::math::student_t_qf(p_scaled, nu_scaled, mu, 2.5);
+    return stan::math::sum(stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled));
   };
 
   stan::test::expect_ad(f, 0.3, 0.5, 3);
   stan::test::expect_ad(f, 0.8, 0.5, 0.1);
   stan::test::expect_ad(f, 0.1, 3, 3);
 
-  stan::test::expect_ad(f2, 0.3, 0.5, 3);
-  stan::test::expect_ad(f2, 0.8, 0.5, 0.1);
-  stan::test::expect_ad(f2, 0.1, 3, 3);
+  Eigen::VectorXd p(3);
+  p << 0.3, 0.8, 0.1;
+
+  Eigen::VectorXd nu(3);
+  nu << 0.5, 0.5, 3;
+
+  Eigen::VectorXd sigma(3);
+  sigma << 3, 0.1, 3;
+
+  stan::test::expect_ad(f, p, nu, sigma);
 }
diff --git a/test/unit/math/prim/prob/student_t_qf_test.cpp b/test/unit/math/prim/prob/student_t_qf_test.cpp
index 3c7005867ac..67009d0c771 100644
--- a/test/unit/math/prim/prob/student_t_qf_test.cpp
+++ b/test/unit/math/prim/prob/student_t_qf_test.cpp
@@ -1,10 +1,16 @@
 #include <stan/math/prim.hpp>
+#include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 #include <limits>
 
 TEST(MathFunctions, student_t_qf_vals) {
   using stan::math::student_t_qf;
 
+  Eigen::VectorXd p(5);
+  p << 0.1, 0.2, 0.5, 0.8, 0.9;
+
+  Eigen::VectorXd res = student_t_qf(p, 4, 2, 3);
+
   EXPECT_FLOAT_EQ(student_t_qf(0.2, 4, 2, 3), -0.8228937317);
   EXPECT_FLOAT_EQ(student_t_qf(0.8, 4, 2, 3), 4.822893732);
   EXPECT_FLOAT_EQ(student_t_qf(0.8, 3, 0, 1), 0.9784723124);
@@ -12,6 +18,28 @@ TEST(MathFunctions, student_t_qf_vals) {
   EXPECT_FLOAT_EQ(student_t_qf(0.5, 3, 0, 1), 0.0);
 }
 
+TEST(MathFunctions, student_t_qf_vec) {
+  using stan::math::student_t_qf;
+
+  Eigen::VectorXd p(5);
+  p << 0.1, 0.2, 0.5, 0.8, 0.9;
+
+  Eigen::VectorXd nu(5);
+  nu << 4, 4, 3, 3, 3;
+
+  Eigen::VectorXd res(5);
+  for (int i = 0; i < 5; ++i) {
+    res[i] = student_t_qf(p[i], nu[i], 2, 3);
+  }
+  EXPECT_MATRIX_FLOAT_EQ(res, student_t_qf(p, nu, 2, 3));
+
+  Eigen::RowVectorXd p_rowvec = p.transpose();
+  Eigen::RowVectorXd nu_rowvec = nu.transpose();
+  Eigen::RowVectorXd res_rowvec = res.transpose();
+
+  EXPECT_MATRIX_FLOAT_EQ(res_rowvec, student_t_qf(p_rowvec, nu_rowvec, 2, 3));
+}
+
 TEST(MathFunctions, student_t_qf_inf) {
   using stan::math::student_t_qf;
   long double log_p = std::numeric_limits<long double>::min();

From 774cd41c504ffca48a6c2381a1921a4316511a9b Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Thu, 19 Jun 2025 22:12:32 +0800
Subject: [PATCH 6/9] Update doc

---
 stan/math/fwd/prob/student_t_qf.hpp           | 14 +++++----
 stan/math/prim/meta/common_container_type.hpp | 12 +++++--
 stan/math/prim/prob/student_t_qf.hpp          | 31 +++++++++++++++++++
 stan/math/rev/prob/student_t_qf.hpp           | 30 ++++++++++++------
 4 files changed, 69 insertions(+), 18 deletions(-)

diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
index 2b27a08259e..a03c052ee6f 100644
--- a/stan/math/fwd/prob/student_t_qf.hpp
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -35,7 +35,6 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     return {mu_val, 0.0};
   }
 
-
   const T_partials p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
   const double p_sign = value_of_rec(p_val) < 0.5 ? -1.0 : 1.0;
 
@@ -53,13 +52,16 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     const T_partials half_nu = nu_val / 2.0;
     Eigen::Matrix<T_partials, 3, 1> hyper_arg_a{0.5, half_nu, half_nu};
     Eigen::Matrix<T_partials, 2, 1> hyper_arg_b{1 + half_nu, 1 + half_nu};
-    const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
-    const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
+    const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
+                                                    ibeta_arg);
+    const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2,
+                                                   (2 + nu_val) / 2, ibeta_arg);
     const T_partials digamma_a1 = digamma(half_nu);
     const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
 
     const T_partials arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
-    const T_partials arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg) * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+    const T_partials arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg)
+                                * (log(ibeta_arg) - digamma_a1 + digamma_a2);
 
     const T_partials num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
     const T_partials den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
@@ -71,11 +73,11 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
   }
 
   if constexpr (!is_constant<T_sigma>::value) {
-    rtn.d_ += sigma.d_ * p_sign * sqrt(nu_val) * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+    rtn.d_ += sigma.d_ * p_sign * sqrt(nu_val)
+                       * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
   return rtn;
-
 }
 }
 }
diff --git a/stan/math/prim/meta/common_container_type.hpp b/stan/math/prim/meta/common_container_type.hpp
index 453c1ca0941..fcfd09d45dd 100644
--- a/stan/math/prim/meta/common_container_type.hpp
+++ b/stan/math/prim/meta/common_container_type.hpp
@@ -55,12 +55,20 @@ struct common_container_type;
 template <typename T>
 struct common_container_type<T> {
   using type = typename internal::common_container_type_impl<
-      T, double>::type;  // Use double as a fallback for scalar types
+      T, double>::type;  // Use double for base case
 };
 
+/**
+ * Determine the common container type for a variadic list of types.
+ * If all types are scalars, then the common scalar type is returned.
+ * If all container types the same, but not necessarily the same scalar type,
+ * the common container type with the common scalar type is returned.
+ *
+ * If different container types are present, the result is `void`.
+ */
 template <typename T1, typename... Ts>
 struct common_container_type<T1, Ts...> {
-  using type =  typename internal::common_container_type_impl<
+  using type = typename internal::common_container_type_impl<
       T1, typename common_container_type<Ts...>::type>::type;
 };
 
diff --git a/stan/math/prim/prob/student_t_qf.hpp b/stan/math/prim/prob/student_t_qf.hpp
index abbeb11f1ab..f4f018aefd2 100644
--- a/stan/math/prim/prob/student_t_qf.hpp
+++ b/stan/math/prim/prob/student_t_qf.hpp
@@ -10,6 +10,21 @@
 namespace stan {
 namespace math {
 
+/**
+ * The quantile function of the Student's t-distribution.
+ *
+ * @tparam T_p type of the probability parameter
+ * @tparam T_nu type of the degrees of freedom parameter
+ * @tparam T_mu type of the location parameter
+ * @tparam T_sigma type of the scale parameter
+ * @param p Probability in the range [0, 1].
+ * @param nu Degrees of freedom, must be non-negative.
+ * @param mu Location parameter.
+ * @param sigma Scale parameter, must be positive.
+ * @return Quantile function value.
+ * @throw std::domain_error if `nu` is negative or `sigma` is not positive,
+ * or if `p` is not in [0, 1].
+ */
 template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
           require_all_stan_scalar_t<T_p, T_nu, T_mu, T_sigma>* = nullptr,
           require_all_arithmetic_t<T_p, T_nu, T_mu, T_sigma>* = nullptr>
@@ -35,6 +50,22 @@ double student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
   return mu + p_sign * sigma * sqrt(nu) * sqrt(-1.0 + 1.0 / ibeta_arg);
 }
 
+/**
+ * A vectorized version of the Student's t quantile function that accepts
+ * std::vectors, Eigen Matrix/Array objects, or expressions, and containers of
+ * these.
+ *
+ * @tparam T_p type of the probability parameter
+ * @tparam T_nu type of the degrees of freedom parameter
+ * @tparam T_mu type of the location parameter
+ * @tparam T_sigma type of the scale parameter
+ * @tparam T_container type of the container to hold results
+ * @param p Probability in the range [0, 1].
+ * @param nu Degrees of freedom, must be non-negative.
+ * @param mu Location parameter.
+ * @param sigma Scale parameter, must be positive.
+ * @return Container with quantile function values for each input.
+ */
 template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
           typename T_container = common_container_t<T_p, T_nu, T_mu, T_sigma>,
           require_any_vector_t<T_p, T_nu, T_mu, T_sigma>* = nullptr,
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index 96917ef0bd1..3416abb3a35 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -34,6 +34,7 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
   } else if (p == 0.5) {
     ret_val = mu_val;
   } else {
+    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) -1.0);
     const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
     const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
     ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
@@ -42,49 +43,58 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
                   * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
-  return make_callback_var(ret_val, [p, mu, sigma, nu, ibeta_arg](auto& vi) mutable {
+  return make_callback_var(ret_val,
+    [p, mu, sigma, nu, ibeta_arg](auto& vi) mutable {
     const double p_val = value_of(p);
     const double mu_val = value_of(mu);
     const double sigma_val = value_of(sigma);
     const double nu_val = value_of(nu);
     const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
     const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
+    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) -1.0);
 
     if constexpr (!is_constant<T_p>::value) {
       if (p.val() > 0.0 && p.val() < 1.0) {
-        p.adj() += vi.adj() * exp(-student_t_lpdf(vi.val(), nu_val, mu_val, sigma_val));
+        p.adj() += vi.adj()
+                    * exp(-student_t_lpdf(vi.val(), nu_val, mu_val, sigma_val));
       }
     }
 
     if constexpr (!is_constant<T_nu>::value) {
-      const double sigma_val = value_of(sigma);
       const double half_nu = nu_val / 2.0;
+      const double nu_p1_div_2 = (1 + nu_val) / 2.0;
       Eigen::Vector3d hyper_arg_a{0.5, half_nu, half_nu};
       Eigen::Vector2d hyper_arg_b{1 + half_nu, 1 + half_nu};
-      const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
-      const double hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
+      const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
+                                                  ibeta_arg);
+
+      const double hyper2f1 = hypergeometric_2F1(1.0, nu_p1_div_2,
+                                                (2.0 + nu_val) / 2.0,
+                                                ibeta_arg);
       const double digamma_a1 = digamma(half_nu);
-      const double digamma_a2 = digamma((1 + nu_val) / 2);
+      const double digamma_a2 = digamma(nu_p1_div_2);
 
       const double arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
-      const double arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg) * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+      const double arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg)
+                              * (log(ibeta_arg) - digamma_a1 + digamma_a2);
 
       const double num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
-      const double den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
+      const double den1 = 4 * sqrt(nu_val) * sqrt_inv_ibeta_m1;
       nu.adj() += vi.adj() * p_sign * num1 / den1;
     }
+
     if constexpr (!is_constant<T_mu>::value) {
       if (p_val > 0.0 && p_val < 1.0) {
         mu.adj() += vi.adj();
       }
     }
+
     if constexpr (!is_constant<T_sigma>::value) {
       if (p_val > 0.0 && p_val < 1.0) {
-        sigma.adj() += vi.adj() * p_sign * sqrt(value_of(nu)) * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+        sigma.adj() += vi.adj() * p_sign * sqrt(nu_val) * sqrt_inv_ibeta_m1;
       }
     }
   });
-
 }
 }
 }

From f0174e65f045617000a1c7a2c5b8fa8122e29c0c Mon Sep 17 00:00:00 2001
From: Stan Jenkins <mc.stanislaw@gmail.com>
Date: Thu, 19 Jun 2025 10:25:03 -0400
Subject: [PATCH 7/9] [Jenkins] auto-formatting by clang-format version
 10.0.0-4ubuntu1

---
 stan/math/fwd/prob/student_t_qf.hpp           | 23 ++++++-------
 stan/math/prim/meta/common_container_type.hpp | 11 +++----
 stan/math/prim/prob/student_t_qf.hpp          |  4 +--
 stan/math/rev/prob/student_t_qf.hpp           | 32 +++++++++----------
 test/unit/math/mix/prob/student_t_qf_test.cpp |  3 +-
 5 files changed, 36 insertions(+), 37 deletions(-)

diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
index a03c052ee6f..960c90ba489 100644
--- a/stan/math/fwd/prob/student_t_qf.hpp
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -14,7 +14,7 @@ template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
           require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
           require_fvar_t<FvarT>* = nullptr>
 FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
-                    const T_sigma& sigma) {
+                   const T_sigma& sigma) {
   static constexpr const char* function = "student_t_qf";
   using T_partials = partials_type_t<FvarT>;
 
@@ -39,8 +39,9 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
   const double p_sign = value_of_rec(p_val) < 0.5 ? -1.0 : 1.0;
 
   T_partials ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
-  T_partials rtn_val = mu_val + p_sign * sigma_val * sqrt(nu_val)
-                * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+  T_partials rtn_val = mu_val
+                       + p_sign * sigma_val * sqrt(nu_val)
+                             * math::sqrt(-1.0 + 1.0 / ibeta_arg);
 
   FvarT rtn(rtn_val, 0.0);
 
@@ -52,16 +53,16 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     const T_partials half_nu = nu_val / 2.0;
     Eigen::Matrix<T_partials, 3, 1> hyper_arg_a{0.5, half_nu, half_nu};
     Eigen::Matrix<T_partials, 2, 1> hyper_arg_b{1 + half_nu, 1 + half_nu};
-    const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
-                                                    ibeta_arg);
-    const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2,
-                                                   (2 + nu_val) / 2, ibeta_arg);
+    const T_partials hyper_arg
+        = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
+    const T_partials hyper2f1
+        = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
     const T_partials digamma_a1 = digamma(half_nu);
     const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
 
     const T_partials arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
     const T_partials arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg)
-                                * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+                             * (log(ibeta_arg) - digamma_a1 + digamma_a2);
 
     const T_partials num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
     const T_partials den1 = 4 * sqrt(nu_val) * sqrt(-1 + 1 / ibeta_arg);
@@ -74,12 +75,12 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
 
   if constexpr (!is_constant<T_sigma>::value) {
     rtn.d_ += sigma.d_ * p_sign * sqrt(nu_val)
-                       * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+              * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
   return rtn;
 }
-}
-}
+}  // namespace math
+}  // namespace stan
 
 #endif
diff --git a/stan/math/prim/meta/common_container_type.hpp b/stan/math/prim/meta/common_container_type.hpp
index fcfd09d45dd..df74fddff5a 100644
--- a/stan/math/prim/meta/common_container_type.hpp
+++ b/stan/math/prim/meta/common_container_type.hpp
@@ -28,12 +28,9 @@ struct common_container_type_impl<T1, T2, require_container_t<T1>,
   using return_t = return_type_t<T1, T2>;
   using container_type_1 = math::promote_scalar_t<return_t, plain_type_t<T1>>;
   using container_type_2 = math::promote_scalar_t<return_t, plain_type_t<T2>>;
-  using type =
-    std::conditional_t<
-      std::is_same<container_type_1, container_type_2>::value,
-      container_type_1,
-      void
-    >;
+  using type = std::conditional_t<
+      std::is_same<container_type_1, container_type_2>::value, container_type_1,
+      void>;
 };
 
 template <typename T1, typename T2>
@@ -47,7 +44,7 @@ struct common_container_type_impl<T1, T2, require_container_t<T1>,
                                   require_stan_scalar_t<T2>> {
   using type = math::promote_scalar_t<return_type_t<T1, T2>, plain_type_t<T1>>;
 };
-}
+}  // namespace internal
 
 template <typename... Ts>
 struct common_container_type;
diff --git a/stan/math/prim/prob/student_t_qf.hpp b/stan/math/prim/prob/student_t_qf.hpp
index f4f018aefd2..8a134498d75 100644
--- a/stan/math/prim/prob/student_t_qf.hpp
+++ b/stan/math/prim/prob/student_t_qf.hpp
@@ -93,7 +93,7 @@ auto student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
   return result;
 }
 
-}
-}
+}  // namespace math
+}  // namespace stan
 
 #endif
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index 3416abb3a35..cab5550082f 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -13,7 +13,7 @@ template <typename T_p, typename T_nu, typename T_mu, typename T_sigma,
           require_all_stan_scalar_t<T_p, T_mu, T_sigma, T_nu>* = nullptr,
           require_any_var_t<T_p, T_mu, T_sigma, T_nu>* = nullptr>
 var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
-                    const T_sigma& sigma) {
+                 const T_sigma& sigma) {
   static constexpr const char* function = "student_t_qf";
   const double p_val = value_of(p);
   const double nu_val = value_of(nu);
@@ -34,29 +34,30 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
   } else if (p == 0.5) {
     ret_val = mu_val;
   } else {
-    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) -1.0);
+    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) - 1.0);
     const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
     const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
     ibeta_arg = inv_inc_beta(0.5 * nu_val, 0.5, 2 * p_val_flip);
 
-    ret_val = mu_val + p_sign * sigma_val * math::sqrt(nu_val)
-                  * math::sqrt(-1.0 + 1.0 / ibeta_arg);
+    ret_val = mu_val
+              + p_sign * sigma_val * math::sqrt(nu_val)
+                    * math::sqrt(-1.0 + 1.0 / ibeta_arg);
   }
 
-  return make_callback_var(ret_val,
-    [p, mu, sigma, nu, ibeta_arg](auto& vi) mutable {
+  return make_callback_var(ret_val, [p, mu, sigma, nu,
+                                     ibeta_arg](auto& vi) mutable {
     const double p_val = value_of(p);
     const double mu_val = value_of(mu);
     const double sigma_val = value_of(sigma);
     const double nu_val = value_of(nu);
     const double p_val_flip = p_val < 0.5 ? p_val : 1.0 - p_val;
     const double p_sign = p_val < 0.5 ? -1.0 : 1.0;
-    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) -1.0);
+    const double sqrt_inv_ibeta_m1 = sqrt(inv(ibeta_arg) - 1.0);
 
     if constexpr (!is_constant<T_p>::value) {
       if (p.val() > 0.0 && p.val() < 1.0) {
         p.adj() += vi.adj()
-                    * exp(-student_t_lpdf(vi.val(), nu_val, mu_val, sigma_val));
+                   * exp(-student_t_lpdf(vi.val(), nu_val, mu_val, sigma_val));
       }
     }
 
@@ -65,18 +66,17 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
       const double nu_p1_div_2 = (1 + nu_val) / 2.0;
       Eigen::Vector3d hyper_arg_a{0.5, half_nu, half_nu};
       Eigen::Vector2d hyper_arg_b{1 + half_nu, 1 + half_nu};
-      const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
-                                                  ibeta_arg);
+      const double hyper_arg
+          = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
 
-      const double hyper2f1 = hypergeometric_2F1(1.0, nu_p1_div_2,
-                                                (2.0 + nu_val) / 2.0,
-                                                ibeta_arg);
+      const double hyper2f1 = hypergeometric_2F1(
+          1.0, nu_p1_div_2, (2.0 + nu_val) / 2.0, ibeta_arg);
       const double digamma_a1 = digamma(half_nu);
       const double digamma_a2 = digamma(nu_p1_div_2);
 
       const double arg_1 = (4 * hyper_arg * sqrt(1 - ibeta_arg)) / nu_val;
       const double arg_2 = -2 * hyper2f1 * (-1 + ibeta_arg)
-                              * (log(ibeta_arg) - digamma_a1 + digamma_a2);
+                           * (log(ibeta_arg) - digamma_a1 + digamma_a2);
 
       const double num1 = sigma_val * (-2 + (2 - arg_1 + arg_2) / ibeta_arg);
       const double den1 = 4 * sqrt(nu_val) * sqrt_inv_ibeta_m1;
@@ -96,7 +96,7 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     }
   });
 }
-}
-}
+}  // namespace math
+}  // namespace stan
 
 #endif
diff --git a/test/unit/math/mix/prob/student_t_qf_test.cpp b/test/unit/math/mix/prob/student_t_qf_test.cpp
index b82d81504c8..1e5dc213b3c 100644
--- a/test/unit/math/mix/prob/student_t_qf_test.cpp
+++ b/test/unit/math/mix/prob/student_t_qf_test.cpp
@@ -8,7 +8,8 @@ TEST_F(AgradRev, mathMixProb_student_t_qf) {
     const auto& p_scaled = stan::math::inv_logit(p);
     const auto& nu_scaled = stan::math::exp(nu);
     const auto& sigma_scaled = stan::math::exp(sigma);
-    return stan::math::sum(stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled));
+    return stan::math::sum(
+        stan::math::student_t_qf(p_scaled, nu_scaled, 0, sigma_scaled));
   };
 
   stan::test::expect_ad(f, 0.3, 0.5, 3);

From 63004a792b09398bd2cbf27c5a8639cd82e1380d Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Fri, 20 Jun 2025 09:40:45 +0800
Subject: [PATCH 8/9] Dynamic

---
 stan/math/fwd/prob/student_t_qf.hpp | 14 ++++++++------
 stan/math/rev/prob/student_t_qf.hpp | 10 ++++++----
 2 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
index 960c90ba489..b3fbf78d9e3 100644
--- a/stan/math/fwd/prob/student_t_qf.hpp
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -51,12 +51,14 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
 
   if constexpr (!is_constant<T_nu>::value) {
     const T_partials half_nu = nu_val / 2.0;
-    Eigen::Matrix<T_partials, 3, 1> hyper_arg_a{0.5, half_nu, half_nu};
-    Eigen::Matrix<T_partials, 2, 1> hyper_arg_b{1 + half_nu, 1 + half_nu};
-    const T_partials hyper_arg
-        = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
-    const T_partials hyper2f1
-        = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
+    Eigen::Matrix<T_partials, -1, 1> hyper_arg_a(3);
+    hyper_arg_a << 0.5, half_nu, half_nu;
+    Eigen::Matrix<T_partials, -1, 1> hyper_arg_b(2);
+    hyper_arg_b << 1 + half_nu, 1 + half_nu;
+    const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
+                                                    ibeta_arg);
+    const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2,
+                                                   (2 + nu_val) / 2, ibeta_arg);
     const T_partials digamma_a1 = digamma(half_nu);
     const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
 
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index cab5550082f..b362b1c26a4 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -64,10 +64,12 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     if constexpr (!is_constant<T_nu>::value) {
       const double half_nu = nu_val / 2.0;
       const double nu_p1_div_2 = (1 + nu_val) / 2.0;
-      Eigen::Vector3d hyper_arg_a{0.5, half_nu, half_nu};
-      Eigen::Vector2d hyper_arg_b{1 + half_nu, 1 + half_nu};
-      const double hyper_arg
-          = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
+      Eigen::VectorXd hyper_arg_a(3);
+      hyper_arg_a << 0.5, half_nu, half_nu;
+      Eigen::VectorXd hyper_arg_b(2);
+      hyper_arg_b << 1 + half_nu, 1 + half_nu;
+      const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
+                                                  ibeta_arg);
 
       const double hyper2f1 = hypergeometric_2F1(
           1.0, nu_p1_div_2, (2.0 + nu_val) / 2.0, ibeta_arg);

From 332c74f8365a26a44673cb13855b27c108ad569f Mon Sep 17 00:00:00 2001
From: Stan Jenkins <mc.stanislaw@gmail.com>
Date: Thu, 19 Jun 2025 21:41:48 -0400
Subject: [PATCH 9/9] [Jenkins] auto-formatting by clang-format version
 10.0.0-4ubuntu1

---
 stan/math/fwd/prob/student_t_qf.hpp | 8 ++++----
 stan/math/rev/prob/student_t_qf.hpp | 4 ++--
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/stan/math/fwd/prob/student_t_qf.hpp b/stan/math/fwd/prob/student_t_qf.hpp
index b3fbf78d9e3..0b87c60171d 100644
--- a/stan/math/fwd/prob/student_t_qf.hpp
+++ b/stan/math/fwd/prob/student_t_qf.hpp
@@ -55,10 +55,10 @@ FvarT student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
     hyper_arg_a << 0.5, half_nu, half_nu;
     Eigen::Matrix<T_partials, -1, 1> hyper_arg_b(2);
     hyper_arg_b << 1 + half_nu, 1 + half_nu;
-    const T_partials hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
-                                                    ibeta_arg);
-    const T_partials hyper2f1 = hypergeometric_2F1(1, (1 + nu_val) / 2,
-                                                   (2 + nu_val) / 2, ibeta_arg);
+    const T_partials hyper_arg
+        = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
+    const T_partials hyper2f1
+        = hypergeometric_2F1(1, (1 + nu_val) / 2, (2 + nu_val) / 2, ibeta_arg);
     const T_partials digamma_a1 = digamma(half_nu);
     const T_partials digamma_a2 = digamma((1 + nu_val) / 2);
 
diff --git a/stan/math/rev/prob/student_t_qf.hpp b/stan/math/rev/prob/student_t_qf.hpp
index b362b1c26a4..baaae8377e2 100644
--- a/stan/math/rev/prob/student_t_qf.hpp
+++ b/stan/math/rev/prob/student_t_qf.hpp
@@ -68,8 +68,8 @@ var student_t_qf(const T_p& p, const T_nu& nu, const T_mu& mu,
       hyper_arg_a << 0.5, half_nu, half_nu;
       Eigen::VectorXd hyper_arg_b(2);
       hyper_arg_b << 1 + half_nu, 1 + half_nu;
-      const double hyper_arg = hypergeometric_pFq(hyper_arg_a, hyper_arg_b,
-                                                  ibeta_arg);
+      const double hyper_arg
+          = hypergeometric_pFq(hyper_arg_a, hyper_arg_b, ibeta_arg);
 
       const double hyper2f1 = hypergeometric_2F1(
           1.0, nu_p1_div_2, (2.0 + nu_val) / 2.0, ibeta_arg);