diff --git a/stan/math/opencl/matrix_cl.hpp b/stan/math/opencl/matrix_cl.hpp
index ded00e3a821..5248d719b93 100644
--- a/stan/math/opencl/matrix_cl.hpp
+++ b/stan/math/opencl/matrix_cl.hpp
@@ -272,7 +272,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   matrix_cl(const int rows, const int cols,
             matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(rows), cols_(cols), view_(partial_view) {
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     cl::Context& ctx = opencl_context.context();
@@ -310,7 +310,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
                      matrix_cl_view partial_view = matrix_cl_view::Entire)
       : rows_(A.rows()), cols_(A.cols()), view_(partial_view) {
     using Mat_type = std::decay_t<ref_type_for_opencl_t<Mat>>;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer_no_heap_if<
@@ -489,7 +489,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(const T* A) {
     cl::Event transfer_event;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -509,7 +509,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool in_order = false>
   cl::Event initialize_buffer(T* A) {
     cl::Event transfer_event;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return transfer_event;
     }
     cl::Context& ctx = opencl_context.context();
@@ -548,7 +548,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
    */
   template <bool No_heap, typename U, std::enable_if_t<No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     initialize_buffer(obj.data());
@@ -558,7 +558,7 @@ class matrix_cl<T, require_arithmetic_t<T>> : public matrix_cl_base {
   template <bool No_heap, typename U, std::enable_if_t<!No_heap>* = nullptr>
   void initialize_buffer_no_heap_if(U&& obj) {
     using U_val = std::decay_t<ref_type_for_opencl_t<U>>;
-    if (size() == 0) {
+    if (this->size() == 0) {
       return;
     }
     auto* obj_heap = new U_val(std::move(obj));
diff --git a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
index 76a342408df..47228286ac2 100644
--- a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
@@ -66,13 +66,13 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> bernoulli_logit_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
 
   if (N == 0) {
diff --git a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
index 8db546a1853..accc947ab17 100644
--- a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp
@@ -36,7 +36,7 @@ return_type_t<T_prob_cl> bernoulli_logit_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_lpmf.hpp b/stan/math/opencl/prim/bernoulli_lpmf.hpp
index b6335842b60..ddbfb5e9074 100644
--- a/stan/math/opencl/prim/bernoulli_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lpmf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n,
                          "Probability parameter", theta);
-  const size_t N = is_n_vector ? size(n) : size(theta);
+  const size_t N = is_n_vector ? math::size(n) : math::size(theta);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
index df3ecc809f8..b6be840c7ec 100644
--- a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
@@ -60,10 +60,10 @@ return_type_t<T_x, T_alpha, T_beta> categorical_logit_glm_lpmf(
   static const char* function = "categorical_logit_glm_lpmf";
   if (is_y_vector) {
     check_size_match(function, "Rows of ", "x", N_instances, "size of ", "y",
-                     size(y));
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "beta", N_classes, "size of ",
-                   "alpha", size(alpha));
+                   "alpha", math::size(alpha));
   check_size_match(function, "Columns of ", "x", N_attributes, "Rows of",
                    "beta", beta.rows());
 
diff --git a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
index 1f4aad75e88..35a3e734ab0 100644
--- a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp
@@ -82,17 +82,17 @@ neg_binomial_2_log_glm_lpmf(const T_y_cl& y, const T_x_cl& x,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_phi_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "phi",
-                     size(phi));
+                     math::size(phi));
   }
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (N == 0) {
     return 0;
diff --git a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
index 786a6f6cd2a..2d98f29a163 100644
--- a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
+++ b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
@@ -74,17 +74,17 @@ normal_id_glm_lpdf(const T_y_cl& y, const T_x_cl& x, const T_alpha_cl& alpha,
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_sigma_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "sigma",
-                     size(sigma));
+                     math::size(sigma));
   }
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (!include_summand<propto, T_y_cl, T_x_cl, T_alpha_cl, T_beta_cl,
                        T_sigma_cl>::value) {
diff --git a/stan/math/opencl/prim/num_elements.hpp b/stan/math/opencl/prim/num_elements.hpp
index a5d965290d3..aedae05f616 100644
--- a/stan/math/opencl/prim/num_elements.hpp
+++ b/stan/math/opencl/prim/num_elements.hpp
@@ -3,6 +3,7 @@
 #ifdef STAN_OPENCL
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/opencl/prim/size.hpp>
 
 namespace stan {
 namespace math {
@@ -16,7 +17,7 @@ namespace math {
 template <typename T,
           require_nonscalar_prim_or_rev_kernel_expression_t<T>* = nullptr>
 size_t num_elements(const T& m) {
-  return size(m);
+  return math::size(m);
 }
 
 }  // namespace math
diff --git a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
index 71726353e90..5fb69cd6047 100644
--- a/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/ordered_logistic_glm_lpmf.hpp
@@ -63,19 +63,19 @@ return_type_t<T_x, T_beta, T_cuts> ordered_logistic_glm_lpmf(
 
   const size_t N_instances = x.rows();
   const size_t N_attributes = x.cols();
-  const size_t N_classes = size(cuts) + 1;
+  const size_t N_classes = math::size(cuts) + 1;
 
   if (is_y_vector) {
     check_size_match(function, "Rows of ", "x", N_instances, "rows of ", "y",
-                     size(y));
+                     math::size(y));
   }
   check_size_match(function, "Columns of ", "x", N_attributes, "Size of",
-                   "beta", size(beta));
+                   "beta", math::size(beta));
 
   const auto& cuts_val = eval(value_of(cuts));
   if (N_classes >= 2) {
-    auto cuts_head = block_zero_based(cuts_val, 0, 0, size(cuts) - 1, 1);
-    auto cuts_tail = block_zero_based(cuts_val, 1, 0, size(cuts) - 1, 1);
+    auto cuts_head = block_zero_based(cuts_val, 0, 0, math::size(cuts) - 1, 1);
+    auto cuts_tail = block_zero_based(cuts_val, 1, 0, math::size(cuts) - 1, 1);
     check_cl(function, "Cuts", cuts_head, "ordered and finite")
         = cuts_head < cuts_tail && isfinite(cuts_head) && isfinite(cuts_tail);
   } else {
diff --git a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
index 47ea9dc69ae..16e4774eb32 100644
--- a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
@@ -66,13 +66,13 @@ return_type_t<T_x_cl, T_alpha_cl, T_beta_cl> poisson_log_glm_lpmf(
   const size_t M = x.cols();
 
   if (is_y_vector) {
-    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", size(y));
+    check_size_match(function, "Rows of ", "x", N, "rows of ", "y", math::size(y));
   }
   check_size_match(function, "Columns of ", "x_cl", M, "size of ", "beta",
-                   size(beta));
+                   math::size(beta));
   if (is_alpha_vector) {
     check_size_match(function, "Rows of ", "x", N, "size of ", "alpha",
-                     size(alpha));
+                     math::size(alpha));
   }
   if (N == 0) {
     return 0;
diff --git a/stan/math/opencl/prim/poisson_log_lpmf.hpp b/stan/math/opencl/prim/poisson_log_lpmf.hpp
index 515a311a690..c2fe7d1a6a6 100644
--- a/stan/math/opencl/prim/poisson_log_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_lpmf.hpp
@@ -39,7 +39,7 @@ return_type_t<T_log_rate_cl> poisson_log_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n, "Log rate parameter",
                          alpha);
-  const size_t N = is_n_vector ? size(n) : size(alpha);
+  const size_t N = is_n_vector ? math::size(n) : math::size(alpha);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/poisson_lpmf.hpp b/stan/math/opencl/prim/poisson_lpmf.hpp
index 896ef4d2209..e4b047f0f34 100644
--- a/stan/math/opencl/prim/poisson_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_lpmf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_rate_cl> poisson_lpmf(const T_n_cl& n,
 
   check_consistent_sizes(function, "Random variable", n, "Rate parameter",
                          lambda);
-  const size_t N = is_n_vector ? size(n) : size(lambda);
+  const size_t N = is_n_vector ? math::size(n) : math::size(lambda);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_lpdf.hpp b/stan/math/opencl/prim/std_normal_lpdf.hpp
index b9782db683e..2bf4de436c1 100644
--- a/stan/math/opencl/prim/std_normal_lpdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lpdf.hpp
@@ -35,7 +35,7 @@ inline return_type_t<T_y_cl> std_normal_lpdf(const T_y_cl& y) {
   using std::isfinite;
   using std::isnan;
 
-  const size_t N = size(y);
+  const size_t N = math::size(y);
   if (N == 0) {
     return 0.0;
   }
diff --git a/stan/math/opencl/zeros_strict_tri.hpp b/stan/math/opencl/zeros_strict_tri.hpp
index 66ff85c64f4..abf7349339a 100644
--- a/stan/math/opencl/zeros_strict_tri.hpp
+++ b/stan/math/opencl/zeros_strict_tri.hpp
@@ -40,7 +40,7 @@ inline void matrix_cl<T, require_arithmetic_t<T>>::zeros_strict_tri() try {
         "zeros_strict_tri", "matrix_view",
         "matrix_cl_view::Diagonal is not a valid template parameter value", "");
   }
-  if (size() == 0) {
+  if (this->size() == 0) {
     return;
   }
   this->view_ = both(this->view_, invert(matrix_view));
diff --git a/stan/math/prim/fun/offset_multiplier_constrain.hpp b/stan/math/prim/fun/offset_multiplier_constrain.hpp
index 5f8a984a949..c724610deda 100644
--- a/stan/math/prim/fun/offset_multiplier_constrain.hpp
+++ b/stan/math/prim/fun/offset_multiplier_constrain.hpp
@@ -90,7 +90,7 @@ inline return_type_t<T, M, S> offset_multiplier_constrain(const T& x,
     return mu + x;
   }
   check_positive_finite("offset_multiplier_constrain", "multiplier", sigma);
-  lp += multiply_log(size(x), sigma);
+  lp += multiply_log(math::size(x), sigma);
   return fma(sigma, x, mu);
 }
 
diff --git a/stan/math/prim/fun/size_zero.hpp b/stan/math/prim/fun/size_zero.hpp
index 8af469eee9f..c6aa81da3fa 100644
--- a/stan/math/prim/fun/size_zero.hpp
+++ b/stan/math/prim/fun/size_zero.hpp
@@ -1,6 +1,7 @@
 #ifndef STAN_MATH_PRIM_FUN_SIZE_ZERO_HPP
 #define STAN_MATH_PRIM_FUN_SIZE_ZERO_HPP
 
+#include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <utility>
 
@@ -16,7 +17,7 @@ namespace math {
  */
 template <typename T>
 inline bool size_zero(const T& x) {
-  return !size(x);
+  return !math::size(x);
 }
 
 /**
diff --git a/stan/math/prim/functor/coupled_ode_system.hpp b/stan/math/prim/functor/coupled_ode_system.hpp
index a4dbe69603d..5298b53a570 100644
--- a/stan/math/prim/functor/coupled_ode_system.hpp
+++ b/stan/math/prim/functor/coupled_ode_system.hpp
@@ -67,9 +67,9 @@ struct coupled_ode_system_impl<true, F, T_y0, Args...> {
 
     dz_dt.resize(y.size());
 
-    Eigen::VectorXd f_y_t
-        = apply([&](const Args&... args) { return f_(t, y, msgs_, args...); },
-                args_tuple_);
+    Eigen::VectorXd f_y_t = math::apply(
+        [&](const Args&... args) { return f_(t, y, msgs_, args...); },
+        args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t.size(), "states",
                      y.size());
@@ -110,8 +110,8 @@ struct coupled_ode_system
                      const Eigen::Matrix<T_y0, Eigen::Dynamic, 1>& y0,
                      std::ostream* msgs, const Args&... args)
       : coupled_ode_system_impl<
-            std::is_arithmetic<return_type_t<T_y0, Args...>>::value, F, T_y0,
-            Args...>(f, y0, msgs, args...) {}
+          std::is_arithmetic<return_type_t<T_y0, Args...>>::value, F, T_y0,
+          Args...>(f, y0, msgs, args...) {}
 };
 
 }  // namespace math
diff --git a/stan/math/prim/functor/ode_rk45.hpp b/stan/math/prim/functor/ode_rk45.hpp
index 8a90bd1e2c7..4850b1169c6 100644
--- a/stan/math/prim/functor/ode_rk45.hpp
+++ b/stan/math/prim/functor/ode_rk45.hpp
@@ -79,7 +79,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   std::tuple<ref_type_t<Args>...> args_ref_tuple(args...);
 
-  apply(
+  math::apply(
       [&](const auto&... args_ref) {
         // Code from https://stackoverflow.com/a/17340003
         std::vector<int> unused_temp{
@@ -102,7 +102,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
 
   using return_t = return_type_t<T_y0, T_t0, T_ts, Args...>;
   // creates basic or coupled system by template specializations
-  auto&& coupled_system = apply(
+  auto&& coupled_system = math::apply(
       [&](const auto&... args_ref) {
         return coupled_ode_system<F, T_y0_t0, ref_type_t<Args>...>(f, y0, msgs,
                                                                    args_ref...);
@@ -128,7 +128,7 @@ ode_rk45_tol_impl(const char* function_name, const F& f, const T_y0& y0_arg,
       observer_initial_recorded = true;
       return;
     }
-    apply(
+    math::apply(
         [&](const auto&... args_ref) {
           y.emplace_back(ode_store_sensitivities(
               f, coupled_state, y0, t0, ts[time_index], msgs, args_ref...));
diff --git a/stan/math/prim/functor/reduce_sum.hpp b/stan/math/prim/functor/reduce_sum.hpp
index 868de060f32..ff0b4f5d592 100644
--- a/stan/math/prim/functor/reduce_sum.hpp
+++ b/stan/math/prim/functor/reduce_sum.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 
 #include <tbb/task_arena.h>
 #include <tbb/parallel_reduce.h>
@@ -80,7 +81,7 @@ struct reduce_sum_impl<ReduceFunction, require_arithmetic_t<ReturnType>,
         sub_slice.emplace_back(vmapped_[i]);
       }
 
-      sum_ += apply(
+      sum_ += math::apply(
           [&](auto&&... args) {
             return ReduceFunction()(sub_slice, r.begin(), r.end() - 1, &msgs_,
                                     args...);
diff --git a/stan/math/prim/mat/fun/LDLT_factor.hpp b/stan/math/prim/mat/fun/LDLT_factor.hpp
index ed6f1121c85..f782f289e9b 100644
--- a/stan/math/prim/mat/fun/LDLT_factor.hpp
+++ b/stan/math/prim/mat/fun/LDLT_factor.hpp
@@ -1 +1 @@
-#include "../prim/fun/LDLT_factor.hpp"
+#include "../../fun/LDLT_factor.hpp"
diff --git a/stan/math/prim/prob/bernoulli_lpmf.hpp b/stan/math/prim/prob/bernoulli_lpmf.hpp
index 4e5876f7359..c288a3ddd48 100644
--- a/stan/math/prim/prob/bernoulli_lpmf.hpp
+++ b/stan/math/prim/prob/bernoulli_lpmf.hpp
@@ -58,7 +58,7 @@ return_type_t<T_prob> bernoulli_lpmf(const T_n& n, const T_prob& theta) {
   scalar_seq_view<T_theta_ref> theta_vec(theta_ref);
   size_t N = max_size(n, theta);
 
-  if (size(theta) == 1) {
+  if (math::size(theta) == 1) {
     size_t sum = 0;
     for (size_t n = 0; n < N; n++) {
       sum += n_vec.val(n);
diff --git a/stan/math/prim/prob/beta_proportion_lpdf.hpp b/stan/math/prim/prob/beta_proportion_lpdf.hpp
index f3bf5ff409b..73aa67acab0 100644
--- a/stan/math/prim/prob/beta_proportion_lpdf.hpp
+++ b/stan/math/prim/prob/beta_proportion_lpdf.hpp
@@ -99,7 +99,7 @@ return_type_t<T_y, T_loc, T_prec> beta_proportion_lpdf(const T_y& y,
   size_t N = max_size(y, mu, kappa);
   T_partials_return logp(0);
   if (include_summand<propto, T_prec>::value) {
-    logp += sum(lgamma(kappa_val)) * N / size(kappa);
+    logp += sum(lgamma(kappa_val)) * N / math::size(kappa);
   }
   if (include_summand<propto, T_loc, T_prec>::value) {
     logp -= sum(lgamma(mukappa) + lgamma(kappa_val - mukappa)) * N
diff --git a/stan/math/prim/prob/binomial_logit_lpmf.hpp b/stan/math/prim/prob/binomial_logit_lpmf.hpp
index 3c4993b9495..816d06a75a7 100644
--- a/stan/math/prim/prob/binomial_logit_lpmf.hpp
+++ b/stan/math/prim/prob/binomial_logit_lpmf.hpp
@@ -92,10 +92,10 @@ return_type_t<T_prob> binomial_logit_lpmf(const T_n& n, const T_N& N,
       ops_partials.edge1_.partials_
           = n_val * inv_logit_neg_alpha - (N_val - n_val) * inv_logit_alpha;
     } else {
-      T_partials_return sum_n = sum(n_val) * maximum_size / size(n);
+      T_partials_return sum_n = sum(n_val) * maximum_size / math::size(n);
       ops_partials.edge1_.partials_[0] = forward_as<T_partials_return>(
           sum_n * inv_logit_neg_alpha
-          - (sum(N_val) * maximum_size / size(N) - sum_n) * inv_logit_alpha);
+          - (sum(N_val) * maximum_size / math::size(N) - sum_n) * inv_logit_alpha);
     }
   }
 
diff --git a/stan/math/prim/prob/cauchy_lpdf.hpp b/stan/math/prim/prob/cauchy_lpdf.hpp
index d1da16262b8..302c3c22a9b 100644
--- a/stan/math/prim/prob/cauchy_lpdf.hpp
+++ b/stan/math/prim/prob/cauchy_lpdf.hpp
@@ -92,7 +92,7 @@ return_type_t<T_y, T_loc, T_scale> cauchy_lpdf(const T_y& y, const T_loc& mu,
     logp -= N * LOG_PI;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/chi_square_cdf.hpp b/stan/math/prim/prob/chi_square_cdf.hpp
index ea439f86bf6..8e8f6cf5148 100644
--- a/stan/math/prim/prob/chi_square_cdf.hpp
+++ b/stan/math/prim/prob/chi_square_cdf.hpp
@@ -69,9 +69,9 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lccdf.hpp b/stan/math/prim/prob/chi_square_lccdf.hpp
index 8347406a373..cb6bead9799 100644
--- a/stan/math/prim/prob/chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/chi_square_lccdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> chi_square_lccdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lcdf.hpp b/stan/math/prim/prob/chi_square_lcdf.hpp
index dc7c196a454..5fefad3fe3a 100644
--- a/stan/math/prim/prob/chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/chi_square_lcdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> chi_square_lcdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/chi_square_lpdf.hpp b/stan/math/prim/prob/chi_square_lpdf.hpp
index 7bac96d7e19..f96260950ec 100644
--- a/stan/math/prim/prob/chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/chi_square_lpdf.hpp
@@ -81,12 +81,12 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
   T_partials_return logp(0);
   if (include_summand<propto, T_dof>::value) {
-    logp -= sum(nu_val * HALF_LOG_TWO + lgamma(half_nu)) * N / size(nu);
+    logp -= sum(nu_val * HALF_LOG_TWO + lgamma(half_nu)) * N / math::size(nu);
   }
   logp += sum((half_nu - 1.0) * log_y);
 
   if (include_summand<propto, T_y>::value) {
-    logp -= 0.5 * sum(y_val) * N / size(y);
+    logp -= 0.5 * sum(y_val) * N / math::size(y);
   }
 
   operands_and_partials<T_y_ref, T_nu_ref> ops_partials(y_ref, nu_ref);
diff --git a/stan/math/prim/prob/double_exponential_lpdf.hpp b/stan/math/prim/prob/double_exponential_lpdf.hpp
index 2587a1ed86f..832cafa50d8 100644
--- a/stan/math/prim/prob/double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/double_exponential_lpdf.hpp
@@ -91,7 +91,7 @@ return_type_t<T_y, T_loc, T_scale> double_exponential_lpdf(
     logp -= N * LOG_TWO;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   logp -= sum(scaled_diff);
 
diff --git a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
index d92de214ae5..0ee1f629f11 100644
--- a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
@@ -87,7 +87,7 @@ return_type_t<T_y, T_loc, T_scale, T_inv_scale> exp_mod_normal_lpdf(
     logp -= LOG_TWO * N;
   }
   if (include_summand<propto, T_inv_scale>::value) {
-    logp += sum(log(lambda_val)) * N / size(lambda);
+    logp += sum(log(lambda_val)) * N / math::size(lambda);
   }
   const auto& log_erfc_calc = log(erfc_calc);
   logp
diff --git a/stan/math/prim/prob/exponential_lccdf.hpp b/stan/math/prim/prob/exponential_lccdf.hpp
index f7d74bea979..f1920bbb2dc 100644
--- a/stan/math/prim/prob/exponential_lccdf.hpp
+++ b/stan/math/prim/prob/exponential_lccdf.hpp
@@ -51,7 +51,7 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     using beta_val_array = Eigen::Array<beta_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_y>::value && !is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = T_partials_array::Constant(
-          size(y), -forward_as<beta_val_scalar>(beta_val));
+          math::size(y), -forward_as<beta_val_scalar>(beta_val));
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
@@ -65,7 +65,7 @@ return_type_t<T_y, T_inv_scale> exponential_lccdf(const T_y& y,
     using y_val_array = Eigen::Array<y_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_inv_scale>::value && !is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = T_partials_array::Constant(
-          size(beta), -forward_as<y_val_scalar>(y_val));
+          math::size(beta), -forward_as<y_val_scalar>(y_val));
     } else if (is_vector<T_y>::value) {
       ops_partials.edge2_.partials_ = -forward_as<y_val_array>(y_val);
     } else {
diff --git a/stan/math/prim/prob/exponential_lpdf.hpp b/stan/math/prim/prob/exponential_lpdf.hpp
index 6f75b1a31ae..2cbb4e7e06a 100644
--- a/stan/math/prim/prob/exponential_lpdf.hpp
+++ b/stan/math/prim/prob/exponential_lpdf.hpp
@@ -92,7 +92,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
     using beta_val_array = Eigen::Array<beta_val_scalar, Eigen::Dynamic, 1>;
     if (is_vector<T_y>::value && !is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = T_partials_array::Constant(
-          size(y), -forward_as<beta_val_scalar>(beta_val));
+          math::size(y), -forward_as<beta_val_scalar>(beta_val));
     } else if (is_vector<T_inv_scale>::value) {
       ops_partials.edge1_.partials_ = -forward_as<beta_val_array>(beta_val);
     } else {
diff --git a/stan/math/prim/prob/frechet_lpdf.hpp b/stan/math/prim/prob/frechet_lpdf.hpp
index a6658972221..56bc37b2b5a 100644
--- a/stan/math/prim/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/prob/frechet_lpdf.hpp
@@ -79,7 +79,7 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, sigma);
   T_partials_return logp = -sum(sigma_div_y_pow_alpha);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp -= sum((alpha_val + 1.0) * log_y) * N / max_size(y, alpha);
diff --git a/stan/math/prim/prob/gamma_cdf.hpp b/stan/math/prim/prob/gamma_cdf.hpp
index 2918ebb9f22..fb9f88f788c 100644
--- a/stan/math/prim/prob/gamma_cdf.hpp
+++ b/stan/math/prim/prob/gamma_cdf.hpp
@@ -81,9 +81,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
   using std::pow;
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp
index 12fbd24a9f2..c031c1a7e63 100644
--- a/stan/math/prim/prob/gamma_lccdf.hpp
+++ b/stan/math/prim/prob/gamma_lccdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lcdf.hpp b/stan/math/prim/prob/gamma_lcdf.hpp
index b40c75bc7f3..19c15050873 100644
--- a/stan/math/prim/prob/gamma_lcdf.hpp
+++ b/stan/math/prim/prob/gamma_lcdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/gamma_lpdf.hpp b/stan/math/prim/prob/gamma_lpdf.hpp
index e23ed36ab8f..05aa0a13d3e 100644
--- a/stan/math/prim/prob/gamma_lpdf.hpp
+++ b/stan/math/prim/prob/gamma_lpdf.hpp
@@ -99,7 +99,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, beta);
   T_partials_return logp(0.0);
   if (include_summand<propto, T_shape>::value) {
-    logp = -sum(lgamma(alpha_val)) * N / size(alpha);
+    logp = -sum(lgamma(alpha_val)) * N / math::size(alpha);
   }
   const auto& log_y = to_ref_if<is_constant_all<T_shape>::value>(log(y_val));
   if (include_summand<propto, T_shape, T_inv_scale>::value) {
diff --git a/stan/math/prim/prob/gumbel_lpdf.hpp b/stan/math/prim/prob/gumbel_lpdf.hpp
index a0f97e7d111..ad02f4ac088 100644
--- a/stan/math/prim/prob/gumbel_lpdf.hpp
+++ b/stan/math/prim/prob/gumbel_lpdf.hpp
@@ -85,7 +85,7 @@ return_type_t<T_y, T_loc, T_scale> gumbel_lpdf(const T_y& y, const T_loc& mu,
   size_t N = max_size(y, mu, beta);
   T_partials_return logp = -sum(y_minus_mu_over_beta + exp_y_m_mu_over_beta);
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(beta_val)) * N / size(beta);
+    logp -= sum(log(beta_val)) * N / math::size(beta);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/inv_chi_square_cdf.hpp b/stan/math/prim/prob/inv_chi_square_cdf.hpp
index 1f1ff29c386..429dbfadcc7 100644
--- a/stan/math/prim/prob/inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_cdf.hpp
@@ -69,9 +69,9 @@ return_type_t<T_y, T_dof> inv_chi_square_cdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lccdf.hpp b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
index 394e7956f97..e7c6a1e9c88 100644
--- a/stan/math/prim/prob/inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> inv_chi_square_lccdf(const T_y& y, const T_dof& nu) {
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/inv_chi_square_lpdf.hpp b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
index 68ed6bfb274..1b41ceda34a 100644
--- a/stan/math/prim/prob/inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
@@ -87,11 +87,11 @@ return_type_t<T_y, T_dof> inv_chi_square_lpdf(const T_y& y, const T_dof& nu) {
   size_t N = max_size(y, nu);
   T_partials_return logp = -sum((half_nu + 1.0) * log_y);
   if (include_summand<propto, T_dof>::value) {
-    logp -= (sum(nu_val) * HALF_LOG_TWO + sum(lgamma(half_nu))) * N / size(nu);
+    logp -= (sum(nu_val) * HALF_LOG_TWO + sum(lgamma(half_nu))) * N / math::size(nu);
   }
   if (include_summand<propto, T_y>::value) {
     const auto& inv_y = to_ref_if<!is_constant_all<T_y>::value>(inv(y_val));
-    logp -= 0.5 * sum(inv_y) * N / size(y);
+    logp -= 0.5 * sum(inv_y) * N / math::size(y);
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_ = (0.5 * inv_y - half_nu - 1.0) * inv_y;
     }
diff --git a/stan/math/prim/prob/inv_gamma_cdf.hpp b/stan/math/prim/prob/inv_gamma_cdf.hpp
index 89271caa42d..3ac21f757d4 100644
--- a/stan/math/prim/prob/inv_gamma_cdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_cdf.hpp
@@ -80,9 +80,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_cdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lccdf.hpp b/stan/math/prim/prob/inv_gamma_lccdf.hpp
index 22212e5adff..c731eccd0e3 100644
--- a/stan/math/prim/prob/inv_gamma_lccdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lccdf.hpp
@@ -66,9 +66,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lccdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lcdf.hpp b/stan/math/prim/prob/inv_gamma_lcdf.hpp
index 27e6800b444..7ad55a35302 100644
--- a/stan/math/prim/prob/inv_gamma_lcdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lcdf.hpp
@@ -66,9 +66,9 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lcdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+      gamma_vec(math::size(alpha));
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+      digamma_vec(math::size(alpha));
 
   if (!is_constant_all<T_shape>::value) {
     for (size_t i = 0; i < stan::math::size(alpha); i++) {
diff --git a/stan/math/prim/prob/inv_gamma_lpdf.hpp b/stan/math/prim/prob/inv_gamma_lpdf.hpp
index 96af312d941..69f20d8aaf1 100644
--- a/stan/math/prim/prob/inv_gamma_lpdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lpdf.hpp
@@ -90,7 +90,7 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lpdf(const T_y& y,
 
   size_t N = max_size(y, alpha, beta);
   if (include_summand<propto, T_shape>::value) {
-    logp -= sum(lgamma(alpha_val)) * N / size(alpha);
+    logp -= sum(lgamma(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_shape, T_scale>::value) {
     const auto& log_beta
diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp
index 4016e423324..46854bd3be6 100644
--- a/stan/math/prim/prob/logistic_lpdf.hpp
+++ b/stan/math/prim/prob/logistic_lpdf.hpp
@@ -73,7 +73,7 @@ return_type_t<T_y, T_loc, T_scale> logistic_lpdf(const T_y& y, const T_loc& mu,
   T_partials_return logp = -sum(y_minus_mu_div_sigma)
                            - 2.0 * sum(log1p(exp(-y_minus_mu_div_sigma)));
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_scale>::value) {
diff --git a/stan/math/prim/prob/lognormal_lpdf.hpp b/stan/math/prim/prob/lognormal_lpdf.hpp
index edb61e3688d..facf5e3e1e5 100644
--- a/stan/math/prim/prob/lognormal_lpdf.hpp
+++ b/stan/math/prim/prob/lognormal_lpdf.hpp
@@ -81,10 +81,10 @@ return_type_t<T_y, T_loc, T_scale> lognormal_lpdf(const T_y& y, const T_loc& mu,
   T_partials_return logp
       = N * NEG_LOG_SQRT_TWO_PI - 0.5 * sum(square(logy_m_mu) * inv_sigma_sq);
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y>::value) {
-    logp -= sum(log_y) * N / size(y);
+    logp -= sum(log_y) * N / math::size(y);
   }
 
   if (!is_constant_all<T_y, T_loc, T_scale>::value) {
diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp
index 31ca456c4ec..f0b92624450 100644
--- a/stan/math/prim/prob/normal_lpdf.hpp
+++ b/stan/math/prim/prob/normal_lpdf.hpp
@@ -92,7 +92,7 @@ inline return_type_t<T_y, T_loc, T_scale> normal_lpdf(const T_y& y,
     logp += NEG_LOG_SQRT_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_scale, T_loc>::value) {
diff --git a/stan/math/prim/prob/normal_sufficient_lpdf.hpp b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
index f4087fdd158..f741d72c11d 100644
--- a/stan/math/prim/prob/normal_sufficient_lpdf.hpp
+++ b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
@@ -117,7 +117,7 @@ return_type_t<T_y, T_s, T_loc, T_scale> normal_sufficient_lpdf(
   size_t N = max_size(y_bar, s_squared, n_obs, mu, sigma);
   T_partials_return logp = -sum(cons_expr / (2 * sigma_squared));
   if (include_summand<propto>::value) {
-    logp += NEG_LOG_SQRT_TWO_PI * sum(n_obs_val) * N / size(n_obs);
+    logp += NEG_LOG_SQRT_TWO_PI * sum(n_obs_val) * N / math::size(n_obs);
   }
   if (include_summand<propto, T_scale>::value) {
     logp -= sum(n_obs_val * log(sigma_val)) * N / max_size(n_obs, sigma);
@@ -151,7 +151,7 @@ return_type_t<T_y, T_s, T_loc, T_scale> normal_sufficient_lpdf(
       } else {
         forward_as<internal::broadcast_array<T_partials_return>>(
             ops_partials.edge2_.partials_)
-            = -0.5 / sigma_squared * N / size(sigma);
+            = -0.5 / sigma_squared * N / math::size(sigma);
       }
     }
   }
diff --git a/stan/math/prim/prob/ordered_logistic_lpmf.hpp b/stan/math/prim/prob/ordered_logistic_lpmf.hpp
index b9a90dabb1f..35516fcffcf 100644
--- a/stan/math/prim/prob/ordered_logistic_lpmf.hpp
+++ b/stan/math/prim/prob/ordered_logistic_lpmf.hpp
@@ -86,7 +86,7 @@ return_type_t<T_loc, T_cut> ordered_logistic_lpmf(const T_y& y,
   T_cut_ref c_ref = c;
   vector_seq_view<T_cut_ref> c_vec(c_ref);
   int K = c_vec[0].size() + 1;
-  int N = size(lambda);
+  int N = math::size(lambda);
   int C_l = size_mvt(c);
 
   check_consistent_sizes(function, "Integers", y, "Locations", lambda);
diff --git a/stan/math/prim/prob/ordered_probit_lpmf.hpp b/stan/math/prim/prob/ordered_probit_lpmf.hpp
index b7a89401737..18bbfdd0174 100644
--- a/stan/math/prim/prob/ordered_probit_lpmf.hpp
+++ b/stan/math/prim/prob/ordered_probit_lpmf.hpp
@@ -54,7 +54,7 @@ return_type_t<T_loc, T_cut> ordered_probit_lpmf(const T_y& y,
   static const char* function = "ordered_probit";
 
   check_nonzero_size(function, "Cut-points", c);
-  int N = size(lambda);
+  int N = math::size(lambda);
   int C_l = size_mvt(c);
   vector_seq_view<T_cut> c_vec(c);
   int K = c_vec[0].size() + 1;
diff --git a/stan/math/prim/prob/pareto_lpdf.hpp b/stan/math/prim/prob/pareto_lpdf.hpp
index 9ef1b53ded0..c39b86f6a14 100644
--- a/stan/math/prim/prob/pareto_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_lpdf.hpp
@@ -70,7 +70,7 @@ return_type_t<T_y, T_scale, T_shape> pareto_lpdf(const T_y& y,
   size_t N = max_size(y, y_min, alpha);
   T_partials_return logp(0);
   if (include_summand<propto, T_shape>::value) {
-    logp = sum(log(alpha_val)) * N / size(alpha);
+    logp = sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp -= sum(alpha_val * log_y + log_y) * N / max_size(alpha, y);
diff --git a/stan/math/prim/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
index efb213cc3c3..b6beef26c1d 100644
--- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
@@ -74,10 +74,10 @@ return_type_t<T_y, T_loc, T_scale, T_shape> pareto_type_2_lpdf(
   size_t N = max_size(y, mu, lambda, alpha);
   T_partials_return logp(0.0);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(lambda_val)) * N / size(lambda);
+    logp -= sum(log(lambda_val)) * N / math::size(lambda);
   }
   if (include_summand<propto, T_y, T_loc, T_scale, T_shape>::value) {
     logp -= sum((alpha_val + 1.0) * log1p_scaled_diff);
diff --git a/stan/math/prim/prob/poisson_log_lpmf.hpp b/stan/math/prim/prob/poisson_log_lpmf.hpp
index 0568c1dcc89..1f927685252 100644
--- a/stan/math/prim/prob/poisson_log_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_log_lpmf.hpp
@@ -78,10 +78,10 @@ return_type_t<T_log_rate> poisson_log_lpmf(const T_n& n,
 
   T_partials_return logp = sum(n_val * alpha_val);
   if (include_summand<propto, T_log_rate>::value) {
-    logp -= sum(exp_alpha) * N / size(alpha);
+    logp -= sum(exp_alpha) * N / math::size(alpha);
   }
   if (include_summand<propto>::value) {
-    logp -= sum(lgamma(n_val + 1.0)) * N / size(n);
+    logp -= sum(lgamma(n_val + 1.0)) * N / math::size(n);
   }
 
   if (!is_constant_all<T_log_rate>::value) {
diff --git a/stan/math/prim/prob/poisson_lpmf.hpp b/stan/math/prim/prob/poisson_lpmf.hpp
index f9bd2320a81..bf11e44ef0c 100644
--- a/stan/math/prim/prob/poisson_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_lpmf.hpp
@@ -72,10 +72,10 @@ return_type_t<T_rate> poisson_lpmf(const T_n& n, const T_rate& lambda) {
 
   T_partials_return logp = stan::math::sum(multiply_log(n_val, lambda_val));
   if (include_summand<propto, T_rate>::value) {
-    logp -= sum(lambda_val) * N / size(lambda);
+    logp -= sum(lambda_val) * N / math::size(lambda);
   }
   if (include_summand<propto>::value) {
-    logp -= sum(lgamma(n_val + 1.0)) * N / size(n);
+    logp -= sum(lgamma(n_val + 1.0)) * N / math::size(n);
   }
 
   if (!is_constant_all<T_rate>::value) {
diff --git a/stan/math/prim/prob/rayleigh_lpdf.hpp b/stan/math/prim/prob/rayleigh_lpdf.hpp
index 7582d3dd7ca..7d43a4a12a9 100644
--- a/stan/math/prim/prob/rayleigh_lpdf.hpp
+++ b/stan/math/prim/prob/rayleigh_lpdf.hpp
@@ -61,10 +61,10 @@ return_type_t<T_y, T_scale> rayleigh_lpdf(const T_y& y, const T_scale& sigma) {
   size_t N = max_size(y, sigma);
   T_partials_return logp = -0.5 * sum(square(y_over_sigma));
   if (include_summand<propto, T_scale>::value) {
-    logp -= 2.0 * sum(log(sigma_val)) * N / size(sigma);
+    logp -= 2.0 * sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y>::value) {
-    logp += sum(log(y_val)) * N / size(y);
+    logp += sum(log(y_val)) * N / math::size(y);
   }
 
   if (!is_constant_all<T_y, T_scale>::value) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
index f441852d649..7ec39b46e56 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
@@ -78,9 +78,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_cdf(const T_y& y,
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
index 83b31549806..52c122b6f4a 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lccdf(
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
index c25413b2831..a5f480a438c 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
@@ -64,9 +64,9 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lcdf(
   }
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
+      gamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     for (size_t i = 0; i < stan::math::size(nu); i++) {
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
index 5187e4b03f5..8c03961ed3f 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
@@ -86,7 +86,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
                 T_partials_return, T_dof>
-      half_nu(size(nu));
+      half_nu(math::size(nu));
   for (size_t i = 0; i < stan::math::size(nu); i++) {
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
       half_nu[i] = 0.5 * nu_vec.val(i);
@@ -95,7 +95,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
+      log_y(math::size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (include_summand<propto, T_dof, T_y>::value) {
       log_y[i] = log(y_vec.val(i));
@@ -104,7 +104,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
                 T_partials_return, T_y>
-      inv_y(size(y));
+      inv_y(math::size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
       inv_y[i] = 1.0 / y_vec.val(i);
@@ -113,7 +113,7 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
 
   VectorBuilder<include_summand<propto, T_dof, T_scale>::value,
                 T_partials_return, T_scale>
-      log_s(size(s));
+      log_s(math::size(s));
   for (size_t i = 0; i < stan::math::size(s); i++) {
     if (include_summand<propto, T_dof, T_scale>::value) {
       log_s[i] = log(s_vec.val(i));
@@ -121,11 +121,11 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
   }
 
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      log_half_nu(size(nu));
+      log_half_nu(math::size(nu));
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
+      digamma_half_nu_over_two(math::size(nu));
   for (size_t i = 0; i < stan::math::size(nu); i++) {
     if (include_summand<propto, T_dof>::value) {
       lgamma_half_nu[i] = lgamma(half_nu[i]);
diff --git a/stan/math/prim/prob/skew_normal_lpdf.hpp b/stan/math/prim/prob/skew_normal_lpdf.hpp
index e7de3a48dd0..92d89406c48 100644
--- a/stan/math/prim/prob/skew_normal_lpdf.hpp
+++ b/stan/math/prim/prob/skew_normal_lpdf.hpp
@@ -85,7 +85,7 @@ return_type_t<T_y, T_loc, T_scale, T_shape> skew_normal_lpdf(
     logp -= HALF_LOG_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
   if (include_summand<propto, T_y, T_loc, T_scale>::value) {
     logp -= sum(square(y_minus_mu_over_sigma)) * 0.5 * N
diff --git a/stan/math/prim/prob/student_t_cdf.hpp b/stan/math/prim/prob/student_t_cdf.hpp
index 38232b2ea0a..87c50cccf44 100644
--- a/stan/math/prim/prob/student_t_cdf.hpp
+++ b/stan/math/prim/prob/student_t_cdf.hpp
@@ -65,11 +65,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_cdf(const T_y& y,
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lccdf.hpp b/stan/math/prim/prob/student_t_lccdf.hpp
index 5c9e39865d2..cfc6e5608ef 100644
--- a/stan/math/prim/prob/student_t_lccdf.hpp
+++ b/stan/math/prim/prob/student_t_lccdf.hpp
@@ -65,11 +65,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lccdf(
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lcdf.hpp b/stan/math/prim/prob/student_t_lcdf.hpp
index 401e938cc5b..84158993b47 100644
--- a/stan/math/prim/prob/student_t_lcdf.hpp
+++ b/stan/math/prim/prob/student_t_lcdf.hpp
@@ -67,11 +67,11 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lcdf(const T_y& y,
   T_partials_return digammaHalf = 0;
 
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+      digamma_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNu_vec(size(nu));
+      digammaNu_vec(math::size(nu));
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digammaNuPlusHalf_vec(size(nu));
+      digammaNuPlusHalf_vec(math::size(nu));
 
   if (!is_constant_all<T_dof>::value) {
     digammaHalf = digamma(0.5);
diff --git a/stan/math/prim/prob/student_t_lpdf.hpp b/stan/math/prim/prob/student_t_lpdf.hpp
index 3c39c228a42..7f203582914 100644
--- a/stan/math/prim/prob/student_t_lpdf.hpp
+++ b/stan/math/prim/prob/student_t_lpdf.hpp
@@ -120,10 +120,10 @@ return_type_t<T_y, T_dof, T_loc, T_scale> student_t_lpdf(const T_y& y,
   if (include_summand<propto, T_dof>::value) {
     logp += (sum(lgamma(half_nu + 0.5)) - sum(lgamma(half_nu))
              - 0.5 * sum(log(nu_val)))
-            * N / size(nu);
+            * N / math::size(nu);
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log(sigma_val)) * N / size(sigma);
+    logp -= sum(log(sigma_val)) * N / math::size(sigma);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/uniform_lcdf.hpp b/stan/math/prim/prob/uniform_lcdf.hpp
index 0962208e046..67e9876e3ee 100644
--- a/stan/math/prim/prob/uniform_lcdf.hpp
+++ b/stan/math/prim/prob/uniform_lcdf.hpp
@@ -73,7 +73,7 @@ return_type_t<T_y, T_low, T_high> uniform_lcdf(const T_y& y, const T_low& alpha,
   if (!is_constant_all<T_y>::value) {
     if (!is_vector<T_y>::value && is_vector<T_high>::value
         && !is_vector<T_low>::value) {
-      ops_partials.edge1_.partials_ = size(beta) * inv(y_minus_alpha);
+      ops_partials.edge1_.partials_ = math::size(beta) * inv(y_minus_alpha);
     } else {
       ops_partials.edge1_.partials_ = inv(y_minus_alpha);
     }
@@ -85,7 +85,7 @@ return_type_t<T_y, T_low, T_high> uniform_lcdf(const T_y& y, const T_low& alpha,
   if (!is_constant_all<T_high>::value) {
     if (is_vector<T_y>::value && !is_vector<T_low>::value
         && !is_vector<T_high>::value) {
-      ops_partials.edge3_.partials_ = inv(-b_minus_a) * size(y);
+      ops_partials.edge3_.partials_ = inv(-b_minus_a) * math::size(y);
     } else {
       ops_partials.edge3_.partials_ = inv(-b_minus_a);
     }
diff --git a/stan/math/prim/prob/uniform_lpdf.hpp b/stan/math/prim/prob/uniform_lpdf.hpp
index 73c875a17a3..837b912b4e7 100644
--- a/stan/math/prim/prob/uniform_lpdf.hpp
+++ b/stan/math/prim/prob/uniform_lpdf.hpp
@@ -102,7 +102,7 @@ return_type_t<T_y, T_low, T_high> uniform_lpdf(const T_y& y, const T_low& alpha,
     if (!is_constant_all<T_high>::value) {
       if (is_vector<T_y>::value && !is_vector<T_low>::value
           && !is_vector<T_high>::value) {
-        ops_partials.edge3_.partials_ = -inv_beta_minus_alpha * size(y);
+        ops_partials.edge3_.partials_ = -inv_beta_minus_alpha * math::size(y);
       } else {
         ops_partials.edge3_.partials_ = -inv_beta_minus_alpha;
       }
@@ -110,7 +110,7 @@ return_type_t<T_y, T_low, T_high> uniform_lpdf(const T_y& y, const T_low& alpha,
     if (!is_constant_all<T_low>::value) {
       if (is_vector<T_y>::value && !is_vector<T_low>::value
           && !is_vector<T_high>::value) {
-        ops_partials.edge2_.partials_ = inv_beta_minus_alpha * size(y);
+        ops_partials.edge2_.partials_ = inv_beta_minus_alpha * math::size(y);
       } else {
         ops_partials.edge2_.partials_ = std::move(inv_beta_minus_alpha);
       }
diff --git a/stan/math/prim/prob/von_mises_lpdf.hpp b/stan/math/prim/prob/von_mises_lpdf.hpp
index 8fe6650ce2c..910eef74347 100644
--- a/stan/math/prim/prob/von_mises_lpdf.hpp
+++ b/stan/math/prim/prob/von_mises_lpdf.hpp
@@ -70,7 +70,7 @@ return_type_t<T_y, T_loc, T_scale> von_mises_lpdf(T_y const& y, T_loc const& mu,
     logp -= LOG_TWO_PI * N;
   }
   if (include_summand<propto, T_scale>::value) {
-    logp -= sum(log_modified_bessel_first_kind(0, kappa_val)) * N / size(kappa);
+    logp -= sum(log_modified_bessel_first_kind(0, kappa_val)) * N / math::size(kappa);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/weibull_lpdf.hpp b/stan/math/prim/prob/weibull_lpdf.hpp
index e355beeac2f..72e4bbc0294 100644
--- a/stan/math/prim/prob/weibull_lpdf.hpp
+++ b/stan/math/prim/prob/weibull_lpdf.hpp
@@ -96,7 +96,7 @@ return_type_t<T_y, T_shape, T_scale> weibull_lpdf(const T_y& y,
   size_t N = max_size(y, alpha, sigma);
   T_partials_return logp = -sum(y_div_sigma_pow_alpha);
   if (include_summand<propto, T_shape>::value) {
-    logp += sum(log(alpha_val)) * N / size(alpha);
+    logp += sum(log(alpha_val)) * N / math::size(alpha);
   }
   if (include_summand<propto, T_y, T_shape>::value) {
     logp += sum((alpha_val - 1.0) * log_y) * N / max_size(alpha, y);
diff --git a/stan/math/rev/functor/adj_jac_apply.hpp b/stan/math/rev/functor/adj_jac_apply.hpp
index 66826b7a7b9..5d81bc6ad61 100644
--- a/stan/math/rev/functor/adj_jac_apply.hpp
+++ b/stan/math/rev/functor/adj_jac_apply.hpp
@@ -509,7 +509,7 @@ struct adj_jac_vari : public vari {
     internal::build_y_adj(y_vi_, M_, y_adj);
     auto y_adj_jacs = f_.multiply_adjoint_jacobian(is_var_, y_adj);
 
-    apply([this](auto&&... args) { this->accumulate_adjoints(args...); },
+    math::apply([this](auto&&... args) { this->accumulate_adjoints(args...); },
           y_adj_jacs);
   }
 };
diff --git a/stan/math/rev/functor/coupled_ode_system.hpp b/stan/math/rev/functor/coupled_ode_system.hpp
index 83783d4d03f..33c0596a0a7 100644
--- a/stan/math/rev/functor/coupled_ode_system.hpp
+++ b/stan/math/rev/functor/coupled_ode_system.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/functor/coupled_ode_system.hpp>
 #include <stan/math/rev/functor/cvodes_utils.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
@@ -125,7 +126,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
       y_vars.coeffRef(n) = z[n];
 
     Eigen::Matrix<var, Eigen::Dynamic, 1> f_y_t_vars
-        = apply([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
+        = math::apply([&](auto&&... args) { return f_(t, y_vars, msgs_, args...); },
                 local_args_tuple_);
 
     check_size_match("coupled_ode_system", "dy_dt", f_y_t_vars.size(), "states",
@@ -140,7 +141,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
       // memset was faster than Eigen setZero
       memset(args_adjoints_.data(), 0, sizeof(double) * num_args_vars);
 
-      apply(
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },
@@ -148,7 +149,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
 
       // The vars here do not live on the nested stack so must be zero'd
       // separately
-      apply([&](auto&&... args) { zero_adjoints(args...); }, local_args_tuple_);
+      math::apply([&](auto&&... args) { zero_adjoints(args...); }, local_args_tuple_);
 
       // No need to zero adjoints after last sweep
       if (i + 1 < N_) {
diff --git a/stan/math/rev/functor/cvodes_integrator.hpp b/stan/math/rev/functor/cvodes_integrator.hpp
index 4e40354540b..83c395160c5 100644
--- a/stan/math/rev/functor/cvodes_integrator.hpp
+++ b/stan/math/rev/functor/cvodes_integrator.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/rev/functor/cvodes_utils.hpp>
 #include <stan/math/rev/functor/coupled_ode_system.hpp>
 #include <stan/math/rev/functor/ode_store_sensitivities.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <cvodes/cvodes.h>
@@ -103,7 +104,7 @@ class cvodes_integrator {
     const Eigen::VectorXd y_vec = Eigen::Map<const Eigen::VectorXd>(y, N_);
 
     Eigen::VectorXd dy_dt_vec
-        = apply([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
+        = math::apply([&](auto&&... args) { return f_(t, y_vec, msgs_, args...); },
                 value_of_args_tuple_);
 
     check_size_match("cvodes_integrator", "dy_dt", dy_dt_vec.size(), "states",
@@ -121,7 +122,7 @@ class cvodes_integrator {
     Eigen::MatrixXd Jfy;
 
     auto f_wrapped = [&](const Eigen::Matrix<var, Eigen::Dynamic, 1>& y) {
-      return apply([&](auto&&... args) { return f_(t, y, msgs_, args...); },
+      return math::apply([&](auto&&... args) { return f_(t, y, msgs_, args...); },
                    value_of_args_tuple_);
     };
 
@@ -211,7 +212,7 @@ class cvodes_integrator {
 
     // Code from: https://stackoverflow.com/a/17340003 . Should probably do
     // something better
-    apply(
+    math::apply(
         [&](auto&&... args) {
           std::vector<int> unused_temp{
               0, (check_finite(function_name, "ode parameters and data", args),
@@ -326,7 +327,7 @@ class cvodes_integrator {
           }
         }
 
-        y.emplace_back(apply(
+        y.emplace_back(math::apply(
             [&](auto&&... args) {
               return ode_store_sensitivities(f_, coupled_state_, y0_, t0_,
                                              ts_[n], msgs_, args...);
diff --git a/stan/math/rev/functor/ode_adams.hpp b/stan/math/rev/functor/ode_adams.hpp
index ee0bdafbbd5..500cd2e4c92 100644
--- a/stan/math/rev/functor/ode_adams.hpp
+++ b/stan/math/rev/functor/ode_adams.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/cvodes_integrator.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <ostream>
 #include <vector>
 
@@ -54,7 +55,7 @@ ode_adams_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                    long int max_num_steps,  // NOLINT(runtime/int)
                    std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return apply(
+  return math::apply(
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_ADAMS, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/ode_bdf.hpp b/stan/math/rev/functor/ode_bdf.hpp
index a07af2e3339..8042f848b69 100644
--- a/stan/math/rev/functor/ode_bdf.hpp
+++ b/stan/math/rev/functor/ode_bdf.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/cvodes_integrator.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <ostream>
 #include <vector>
@@ -55,7 +56,7 @@ ode_bdf_tol_impl(const char* function_name, const F& f, const T_y0& y0,
                  long int max_num_steps,  // NOLINT(runtime/int)
                  std::ostream* msgs, const T_Args&... args) {
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  return apply(
+  return math::apply(
       [&](const auto&... args_refs) {
         cvodes_integrator<CV_BDF, F, T_y0, T_t0, T_ts, ref_type_t<T_Args>...>
         integrator(function_name, f, y0, t0, ts, relative_tolerance,
diff --git a/stan/math/rev/functor/reduce_sum.hpp b/stan/math/rev/functor/reduce_sum.hpp
index dfe51dee2a7..2ee3927618a 100644
--- a/stan/math/rev/functor/reduce_sum.hpp
+++ b/stan/math/rev/functor/reduce_sum.hpp
@@ -113,7 +113,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
         // scope. In this case no need for zeroing adjoints, since the
         // fresh copy has all adjoints set to zero.
         local_args_tuple_scope_.stack_.execute([&]() {
-          apply(
+          math::apply(
               [&](auto&&... args) {
                 local_args_tuple_scope_.args_tuple_holder_ = std::make_unique<
                     typename scoped_args_tuple::args_tuple_t>(
@@ -140,7 +140,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
       }
 
       // Perform calculation
-      var sub_sum_v = apply(
+      var sub_sum_v = math::apply(
           [&](auto&&... args) {
             return ReduceFunction()(local_sub_slice, r.begin(), r.end() - 1,
                                     &msgs_, args...);
@@ -158,7 +158,7 @@ struct reduce_sum_impl<ReduceFunction, require_var_t<ReturnType>, ReturnType,
                           std::move(local_sub_slice));
 
       // Accumulate adjoints of shared_arguments
-      apply(
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },