diff --git a/Jenkinsfile b/Jenkinsfile
index 890677505fe..f3d438552fd 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -87,7 +87,7 @@ pipeline {
     }
     environment {
         STAN_NUM_THREADS = 4
-        CLANG_CXX = 'clang++-6.0'
+        CLANG_CXX = 'clang++-7'
         GCC = 'g++'
         MPICXX = 'mpicxx.openmpi'
         N_TESTS = 150
@@ -116,7 +116,7 @@ pipeline {
         stage("Clang-format") {
             agent {
                 docker {
-                    image 'stanorg/ci:gpu'
+                    image 'stanorg/ci:gpu-cpp17'
                     label 'linux'
                 }
             }
@@ -167,7 +167,7 @@ pipeline {
         stage('Linting & Doc checks') {
             agent {
                 docker {
-                    image 'stanorg/ci:gpu'
+                    image 'stanorg/ci:gpu-cpp17'
                     label 'linux'
                 }
             }
@@ -202,7 +202,7 @@ pipeline {
         stage('Verify changes') {
             agent {
                 docker {
-                    image 'stanorg/ci:gpu'
+                    image 'stanorg/ci:gpu-cpp17'
                     label 'linux'
                 }
             }
@@ -224,7 +224,7 @@ pipeline {
         stage('Headers check') {
             agent {
                 docker {
-                    image 'stanorg/ci:gpu'
+                    image 'stanorg/ci:gpu-cpp17'
                     label 'linux'
                 }
             }
@@ -252,7 +252,7 @@ pipeline {
                 stage('Rev/Fwd Unit Tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
                             args '--cap-add SYS_PTRACE'
                         }
@@ -275,7 +275,7 @@ pipeline {
                 stage('Mix Unit Tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
                             args '--cap-add SYS_PTRACE'
                         }
@@ -297,7 +297,7 @@ pipeline {
                 stage('Prim Unit Tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
                             args '--cap-add SYS_PTRACE'
                         }
@@ -313,7 +313,7 @@ pipeline {
                         script {
                             runTests("test/unit/*_test.cpp", false)
                             runTests("test/unit/math/*_test.cpp", false)
-                            runTests("test/unit/math/prim", false)                            
+                            runTests("test/unit/math/prim", false)
                             runTests("test/unit/math/memory", false)
                         }
                     }
@@ -332,7 +332,7 @@ pipeline {
                 stage('MPI tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
                         }
                     }
@@ -352,7 +352,7 @@ pipeline {
                 stage('OpenCL GPU tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'v100'
                             args '--gpus 1'
                         }
@@ -375,9 +375,8 @@ pipeline {
                 stage('Distribution tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
-                            args '--pull always'
                         }
                     }
                     steps {
@@ -409,9 +408,8 @@ pipeline {
                 stage('Expressions test') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
-                            args '--pull always'
                         }
                     }
                     steps {
@@ -445,9 +443,8 @@ pipeline {
                 stage('Threading tests') {
                     agent {
                         docker {
-                            image 'stanorg/ci:gpu'
+                            image 'stanorg/ci:gpu-cpp17'
                             label 'linux'
-                            args '--pull always'
                         }
                     }
                     steps {
@@ -495,7 +492,7 @@ pipeline {
         stage('Upload doxygen') {
             agent {
                 docker {
-                    image 'stanorg/ci:gpu'
+                    image 'stanorg/ci:gpu-cpp17'
                     label 'linux'
                 }
             }
diff --git a/stan/math/opencl/matrix_cl.hpp b/stan/math/opencl/matrix_cl.hpp
index 6f92f6f62b6..6c5470e055a 100644
--- a/stan/math/opencl/matrix_cl.hpp
+++ b/stan/math/opencl/matrix_cl.hpp
@@ -283,7 +283,7 @@ class matrix_cl : 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();
@@ -321,7 +321,7 @@ class matrix_cl : 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<
@@ -457,7 +457,7 @@ class matrix_cl : public matrix_cl_base {
       return *this;
     }
     this->wait_for_read_write_events();
-    if (size() != a.size()) {
+    if (this->size() != a.size()) {
       buffer_cl_ = cl::Buffer(opencl_context.context(), CL_MEM_READ_WRITE,
                               sizeof(T) * a.size());
     }
@@ -518,7 +518,7 @@ class matrix_cl : 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();
@@ -538,7 +538,7 @@ class matrix_cl : 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();
@@ -577,7 +577,7 @@ class matrix_cl : 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());
@@ -587,7 +587,7 @@ class matrix_cl : 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_cdf.hpp b/stan/math/opencl/prim/bernoulli_cdf.hpp
index 219e5785050..38e987f3b4d 100644
--- a/stan/math/opencl/prim/bernoulli_cdf.hpp
+++ b/stan/math/opencl/prim/bernoulli_cdf.hpp
@@ -37,7 +37,7 @@ return_type_t<T_prob_cl> bernoulli_cdf(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 1.0;
   }
diff --git a/stan/math/opencl/prim/bernoulli_lccdf.hpp b/stan/math/opencl/prim/bernoulli_lccdf.hpp
index 44530cbde89..d2fd582b775 100644
--- a/stan/math/opencl/prim/bernoulli_lccdf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lccdf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lccdf(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_lcdf.hpp b/stan/math/opencl/prim/bernoulli_lcdf.hpp
index 2b3dafc778e..30e42538f43 100644
--- a/stan/math/opencl/prim/bernoulli_lcdf.hpp
+++ b/stan/math/opencl/prim/bernoulli_lcdf.hpp
@@ -38,7 +38,7 @@ return_type_t<T_prob_cl> bernoulli_lcdf(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_logit_glm_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
index 038e024e7b1..9de5422b293 100644
--- a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp
@@ -66,13 +66,14 @@ 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 337daf326dd..baf4974f60d 100644
--- a/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/categorical_logit_glm_lpmf.hpp
@@ -61,10 +61,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 4f0ac0b2aec..1bfec8258b6 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,18 @@ 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 b2210f94683..b999c25e535 100644
--- a/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
+++ b/stan/math/opencl/prim/normal_id_glm_lpdf.hpp
@@ -74,17 +74,18 @@ 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 b09e48c990c..5f5912b2328 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/ordered_logistic_lpmf.hpp b/stan/math/opencl/prim/ordered_logistic_lpmf.hpp
index e72276b05aa..9e248340335 100644
--- a/stan/math/opencl/prim/ordered_logistic_lpmf.hpp
+++ b/stan/math/opencl/prim/ordered_logistic_lpmf.hpp
@@ -73,8 +73,8 @@ inline return_type_t<T_y_cl, T_loc_cl, T_cuts_cl> ordered_logistic_lpmf(
   static const char* function = "ordered_logistic_lpmf(OpenCL)";
 
   if (size(y) != 1) {
-    check_size_match(function, "Size of ", "y", size(y), "Size of", "lambda",
-                     size(lambda));
+    check_size_match(function, "Size of ", "y", math::size(y), "Size of",
+                     "lambda", math::size(lambda));
   }
 
   int N_instances = max_size(y, lambda);
diff --git a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
index 541fe38d5be..c6b730e7644 100644
--- a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
+++ b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp
@@ -66,13 +66,14 @@ 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 4d01d570441..31ed3dfaecd 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 6e7ef9dd2aa..2878e9d2849 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_cdf.hpp b/stan/math/opencl/prim/std_normal_cdf.hpp
index 23448b45747..739d496780b 100644
--- a/stan/math/opencl/prim/std_normal_cdf.hpp
+++ b/stan/math/opencl/prim/std_normal_cdf.hpp
@@ -29,7 +29,7 @@ return_type_t<T_y_cl> std_normal_cdf(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 1.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_lccdf.hpp b/stan/math/opencl/prim/std_normal_lccdf.hpp
index c16138d944c..1e531fdaf0f 100644
--- a/stan/math/opencl/prim/std_normal_lccdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lccdf.hpp
@@ -30,7 +30,7 @@ return_type_t<T_y_cl> std_normal_lccdf(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 1.0;
   }
diff --git a/stan/math/opencl/prim/std_normal_lcdf.hpp b/stan/math/opencl/prim/std_normal_lcdf.hpp
index ffee5704877..86966cf7a60 100644
--- a/stan/math/opencl/prim/std_normal_lcdf.hpp
+++ b/stan/math/opencl/prim/std_normal_lcdf.hpp
@@ -183,7 +183,7 @@ return_type_t<T_y_cl> std_normal_lcdf(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 1.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 1c916c0e90c..3ff0b92842a 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>::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 a1fd18ee228..85a58ce4771 100644
--- a/stan/math/prim/fun/offset_multiplier_constrain.hpp
+++ b/stan/math/prim/fun/offset_multiplier_constrain.hpp
@@ -9,6 +9,7 @@
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/sum.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <cmath>
 
 namespace stan {
@@ -106,8 +107,8 @@ inline auto offset_multiplier_constrain(const T& x, const M& mu, const S& sigma,
   check_finite("offset_multiplier_constrain", "offset", value_of_rec(mu_ref));
   check_positive_finite("offset_multiplier_constrain", "multiplier",
                         value_of_rec(sigma_ref));
-  if (size(sigma_ref) == 1) {
-    lp += sum(multiply_log(size(x), sigma_ref));
+  if (math::size(sigma_ref) == 1) {
+    lp += sum(multiply_log(math::size(x), sigma_ref));
   } else {
     lp += sum(log(sigma_ref));
   }
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 90df071de8b..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());
diff --git a/stan/math/prim/functor/ode_ckrk.hpp b/stan/math/prim/functor/ode_ckrk.hpp
index 3f8fd0a0972..10c7f901aa7 100644
--- a/stan/math/prim/functor/ode_ckrk.hpp
+++ b/stan/math/prim/functor/ode_ckrk.hpp
@@ -78,7 +78,7 @@ ode_ckrk_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{
@@ -101,7 +101,7 @@ ode_ckrk_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...);
@@ -127,7 +127,7 @@ ode_ckrk_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/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/prob/bernoulli_lpmf.hpp b/stan/math/prim/prob/bernoulli_lpmf.hpp
index 29fc79d1bf8..9b54f02d839 100644
--- a/stan/math/prim/prob/bernoulli_lpmf.hpp
+++ b/stan/math/prim/prob/bernoulli_lpmf.hpp
@@ -59,7 +59,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 c82ba7d40b5..523fdc1c138 100644
--- a/stan/math/prim/prob/beta_proportion_lpdf.hpp
+++ b/stan/math/prim/prob/beta_proportion_lpdf.hpp
@@ -89,7 +89,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 d530ee15a4b..08e42e2873a 100644
--- a/stan/math/prim/prob/binomial_logit_lpmf.hpp
+++ b/stan/math/prim/prob/binomial_logit_lpmf.hpp
@@ -87,10 +87,11 @@ 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 59dba34dcfc..6e16581e786 100644
--- a/stan/math/prim/prob/cauchy_lpdf.hpp
+++ b/stan/math/prim/prob/cauchy_lpdf.hpp
@@ -84,7 +84,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 d2044d75986..885d9a3c988 100644
--- a/stan/math/prim/prob/chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/chi_square_lpdf.hpp
@@ -74,12 +74,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 e15dfeceebe..3db919362d2 100644
--- a/stan/math/prim/prob/double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/double_exponential_lpdf.hpp
@@ -83,7 +83,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 b2db8f387c9..ba3eddbcbde 100644
--- a/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_lpdf.hpp
@@ -79,7 +79,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 d3e9d09deb4..c48c90ab9d3 100644
--- a/stan/math/prim/prob/exponential_lccdf.hpp
+++ b/stan/math/prim/prob/exponential_lccdf.hpp
@@ -47,7 +47,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 {
@@ -59,7 +59,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 4328c1af2d6..18f3bbfd0da 100644
--- a/stan/math/prim/prob/exponential_lpdf.hpp
+++ b/stan/math/prim/prob/exponential_lpdf.hpp
@@ -76,7 +76,7 @@ return_type_t<T_y, T_inv_scale> exponential_lpdf(const T_y& y,
 
   T_partials_return logp(0.0);
   if (include_summand<propto, T_inv_scale>::value) {
-    logp = sum(log(beta_val)) * max_size(y, beta) / size(beta);
+    logp = sum(log(beta_val)) * max_size(y, beta) / math::size(beta);
   }
   if (include_summand<propto, T_y, T_inv_scale>::value) {
     logp -= sum(beta_val * y_val);
@@ -87,7 +87,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 fe0f6783bef..b305f0a8f20 100644
--- a/stan/math/prim/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/prob/frechet_lpdf.hpp
@@ -72,7 +72,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 ba289902c9d..50e7a6b65c8 100644
--- a/stan/math/prim/prob/gamma_lpdf.hpp
+++ b/stan/math/prim/prob/gamma_lpdf.hpp
@@ -92,7 +92,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 998e4140f37..22814b8982b 100644
--- a/stan/math/prim/prob/gumbel_lpdf.hpp
+++ b/stan/math/prim/prob/gumbel_lpdf.hpp
@@ -78,7 +78,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_lcdf.hpp b/stan/math/prim/prob/inv_chi_square_lcdf.hpp
index abd879cb2c6..b48922fbd63 100644
--- a/stan/math/prim/prob/inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lcdf.hpp
@@ -71,9 +71,9 @@ return_type_t<T_y, T_dof> inv_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/inv_chi_square_lpdf.hpp b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
index aff4dece7f8..d47c4355017 100644
--- a/stan/math/prim/prob/inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
@@ -82,11 +82,12 @@ 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 ceee125329f..0184a6ba572 100644
--- a/stan/math/prim/prob/inv_gamma_lpdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lpdf.hpp
@@ -83,7 +83,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 5ed1b4982ba..f6359a50142 100644
--- a/stan/math/prim/prob/logistic_lpdf.hpp
+++ b/stan/math/prim/prob/logistic_lpdf.hpp
@@ -66,7 +66,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 4afdbae67b5..422b19d65d6 100644
--- a/stan/math/prim/prob/lognormal_lpdf.hpp
+++ b/stan/math/prim/prob/lognormal_lpdf.hpp
@@ -74,10 +74,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 e5a1be8265a..5d688964549 100644
--- a/stan/math/prim/prob/normal_lpdf.hpp
+++ b/stan/math/prim/prob/normal_lpdf.hpp
@@ -85,7 +85,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 475a13aefd1..abf5924b75d 100644
--- a/stan/math/prim/prob/normal_sufficient_lpdf.hpp
+++ b/stan/math/prim/prob/normal_sufficient_lpdf.hpp
@@ -107,7 +107,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);
@@ -141,7 +141,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 ccc92be4bd5..858c6c720ce 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 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 65f7d0e646d..e003450d7d7 100644
--- a/stan/math/prim/prob/pareto_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_lpdf.hpp
@@ -63,7 +63,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 3946f0a29d8..37500c11fd8 100644
--- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp
+++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp
@@ -66,10 +66,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 44c884d8839..85be278124c 100644
--- a/stan/math/prim/prob/poisson_log_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_log_lpmf.hpp
@@ -73,10 +73,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 85732a158ed..01ab3832a95 100644
--- a/stan/math/prim/prob/poisson_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_lpmf.hpp
@@ -68,10 +68,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 aebd139c317..1e447b33956 100644
--- a/stan/math/prim/prob/rayleigh_lpdf.hpp
+++ b/stan/math/prim/prob/rayleigh_lpdf.hpp
@@ -56,10 +56,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_double_exponential_lpdf.hpp b/stan/math/prim/prob/skew_double_exponential_lpdf.hpp
index c14ba95f0c1..d3b6d5cea3e 100644
--- a/stan/math/prim/prob/skew_double_exponential_lpdf.hpp
+++ b/stan/math/prim/prob/skew_double_exponential_lpdf.hpp
@@ -97,10 +97,10 @@ return_type_t<T_y, T_loc, T_scale, T_skewness> skew_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);
   }
   if (include_summand<propto, T_skewness>::value) {
-    logp += sum(log(tau_val) + log1m(tau_val)) * N / size(tau);
+    logp += sum(log(tau_val) + log1m(tau_val)) * N / math::size(tau);
   }
 
   if (!is_constant_all<T_y, T_loc>::value) {
diff --git a/stan/math/prim/prob/skew_normal_lpdf.hpp b/stan/math/prim/prob/skew_normal_lpdf.hpp
index c8481abf6c2..d811f507ae5 100644
--- a/stan/math/prim/prob/skew_normal_lpdf.hpp
+++ b/stan/math/prim/prob/skew_normal_lpdf.hpp
@@ -76,7 +76,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 7df5ea97a22..b6ee5f3f176 100644
--- a/stan/math/prim/prob/student_t_lpdf.hpp
+++ b/stan/math/prim/prob/student_t_lpdf.hpp
@@ -111,10 +111,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 174a847b0db..83e8b22e0d8 100644
--- a/stan/math/prim/prob/uniform_lcdf.hpp
+++ b/stan/math/prim/prob/uniform_lcdf.hpp
@@ -68,7 +68,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);
     }
@@ -80,7 +80,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 cbb2d50261a..e19290bad1c 100644
--- a/stan/math/prim/prob/uniform_lpdf.hpp
+++ b/stan/math/prim/prob/uniform_lpdf.hpp
@@ -95,7 +95,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;
       }
@@ -103,7 +103,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 dec4681dc87..dcaef3036b1 100644
--- a/stan/math/prim/prob/von_mises_lpdf.hpp
+++ b/stan/math/prim/prob/von_mises_lpdf.hpp
@@ -63,7 +63,8 @@ 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 5f69eff760d..b0eeeb43ea7 100644
--- a/stan/math/prim/prob/weibull_lpdf.hpp
+++ b/stan/math/prim/prob/weibull_lpdf.hpp
@@ -89,7 +89,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/algebra_solver_newton.hpp b/stan/math/rev/functor/algebra_solver_newton.hpp
index d77921cd90c..04b426b212d 100644
--- a/stan/math/rev/functor/algebra_solver_newton.hpp
+++ b/stan/math/rev/functor/algebra_solver_newton.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/eval.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/algebra_solver_adapter.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <unsupported/Eigen/NonLinearOptimization>
 #include <iostream>
 #include <string>
@@ -142,7 +143,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
     const int64_t max_num_steps, const T_Args&... args) {
   const auto& x_ref = to_ref(value_of(x));
   auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...));
-  auto args_vals_tuple = apply(
+  auto args_vals_tuple = math::apply(
       [&](const auto&... args) {
         return std::make_tuple(to_ref(value_of(args))...);
       },
@@ -157,7 +158,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
   check_positive("algebra_solver_newton", "max_num_steps", max_num_steps);
 
   // Solve the system
-  Eigen::VectorXd theta_dbl = apply(
+  Eigen::VectorXd theta_dbl = math::apply(
       [&](const auto&... vals) {
         return kinsol_solve(f, x_ref, scaling_step_size, function_tolerance,
                             max_num_steps, 1, 10, KIN_LINESEARCH, msgs,
@@ -166,8 +167,8 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
       args_vals_tuple);
 
   auto f_wrt_x = [&](const auto& x) {
-    return apply([&](const auto&... args) { return f(x, msgs, args...); },
-                 args_vals_tuple);
+    return math::apply([&](const auto&... args) { return f(x, msgs, args...); },
+                       args_vals_tuple);
   };
 
   Eigen::MatrixXd Jf_x;
@@ -190,7 +191,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_newton_impl(
           nested_rev_autodiff rev;
 
           Eigen::VectorXd ret_val = ret.val();
-          auto x_nrad_ = apply(
+          auto x_nrad_ = math::apply(
               [&ret_val, &f, msgs](const auto&... args) {
                 return eval(f(ret_val, msgs, args...));
               },
diff --git a/stan/math/rev/functor/algebra_solver_powell.hpp b/stan/math/rev/functor/algebra_solver_powell.hpp
index 5415f361cad..2d09f97ea8c 100644
--- a/stan/math/rev/functor/algebra_solver_powell.hpp
+++ b/stan/math/rev/functor/algebra_solver_powell.hpp
@@ -136,7 +136,7 @@ Eigen::VectorXd algebra_solver_powell_impl(const F& f, const T& x,
   auto args_vals_tuple = std::make_tuple(to_ref(args)...);
 
   auto f_wrt_x = [&args_vals_tuple, &f, msgs](const auto& x) {
-    return apply(
+    return math::apply(
         [&x, &f, msgs](const auto&... args) { return f(x, msgs, args...); },
         args_vals_tuple);
   };
@@ -338,14 +338,14 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
     const int64_t max_num_steps, const T_Args&... args) {
   auto x_ref = eval(value_of(x));
   auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...));
-  auto args_vals_tuple = apply(
+  auto args_vals_tuple = math::apply(
       [&](const auto&... args) {
         return std::make_tuple(to_ref(value_of(args))...);
       },
       *arena_args_tuple);
 
   auto f_wrt_x = [&args_vals_tuple, &f, msgs](const auto& x) {
-    return apply(
+    return math::apply(
         [&x, &f, msgs](const auto&... args) { return f(x, msgs, args...); },
         args_vals_tuple);
   };
@@ -384,7 +384,7 @@ Eigen::Matrix<var, Eigen::Dynamic, 1> algebra_solver_powell_impl(
     {
       nested_rev_autodiff rev;
       Eigen::VectorXd ret_val = ret.val();
-      auto x_nrad_ = apply(
+      auto x_nrad_ = math::apply(
           [&](const auto&... args) { return eval(f(ret_val, msgs, args...)); },
           *arena_args_tuple);
       x_nrad_.adj() = eta;
diff --git a/stan/math/rev/functor/coupled_ode_system.hpp b/stan/math/rev/functor/coupled_ode_system.hpp
index 1f043fe0e9a..220dda25a4c 100644
--- a/stan/math/rev/functor/coupled_ode_system.hpp
+++ b/stan/math/rev/functor/coupled_ode_system.hpp
@@ -6,6 +6,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/value_of.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stdexcept>
@@ -125,9 +126,9 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
     for (size_t n = 0; n < N_; ++n)
       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...); },
-                local_args_tuple_);
+    Eigen::Matrix<var, Eigen::Dynamic, 1> f_y_t_vars = 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",
                      N_);
@@ -143,7 +144,7 @@ struct coupled_ode_system_impl<false, F, T_y0, Args...> {
                sizeof(double) * args_adjoints_.size());
       }
 
-      apply(
+      math::apply(
           [&](auto&&... args) {
             accumulate_adjoints(args_adjoints_.data(), args...);
           },
diff --git a/stan/math/rev/functor/cvodes_integrator.hpp b/stan/math/rev/functor/cvodes_integrator.hpp
index 637526938e7..b36bae594e0 100644
--- a/stan/math/rev/functor/cvodes_integrator.hpp
+++ b/stan/math/rev/functor/cvodes_integrator.hpp
@@ -6,6 +6,7 @@
 #include <stan/math/rev/functor/coupled_ode_system.hpp>
 #include <stan/math/rev/functor/ode_store_sensitivities.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <sundials/sundials_context.h>
 #include <cvodes/cvodes.h>
@@ -104,9 +105,9 @@ class cvodes_integrator {
   inline void rhs(double t, const double y[], double dy_dt[]) const {
     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...); },
-                value_of_args_tuple_);
+    Eigen::VectorXd dy_dt_vec = 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",
                      N_);
@@ -123,8 +124,9 @@ 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...); },
-                   value_of_args_tuple_);
+      return math::apply(
+          [&](auto&&... args) { return f_(t, y, msgs_, args...); },
+          value_of_args_tuple_);
     };
 
     jacobian(f_wrapped, Eigen::Map<const Eigen::VectorXd>(y, N_), fy, Jfy);
@@ -214,7 +216,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),
@@ -313,7 +315,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/cvodes_integrator_adjoint.hpp b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
index 117b869e89f..0c23f0758bf 100644
--- a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
+++ b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
@@ -7,6 +7,7 @@
 #include <stan/math/rev/functor/ode_store_sensitivities.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
 #include <sundials/sundials_context.h>
 #include <cvodes/cvodes.h>
@@ -544,7 +545,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
   template <typename yT, typename... ArgsT>
   constexpr auto rhs(double t, const yT& y,
                      const std::tuple<ArgsT...>& args_tuple) const {
-    return apply(
+    return math::apply(
         [&](auto&&... args) { return solver_->f_(t, y, msgs_, args...); },
         args_tuple);
   }
@@ -641,7 +642,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
                      f_y_t_vars.size(), "states", N_);
     f_y_t_vars.adj() = -Eigen::Map<Eigen::VectorXd>(NV_DATA_S(yB), N_);
     grad();
-    apply(
+    math::apply(
         [&qBdot](auto&&... args) {
           accumulate_adjoints(NV_DATA_S(qBdot), args...);
         },
diff --git a/stan/math/rev/functor/dae.hpp b/stan/math/rev/functor/dae.hpp
index b3ca6d2c6dc..77538e1f64e 100644
--- a/stan/math/rev/functor/dae.hpp
+++ b/stan/math/rev/functor/dae.hpp
@@ -1,6 +1,7 @@
 #ifndef STAN_MATH_REV_FUNCTOR_DAE_HPP
 #define STAN_MATH_REV_FUNCTOR_DAE_HPP
 
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/functor/idas_integrator.hpp>
 #include <stan/math/rev/functor/dae_system.hpp>
@@ -66,14 +67,14 @@ dae_tol_impl(const char* func, const F& f, const T_yy& yy0, const T_yp& yp0,
   check_positive(func, "max_num_steps", max_num_steps);
 
   const auto& args_ref_tuple = std::make_tuple(to_ref(args)...);
-  apply(
+  math::apply(
       [&](auto&&... args) {
         std::vector<int> unused_temp{
             0, (check_finite("dae", "DAE parameters and data", args), 0)...};
       },
       args_ref_tuple);
 
-  return apply(
+  return math::apply(
       [&](const auto&... args_refs) {
         dae_system<F, T_yy, T_yp, ref_type_t<T_Args>...> dae(f, yy0, yp0, msgs,
                                                              args_refs...);
diff --git a/stan/math/rev/functor/dae_system.hpp b/stan/math/rev/functor/dae_system.hpp
index 4e694524b7f..2a1f97f2bd8 100644
--- a/stan/math/rev/functor/dae_system.hpp
+++ b/stan/math/rev/functor/dae_system.hpp
@@ -99,7 +99,7 @@ class dae_system {
   }
 
   void eval_residual(double t) {
-    dbl_rr = apply(
+    dbl_rr = math::apply(
         [&](auto&&... args) { return f_(t, dbl_yy, dbl_yp, msgs_, args...); },
         dbl_args_tuple_);
   }
@@ -144,7 +144,7 @@ class dae_system {
     }
 
     Eigen::VectorXd g(m);
-    Eigen::Matrix<var, -1, 1> fy = apply(
+    Eigen::Matrix<var, -1, 1> fy = math::apply(
         [&](auto&&... args) {
           return dae->f_(t, yy_var, yp_var, dae->msgs_, args...);
         },
diff --git a/stan/math/rev/functor/integrate_1d.hpp b/stan/math/rev/functor/integrate_1d.hpp
index be8ff106393..60607db7894 100644
--- a/stan/math/rev/functor/integrate_1d.hpp
+++ b/stan/math/rev/functor/integrate_1d.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/integrate_1d.hpp>
 #include <cmath>
 #include <functional>
@@ -57,7 +58,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
 
     double integral = integrate(
         [&](const auto &x, const auto &xc) {
-          return apply(
+          return math::apply(
               [&](auto &&... val_args) { return f(x, xc, msgs, val_args...); },
               args_val_tuple);
         },
@@ -79,7 +80,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
     }
 
     if (is_var<T_a>::value && !is_inf(a)) {
-      *partials_ptr = apply(
+      *partials_ptr = math::apply(
           [&f, a_val, msgs](auto &&... val_args) {
             return -f(a_val, 0.0, msgs, val_args...);
           },
@@ -88,7 +89,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
     }
 
     if (!is_inf(b) && is_var<T_b>::value) {
-      *partials_ptr = apply(
+      *partials_ptr = math::apply(
           [&f, b_val, msgs](auto &&... val_args) {
             return f(b_val, 0.0, msgs, val_args...);
           },
@@ -104,7 +105,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
 
       // Save the varis so it's easy to efficiently access the nth adjoint
       std::vector<vari *> local_varis(num_vars_args);
-      apply(
+      math::apply(
           [&](const auto &... args) {
             save_varis(local_varis.data(), args...);
           },
@@ -118,7 +119,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_impl(
               argument_nest.set_zero_all_adjoints();
 
               nested_rev_autodiff gradient_nest;
-              var fx = apply(
+              var fx = math::apply(
                   [&f, &x, &xc, msgs](auto &&... local_args) {
                     return f(x, xc, msgs, local_args...);
                   },
diff --git a/stan/math/rev/functor/kinsol_data.hpp b/stan/math/rev/functor/kinsol_data.hpp
index 73f74d9e684..a0847cf9fbd 100644
--- a/stan/math/rev/functor/kinsol_data.hpp
+++ b/stan/math/rev/functor/kinsol_data.hpp
@@ -74,7 +74,7 @@ class kinsol_system_data {
 
     Eigen::Map<Eigen::VectorXd> f_eval_map(N_VGetArrayPointer(f_eval),
                                            explicit_system->N_);
-    auto result = apply(
+    auto result = math::apply(
         [&](const auto&... args) {
           return explicit_system->f_(x_eigen, explicit_system->msgs_, args...);
         },
@@ -109,7 +109,7 @@ class kinsol_system_data {
                                            explicit_system->N_);
 
     auto f_wrt_x = [&](const auto& x) {
-      return apply(
+      return math::apply(
           [&](const auto&... args) {
             return explicit_system->f_(x, explicit_system->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..009bfe75c00 100644
--- a/stan/math/rev/functor/ode_bdf.hpp
+++ b/stan/math/rev/functor/ode_bdf.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>
 
@@ -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 8002a320117..da9ac216880 100644
--- a/stan/math/rev/functor/reduce_sum.hpp
+++ b/stan/math/rev/functor/reduce_sum.hpp
@@ -114,7 +114,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>(
@@ -141,7 +141,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...);
@@ -159,7 +159,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...);
           },
diff --git a/test/unit/math/prim/meta/StdVectorBuilder_test.cpp b/test/unit/math/prim/meta/StdVectorBuilder_test.cpp
index 45707477ec7..cf9d082560d 100644
--- a/test/unit/math/prim/meta/StdVectorBuilder_test.cpp
+++ b/test/unit/math/prim/meta/StdVectorBuilder_test.cpp
@@ -9,7 +9,7 @@ TEST(MathMetaPrim, StdVectorBuilder_false_false_scalar) {
 
   double a_double(1);
 
-  StdVectorBuilder<false, double, double> dvv1(size(a_double));
+  StdVectorBuilder<false, double, double> dvv1(stan::math::size(a_double));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -20,7 +20,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_scalar) {
 
   double a_double(1);
 
-  StdVectorBuilder<true, double, double> dvv1(size(a_double));
+  StdVectorBuilder<true, double, double> dvv1(stan::math::size(a_double));
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   EXPECT_FLOAT_EQ(0.0, dvv1[1]);
   EXPECT_FLOAT_EQ(0.0, dvv1[100]);
@@ -47,7 +47,7 @@ TEST(MathMetaPrim, StdVectorBuilder_false_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  StdVectorBuilder<false, double, double> dvv2(size(a_std_vector));
+  StdVectorBuilder<false, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -59,7 +59,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  StdVectorBuilder<true, double, double> dvv2(size(a_std_vector));
+  StdVectorBuilder<true, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_FLOAT_EQ(0.0, dvv2[0]);
   EXPECT_FLOAT_EQ(0.0, dvv2[1]);
   EXPECT_FLOAT_EQ(0.0, dvv2[2]);
@@ -89,11 +89,11 @@ TEST(MathMetaPrim, StdVectorBuilder_false_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  StdVectorBuilder<false, double, double> dvv3(size(a_vector));
+  StdVectorBuilder<false, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  StdVectorBuilder<false, double, double> dvv4(size(a_row_vector));
+  StdVectorBuilder<false, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv4.data(), std::logic_error);
 }
@@ -107,7 +107,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  StdVectorBuilder<true, double, double> dvv3(size(a_vector));
+  StdVectorBuilder<true, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_FLOAT_EQ(0.0, dvv3[0]);
   EXPECT_FLOAT_EQ(0.0, dvv3[1]);
   EXPECT_FLOAT_EQ(0.0, dvv3[2]);
@@ -115,7 +115,7 @@ TEST(MathMetaPrim, StdVectorBuilder_true_false_matrix) {
   EXPECT_NO_THROW(data3 = dvv3.data());
   EXPECT_FLOAT_EQ(0.0, data3);
 
-  StdVectorBuilder<true, double, double> dvv4(size(a_row_vector));
+  StdVectorBuilder<true, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_FLOAT_EQ(0.0, dvv4[0]);
   EXPECT_FLOAT_EQ(0.0, dvv4[1]);
   EXPECT_FLOAT_EQ(0.0, dvv4[2]);
diff --git a/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp b/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp
index 8c0a105f4fe..a8be33072b4 100644
--- a/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp
+++ b/test/unit/math/prim/meta/VectorBuilderHelper_test.cpp
@@ -9,7 +9,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_false_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilderHelper<double, false, false> dvv1(size(a_double));
+  VectorBuilderHelper<double, false, false> dvv1(stan::math::size(a_double));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -20,7 +20,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilderHelper<double, true, false> dvv1(size(a_double));
+  VectorBuilderHelper<double, true, false> dvv1(stan::math::size(a_double));
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   EXPECT_FLOAT_EQ(0.0, dvv1[1]);
   EXPECT_FLOAT_EQ(0.0, dvv1[100]);
@@ -36,7 +36,8 @@ TEST(MathMetaPrim, VectorBuilderHelper_false_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilderHelper<double, false, false> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, false, false> dvv2(
+      stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -48,7 +49,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilderHelper<double, true, false> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, true, false> dvv2(stan::math::size(a_std_vector));
   EXPECT_FLOAT_EQ(0.0, dvv2[0]);
   EXPECT_FLOAT_EQ(0.0, dvv2[1]);
   EXPECT_FLOAT_EQ(0.0, dvv2[2]);
@@ -66,11 +67,12 @@ TEST(MathMetaPrim, VectorBuilderHelper_false_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, false, false> dvv3(size(a_vector));
+  VectorBuilderHelper<double, false, false> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  VectorBuilderHelper<double, false, false> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, false, false> dvv4(
+      stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv4.data(), std::logic_error);
 }
@@ -84,7 +86,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, true, false> dvv3(size(a_vector));
+  VectorBuilderHelper<double, true, false> dvv3(stan::math::size(a_vector));
   EXPECT_FLOAT_EQ(0.0, dvv3[0]);
   EXPECT_FLOAT_EQ(0.0, dvv3[1]);
   EXPECT_FLOAT_EQ(0.0, dvv3[2]);
@@ -92,7 +94,7 @@ TEST(MathMetaPrim, VectorBuilderHelper_true_false_matrix) {
   EXPECT_NO_THROW(data3 = dvv3.data());
   EXPECT_FLOAT_EQ(0.0, data3);
 
-  VectorBuilderHelper<double, true, false> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, true, false> dvv4(stan::math::size(a_row_vector));
   EXPECT_FLOAT_EQ(0.0, dvv4[0]);
   EXPECT_FLOAT_EQ(0.0, dvv4[1]);
   EXPECT_FLOAT_EQ(0.0, dvv4[2]);
diff --git a/test/unit/math/prim/meta/VectorBuilder_test.cpp b/test/unit/math/prim/meta/VectorBuilder_test.cpp
index a1dca9ae971..e48eb3dbc74 100644
--- a/test/unit/math/prim/meta/VectorBuilder_test.cpp
+++ b/test/unit/math/prim/meta/VectorBuilder_test.cpp
@@ -9,7 +9,7 @@ TEST(MathMetaPrim, VectorBuilder_false_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilder<false, double, double> dvv1(size(a_double));
+  VectorBuilder<false, double, double> dvv1(stan::math::size(a_double));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -20,7 +20,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_scalar) {
 
   double a_double(1);
 
-  VectorBuilder<true, double, double> dvv1(size(a_double));
+  VectorBuilder<true, double, double> dvv1(stan::math::size(a_double));
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   EXPECT_FLOAT_EQ(0.0, dvv1[1]);
   EXPECT_FLOAT_EQ(0.0, dvv1[100]);
@@ -36,7 +36,7 @@ TEST(MathMetaPrim, VectorBuilder_false_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilder<false, double, double> dvv2(size(a_std_vector));
+  VectorBuilder<false, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -48,7 +48,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_vector) {
 
   std::vector<double> a_std_vector(3);
 
-  VectorBuilder<true, double, double> dvv2(size(a_std_vector));
+  VectorBuilder<true, double, double> dvv2(stan::math::size(a_std_vector));
   EXPECT_FLOAT_EQ(0.0, dvv2[0]);
   EXPECT_FLOAT_EQ(0.0, dvv2[1]);
   EXPECT_FLOAT_EQ(0.0, dvv2[2]);
@@ -66,11 +66,11 @@ TEST(MathMetaPrim, VectorBuilder_false_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<false, double, double> dvv3(size(a_vector));
+  VectorBuilder<false, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  VectorBuilder<false, double, double> dvv4(size(a_row_vector));
+  VectorBuilder<false, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv4.data(), std::logic_error);
 }
@@ -84,7 +84,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_matrix) {
   Matrix<double, Dynamic, 1> a_vector(4);
   Matrix<double, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<true, double, double> dvv3(size(a_vector));
+  VectorBuilder<true, double, double> dvv3(stan::math::size(a_vector));
   EXPECT_FLOAT_EQ(0.0, dvv3[0]);
   EXPECT_FLOAT_EQ(0.0, dvv3[1]);
   EXPECT_FLOAT_EQ(0.0, dvv3[2]);
@@ -92,7 +92,7 @@ TEST(MathMetaPrim, VectorBuilder_true_false_matrix) {
   EXPECT_NO_THROW(data3 = dvv3.data());
   EXPECT_FLOAT_EQ(0.0, data3);
 
-  VectorBuilder<true, double, double> dvv4(size(a_row_vector));
+  VectorBuilder<true, double, double> dvv4(stan::math::size(a_row_vector));
   EXPECT_FLOAT_EQ(0.0, dvv4[0]);
   EXPECT_FLOAT_EQ(0.0, dvv4[1]);
   EXPECT_FLOAT_EQ(0.0, dvv4[2]);
diff --git a/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp b/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp
index ac845bbfe35..af0d0234e78 100644
--- a/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp
+++ b/test/unit/math/rev/meta/VectorBuilderHelper_test.cpp
@@ -11,7 +11,7 @@ TEST(MetaTraitsRevScal, VectorBuilderHelper_false_true) {
 
   var a_var(1);
 
-  VectorBuilderHelper<double, false, true> dvv1(size(a_var));
+  VectorBuilderHelper<double, false, true> dvv1(stan::math::size(a_var));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 }
@@ -24,7 +24,7 @@ TEST(MetaTraitsRevArr, VectorBuilderHelper_false_true) {
 
   std::vector<var> a_std_vector(3);
 
-  VectorBuilderHelper<double, false, true> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, false, true> dvv2(stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -37,14 +37,14 @@ TEST(MetaTraitsRevArr, VectorBuilderHelper_true_true) {
 
   var a_var(1);
   std::vector<var> a_std_vector(3);
-  VectorBuilderHelper<double, true, true> dvv1(size(a_var));
+  VectorBuilderHelper<double, true, true> dvv1(stan::math::size(a_var));
   dvv1[0] = 0.0;
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   std::vector<double> data1;
   EXPECT_NO_THROW(data1 = dvv1.data());
-  EXPECT_EQ(size(a_var), data1.size());
+  EXPECT_EQ(stan::math::size(a_var), data1.size());
 
-  VectorBuilderHelper<double, true, true> dvv2(size(a_std_vector));
+  VectorBuilderHelper<double, true, true> dvv2(stan::math::size(a_std_vector));
   dvv2[0] = 0.0;
   dvv2[1] = 1.0;
   dvv2[2] = 2.0;
@@ -53,7 +53,7 @@ TEST(MetaTraitsRevArr, VectorBuilderHelper_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv2[2]);
   std::vector<double> data2;
   EXPECT_NO_THROW(data2 = dvv2.data());
-  EXPECT_EQ(size(a_std_vector), data2.size());
+  EXPECT_EQ(stan::math::size(a_std_vector), data2.size());
 }
 
 TEST(MetaTraitsRevMat, VectorBuilderHelper_false_true) {
@@ -66,11 +66,11 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_false_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, false, true> dvv3(size(a_vector));
+  VectorBuilderHelper<double, false, true> dvv3(stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
-  VectorBuilderHelper<double, false, true> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, false, true> dvv4(stan::math::size(a_row_vector));
   EXPECT_THROW(dvv4[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 }
@@ -85,7 +85,7 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_true_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilderHelper<double, true, true> dvv3(size(a_vector));
+  VectorBuilderHelper<double, true, true> dvv3(stan::math::size(a_vector));
   dvv3[0] = 0.0;
   dvv3[1] = 1.0;
   dvv3[2] = 2.0;
@@ -94,9 +94,9 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv3[2]);
   std::vector<double> data3;
   EXPECT_NO_THROW(data3 = dvv3.data());
-  EXPECT_EQ(size(a_vector), data3.size());
+  EXPECT_EQ(stan::math::size(a_vector), data3.size());
 
-  VectorBuilderHelper<double, true, true> dvv4(size(a_row_vector));
+  VectorBuilderHelper<double, true, true> dvv4(stan::math::size(a_row_vector));
   dvv4[0] = 0.0;
   dvv4[1] = 1.0;
   dvv4[2] = 2.0;
@@ -105,5 +105,5 @@ TEST(MetaTraitsRevMat, VectorBuilderHelper_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv4[2]);
   std::vector<double> data4;
   EXPECT_NO_THROW(data4 = dvv4.data());
-  EXPECT_EQ(size(a_row_vector), data4.size());
+  EXPECT_EQ(stan::math::size(a_row_vector), data4.size());
 }
diff --git a/test/unit/math/rev/meta/VectorBuilder_test.cpp b/test/unit/math/rev/meta/VectorBuilder_test.cpp
index 481870f24a8..edd3df52fdc 100644
--- a/test/unit/math/rev/meta/VectorBuilder_test.cpp
+++ b/test/unit/math/rev/meta/VectorBuilder_test.cpp
@@ -13,11 +13,12 @@ TEST(MetaTraitsRevArr, VectorBuilder_false_true) {
   var a_var(1);
   std::vector<var> a_std_vector(3);
 
-  VectorBuilder<false, double, std::vector<var> > dvv1(size(a_var));
+  VectorBuilder<false, double, std::vector<var> > dvv1(stan::math::size(a_var));
   EXPECT_THROW(dvv1[0], std::logic_error);
   EXPECT_THROW(dvv1.data(), std::logic_error);
 
-  VectorBuilder<false, double, std::vector<var> > dvv2(size(a_std_vector));
+  VectorBuilder<false, double, std::vector<var> > dvv2(
+      stan::math::size(a_std_vector));
   EXPECT_THROW(dvv2[0], std::logic_error);
   EXPECT_THROW(dvv2.data(), std::logic_error);
 }
@@ -31,14 +32,15 @@ TEST(MetaTraitsRevArr, VectorBuilder_true_true) {
   var a_var(1);
   std::vector<var> a_std_vector(3);
 
-  VectorBuilder<true, double, std::vector<var> > dvv1(size(a_var));
+  VectorBuilder<true, double, std::vector<var> > dvv1(stan::math::size(a_var));
   dvv1[0] = 0.0;
   EXPECT_FLOAT_EQ(0.0, dvv1[0]);
   std::vector<double> data1;
   EXPECT_NO_THROW(data1 = dvv1.data());
-  EXPECT_EQ(size(a_var), data1.size());
+  EXPECT_EQ(stan::math::size(a_var), data1.size());
 
-  VectorBuilder<true, double, std::vector<var> > dvv2(size(a_std_vector));
+  VectorBuilder<true, double, std::vector<var> > dvv2(
+      stan::math::size(a_std_vector));
   dvv2[0] = 0.0;
   dvv2[1] = 1.0;
   dvv2[2] = 2.0;
@@ -47,7 +49,7 @@ TEST(MetaTraitsRevArr, VectorBuilder_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv2[2]);
   std::vector<double> data2;
   EXPECT_NO_THROW(data2 = dvv2.data());
-  EXPECT_EQ(size(a_std_vector), data2.size());
+  EXPECT_EQ(stan::math::size(a_std_vector), data2.size());
 }
 
 TEST(MetaTraitsRevMat, VectorBuilder_false_true) {
@@ -60,7 +62,8 @@ TEST(MetaTraitsRevMat, VectorBuilder_false_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<false, double, Matrix<var, Dynamic, 1> > dvv3(size(a_vector));
+  VectorBuilder<false, double, Matrix<var, Dynamic, 1> > dvv3(
+      stan::math::size(a_vector));
   EXPECT_THROW(dvv3[0], std::logic_error);
   EXPECT_THROW(dvv3.data(), std::logic_error);
 
@@ -80,7 +83,8 @@ TEST(MetaTraitsRevMat, VectorBuilder_true_true) {
   Matrix<var, Dynamic, 1> a_vector(4);
   Matrix<var, 1, Dynamic> a_row_vector(5);
 
-  VectorBuilder<true, double, Matrix<var, Dynamic, 1> > dvv3(size(a_vector));
+  VectorBuilder<true, double, Matrix<var, Dynamic, 1> > dvv3(
+      stan::math::size(a_vector));
   dvv3[0] = 0.0;
   dvv3[1] = 1.0;
   dvv3[2] = 2.0;
@@ -89,7 +93,7 @@ TEST(MetaTraitsRevMat, VectorBuilder_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv3[2]);
   std::vector<double> data3;
   EXPECT_NO_THROW(data3 = dvv3.data());
-  EXPECT_EQ(size(a_vector), data3.size());
+  EXPECT_EQ(stan::math::size(a_vector), data3.size());
 
   VectorBuilder<true, double, Matrix<var, 1, Dynamic> > dvv4(
       stan::math::size(a_row_vector));
@@ -101,5 +105,5 @@ TEST(MetaTraitsRevMat, VectorBuilder_true_true) {
   EXPECT_FLOAT_EQ(2.0, dvv4[2]);
   std::vector<double> data4;
   EXPECT_NO_THROW(data4 = dvv4.data());
-  EXPECT_EQ(size(a_row_vector), data4.size());
+  EXPECT_EQ(stan::math::size(a_row_vector), data4.size());
 }