diff --git a/make/compiler_flags b/make/compiler_flags
index d6f41aec41b..26e15a0e989 100644
--- a/make/compiler_flags
+++ b/make/compiler_flags
@@ -119,33 +119,9 @@ INC_GTEST ?= -I $(GTEST)/include -I $(GTEST)
 ## setup precompiler options
 CPPFLAGS_BOOST ?= -DBOOST_DISABLE_ASSERTS
 CPPFLAGS_SUNDIALS ?= -DNO_FPRINTF_OUTPUT $(CPPFLAGS_OPTIM_SUNDIALS) $(CXXFLAGS_FLTO_SUNDIALS)
-#CPPFLAGS_GTEST ?=
-STAN_HAS_CXX17 ?= false
-ifeq ($(CXX_TYPE), gcc)
-  GCC_GE_73 := $(shell [ $(CXX_MAJOR) -gt 7 -o \( $(CXX_MAJOR) -eq 7 -a $(CXX_MINOR) -ge 1 \) ] && echo true)
-  ifeq ($(GCC_GE_73),true)
-    STAN_HAS_CXX17 := true
-  endif
-else ifeq ($(CXX_TYPE), clang)
-  CLANG_GE_5 := $(shell [ $(CXX_MAJOR) -gt 5 -o \( $(CXX_MAJOR) -eq 5 -a $(CXX_MINOR) -ge 0 \) ] && echo true)
-  ifeq ($(CLANG_GE_5),true)
-    STAN_HAS_CXX17 := true
-  endif
-else ifeq ($(CXX_TYPE), mingw32-gcc)
-  MINGW_GE_50 := $(shell [ $(CXX_MAJOR) -gt 5 -o \( $(CXX_MAJOR) -eq 5 -a $(CXX_MINOR) -ge 0 \) ] && echo true)
-  ifeq ($(MINGW_GE_50),true)
-    STAN_HAS_CXX17 := true
-  endif
-endif
 
-ifeq ($(STAN_HAS_CXX17), true)
-  CXXFLAGS_LANG ?= -std=c++17
-  CXXFLAGS_STANDARD ?= c++17
-else
-  $(warning "Stan cannot detect if your compiler has the C++17 standard. If it does, please set STAN_HAS_CXX17=true in your make/local file. C++17 support is mandatory in the next release of Stan. Defaulting to C++14")
-  CXXFLAGS_LANG ?= -std=c++1y
-  CXXFLAGS_STANDARD ?= c++1y
-endif
+CXXFLAGS_LANG ?= -std=c++17
+CXXFLAGS_STANDARD ?= c++17
 #CXXFLAGS_BOOST ?=
 CXXFLAGS_SUNDIALS ?= -pipe $(CXXFLAGS_OPTIM_SUNDIALS) $(CPPFLAGS_FLTO_SUNDIALS)
 #CXXFLAGS_GTEST
diff --git a/stan/math/fwd/fun/accumulator.hpp b/stan/math/fwd/fun/accumulator.hpp
index e1bbaea5d3c..a5ceb9f6dbc 100644
--- a/stan/math/fwd/fun/accumulator.hpp
+++ b/stan/math/fwd/fun/accumulator.hpp
@@ -6,7 +6,97 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun/sum.hpp>
 #include <stan/math/prim/fun/accumulator.hpp>
+
 #include <vector>
 #include <type_traits>
 
+namespace stan {
+namespace math {
+template <typename T, typename>
+class accumulator;
+/**
+ * Class to accumulate values and eventually return their sum.  If
+ * no values are ever added, the return value is 0.
+ *
+ * This class is useful for speeding up autodiff of long sums
+ * because it uses the <code>sum()</code> operation (either from
+ * <code>stan::math</code> or one defined by argument-dependent lookup.
+ *
+ * @tparam T Type of scalar added
+ */
+template <typename T>
+class accumulator<T, require_fvar_t<T>> {
+ private:
+  std::vector<T> buf_;
+
+ public:
+  /**
+   * Add the specified arithmetic type value to the buffer after
+   * static casting it to the class type <code>T</code>.
+   *
+   * <p>See the std library doc for <code>std::is_arithmetic</code>
+   * for information on what counts as an arithmetic type.
+   *
+   * @tparam S Type of argument
+   * @param x Value to add
+   */
+  template <typename S, typename = require_stan_scalar_t<S>>
+  inline void add(S x) {
+    buf_.push_back(x);
+  }
+
+  /**
+   * Add each entry in the specified matrix, vector, or row vector
+   * of values to the buffer.
+   *
+   * @tparam S type of the matrix
+   * @param m Matrix of values to add
+   */
+  template <typename S, require_matrix_t<S>* = nullptr>
+  inline void add(const S& m) {
+    buf_.push_back(stan::math::sum(m));
+  }
+
+  /**
+   * Recursively add each entry in the specified standard vector
+   * to the buffer.  This will allow vectors of primitives,
+   * autodiff variables to be added; if the vector entries
+   * are collections, their elements are recursively added.
+   *
+   * @tparam S Type of value to recursively add.
+   * @param xs Vector of entries to add
+   */
+  template <typename S>
+  inline void add(const std::vector<S>& xs) {
+    for (size_t i = 0; i < xs.size(); ++i) {
+      this->add(xs[i]);
+    }
+  }
+
+#ifdef STAN_OPENCL
+
+  /**
+   * Sum each entry and then push to the buffer.
+   * @tparam S A Type inheriting from `matrix_cl_base`
+   * @param xs An OpenCL matrix
+   */
+  template <typename S,
+            require_all_kernel_expressions_and_none_scalar_t<S>* = nullptr>
+  inline void add(const S& xs) {
+    buf_.push_back(stan::math::sum(xs));
+  }
+
+#endif
+
+  /**
+   * Return the sum of the accumulated values.
+   *
+   * @return Sum of accumulated values.
+   */
+  inline T sum() const { return stan::math::sum(buf_); }
+};
+
+}  // namespace math
+}  // namespace stan
+
 #endif
diff --git a/stan/math/fwd/fun/pow.hpp b/stan/math/fwd/fun/pow.hpp
index c181df79b62..c2056c9f5cf 100644
--- a/stan/math/fwd/fun/pow.hpp
+++ b/stan/math/fwd/fun/pow.hpp
@@ -3,10 +3,12 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
-#include <stan/math/fwd/fun/sqrt.hpp>
 #include <stan/math/fwd/fun/inv.hpp>
 #include <stan/math/fwd/fun/inv_sqrt.hpp>
 #include <stan/math/fwd/fun/inv_square.hpp>
+#include <stan/math/fwd/fun/log.hpp>
+#include <stan/math/fwd/fun/sqrt.hpp>
+#include <stan/math/fwd/fun/square.hpp>
 #include <stan/math/prim/fun/pow.hpp>
 #include <cmath>
 #include <complex>
@@ -14,204 +16,74 @@
 
 namespace stan {
 namespace math {
-
-template <typename T>
-inline fvar<T> pow(const fvar<T>& x1, const fvar<T>& x2) {
-  using std::log;
-  using std::pow;
-  T pow_x1_x2(pow(x1.val_, x2.val_));
-  return fvar<T>(pow_x1_x2, (x2.d_ * log(x1.val_) + x2.val_ * x1.d_ / x1.val_)
-                                * pow_x1_x2);
-}
-
-template <typename T, typename U, typename = require_arithmetic_t<U>>
-inline fvar<T> pow(U x1, const fvar<T>& x2) {
+/*
+ *
+ * @tparam T1 Either an `fvar`, `arithmetic`, or `complex` type with an inner
+ * `fvar` or `arithmetic` type.
+ * @tparam T2 Either a `fvar`, `arithmetic`, or `complex` type with an inner
+ * `fvar` or `arithmetic` type.
+ * @param x1 Base variable.
+ * @param x2 Exponent variable.
+ * @return Base raised to the exponent.
+ */
+template <typename T1, typename T2,
+          require_any_fvar_t<base_type_t<T1>, base_type_t<T2>>* = nullptr,
+          require_all_stan_scalar_t<T1, T2>* = nullptr>
+inline auto pow(const T1& x1, const T2& x2) {
   using std::log;
   using std::pow;
-  T u = pow(x1, x2.val_);
-  return fvar<T>(u, x2.d_ * log(x1) * u);
-}
-
-template <typename T, typename U, typename = require_arithmetic_t<U>>
-inline fvar<T> pow(const fvar<T>& x1, U x2) {
-  using std::pow;
-  using std::sqrt;
-  if (x2 == -2) {
-    return inv_square(x1);
-  }
-  if (x2 == -1) {
-    return inv(x1);
-  }
-  if (x2 == -0.5) {
-    return inv_sqrt(x1);
-  }
-  if (x2 == 0.5) {
-    return sqrt(x1);
-  }
-  if (x2 == 1.0) {
-    return x1;
+  if constexpr (is_complex<T1>::value || is_complex<T2>::value) {
+    return internal::complex_pow(x1, x2);
+  } else if constexpr (is_fvar<T1>::value && is_fvar<T2>::value) {
+    auto pow_x1_x2(stan::math::pow(x1.val_, x2.val_));
+    return T1(pow_x1_x2,
+              (x2.d_ * stan::math::log(x1.val_) + x2.val_ * x1.d_ / x1.val_)
+                  * pow_x1_x2);
+  } else if constexpr (is_fvar<T2>::value) {
+    auto u = stan::math::pow(x1, x2.val_);
+    return T2(u, x2.d_ * stan::math::log(x1) * u);
+  } else {
+    using std::sqrt;
+    if (x2 == -2) {
+      return stan::math::inv_square(x1);
+    }
+    if (x2 == -1) {
+      return stan::math::inv(x1);
+    }
+    if (x2 == -0.5) {
+      return stan::math::inv_sqrt(x1);
+    }
+    if (x2 == 0.5) {
+      return stan::math::sqrt(x1);
+    }
+    if (x2 == 1.0) {
+      return x1;
+    }
+    if (x2 == 2.0) {
+      return stan::math::square(x1);
+    }
+    return T1(stan::math::pow(x1.val_, x2),
+              x1.d_ * x2 * stan::math::pow(x1.val_, x2 - 1));
   }
-  if (x2 == 2.0) {
-    return square(x1);
-  }
-  return fvar<T>(pow(x1.val_, x2), x1.d_ * x2 * pow(x1.val_, x2 - 1));
-}
-
-// must uniquely match all pairs of:
-//    { complex<fvar<V>>, complex<T>, fvar<V>, T }
-// with at least one fvar<V> and at least one complex, where T is arithmetic:
-// 1) complex<fvar<V>>, complex<fvar<V>>
-// 2) complex<fvar<V>>, complex<T>
-// 3) complex<fvar<V>>, fvar<V>
-// 4) complex<fvar<V>>, T
-// 5) complex<T>, complex<fvar<V>>
-// 6) complex<T>, fvar<V>
-// 7) fvar<V>, complex<fvar<V>>
-// 8) fvar<V>, complex<T>
-// 9) T, complex<fvar<V>>
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V>
-inline std::complex<fvar<V>> pow(const std::complex<fvar<V>>& x,
-                                 const std::complex<fvar<V>>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V, typename T, typename = require_arithmetic_t<T>>
-inline std::complex<fvar<V>> pow(const std::complex<fvar<V>>& x,
-                                 const std::complex<T>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V>
-inline std::complex<fvar<V>> pow(const std::complex<fvar<V>>& x,
-                                 const fvar<V>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V, typename T, typename = require_arithmetic_t<T>>
-inline std::complex<fvar<V>> pow(const std::complex<fvar<V>>& x, const T& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V, typename T, typename = require_arithmetic_t<T>>
-inline std::complex<fvar<V>> pow(const std::complex<T>& x,
-                                 const std::complex<fvar<V>>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V, typename T, typename = require_arithmetic_t<T>>
-inline std::complex<fvar<V>> pow(const std::complex<T>& x, const fvar<V>& y) {
-  return internal::complex_pow(x, y);
 }
 
 /**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V>
-inline std::complex<fvar<V>> pow(const fvar<V>& x,
-                                 const std::complex<fvar<V>>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename V, typename T, typename = require_arithmetic_t<T>>
-inline std::complex<fvar<V>> pow(const fvar<V>& x, const std::complex<T>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam V autodiff value type
- * @tparam T real type (`fvar<V>` or arithmetic)
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename V, typename = require_arithmetic_t<T>>
-inline std::complex<fvar<V>> pow(T x, const std::complex<fvar<V>>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * Note: this overload is required because gcc still provides the
- * C++99 template function `pow(complex<T>, int)`, which introduces
- * an ambiguity.
+ * Returns the elementwise raising of the first argument to the power of the
+ * second argument.
  *
- * @tparam T autodiff value type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
+ * @tparam T1 type of first argument
+ * @tparam T2 type of second argument
+ * @param a first argument
+ * @param b second argument
+ * @return the elementwise raising of the first argument to the power of the
+ * second argument.
  */
-template <typename T>
-inline std::complex<fvar<T>> pow(const std::complex<fvar<T>>& x, int y) {
-  return internal::complex_pow(x, y);
+template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr,
+          require_all_not_matrix_st<is_var, T1, T2>* = nullptr,
+          require_any_fvar_t<base_type_t<T1>, base_type_t<T2>>* = nullptr>
+inline auto pow(const T1& a, const T2& b) {
+  return apply_scalar_binary(
+      a, b, [](const auto& c, const auto& d) { return stan::math::pow(c, d); });
 }
 
 }  // namespace math
diff --git a/stan/math/fwd/fun/sum.hpp b/stan/math/fwd/fun/sum.hpp
index 2ae7887c1ca..36eef6ad687 100644
--- a/stan/math/fwd/fun/sum.hpp
+++ b/stan/math/fwd/fun/sum.hpp
@@ -1,10 +1,11 @@
 #ifndef STAN_MATH_FWD_FUN_SUM_HPP
 #define STAN_MATH_FWD_FUN_SUM_HPP
 
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/sum.hpp>
-#include <stan/math/fwd/core.hpp>
 #include <vector>
 
 namespace stan {
@@ -18,18 +19,18 @@ namespace math {
  * @param m Vector.
  * @return Sum of vector entries.
  */
-template <typename T>
-inline fvar<T> sum(const std::vector<fvar<T>>& m) {
+template <typename T, require_fvar_t<T>* = nullptr>
+inline auto sum(const std::vector<T>& m) {
   if (m.size() == 0) {
-    return 0.0;
+    return T(0.0);
   }
-  std::vector<T> vals(m.size());
-  std::vector<T> tans(m.size());
+  std::vector<partials_type_t<T>> vals(m.size());
+  std::vector<partials_type_t<T>> tans(m.size());
   for (size_t i = 0; i < m.size(); ++i) {
     vals[i] = m[i].val();
     tans[i] = m[i].d();
   }
-  return fvar<T>(sum(vals), sum(tans));
+  return T(sum(vals), sum(tans));
 }
 
 /**
diff --git a/stan/math/mix.hpp b/stan/math/mix.hpp
index bb7ff7c0b1f..876916443ce 100644
--- a/stan/math/mix.hpp
+++ b/stan/math/mix.hpp
@@ -5,6 +5,13 @@
 #include <stan/math/mix/fun.hpp>
 #include <stan/math/mix/functor.hpp>
 
+#include <stan/math/rev/constraint.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun.hpp>
+#include <stan/math/rev/functor.hpp>
+#include <stan/math/rev/prob.hpp>
+
 #include <stan/math/fwd/constraint.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/meta.hpp>
@@ -17,13 +24,6 @@
 #include <stan/math/opencl/rev_constraint.hpp>
 #endif
 
-#include <stan/math/rev/constraint.hpp>
-#include <stan/math/rev/core.hpp>
-#include <stan/math/rev/meta.hpp>
-#include <stan/math/rev/fun.hpp>
-#include <stan/math/rev/functor.hpp>
-#include <stan/math/rev/prob.hpp>
-
 #include <stan/math/prim.hpp>
 
 #endif
diff --git a/stan/math/prim/fun/pow.hpp b/stan/math/prim/fun/pow.hpp
index 2ed54148042..95809bd8830 100644
--- a/stan/math/prim/fun/pow.hpp
+++ b/stan/math/prim/fun/pow.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/functor/apply_scalar_binary.hpp>
 #include <cmath>
 #include <complex>
@@ -39,14 +40,29 @@ inline complex_return_t<U, V> complex_pow(const U& x, const V& y) {
  * @return the first argument raised to the power of the second
  * argument.
  */
-template <typename T1, typename T2,
-          require_all_t<
-              disjunction<is_complex<T1>, std::is_arithmetic<T1>>,
-              disjunction<is_complex<T2>, std::is_arithmetic<T2>>>* = nullptr>
-inline auto pow(const T1& a, const T2& b) {
+template <typename T1, typename T2, require_arithmetic_t<T1>* = nullptr,
+          require_arithmetic_t<T2>* = nullptr>
+inline auto pow(const std::complex<T1>& a, const std::complex<T2>& b) {
+  return std::pow(a, b);
+}
+
+template <typename T1, typename T2, require_arithmetic_t<T1>* = nullptr,
+          require_arithmetic_t<T2>* = nullptr>
+inline auto pow(const T1& a, const std::complex<T2>& b) {
+  return std::pow(a, b);
+}
+
+template <typename T1, typename T2, require_arithmetic_t<T1>* = nullptr,
+          require_arithmetic_t<T2>* = nullptr>
+inline auto pow(const std::complex<T1>& a, const T2& b) {
   return std::pow(a, b);
 }
 
+template <typename T1, typename T2, require_arithmetic_t<T1>* = nullptr,
+          require_arithmetic_t<T2>* = nullptr>
+inline auto pow(const T1& a, const T2& b) {
+  return std::pow(a, b);
+}
 /**
  * Returns the elementwise raising of the first argument to the power of the
  * second argument.
@@ -59,12 +75,11 @@ inline auto pow(const T1& a, const T2& b) {
  * second argument.
  */
 template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr,
-          require_all_not_matrix_st<is_var, T1, T2>* = nullptr>
+          require_all_not_matrix_st<is_var, T1, T2>* = nullptr,
+          require_all_arithmetic_t<base_type_t<T1>, base_type_t<T2>>* = nullptr>
 inline auto pow(const T1& a, const T2& b) {
-  return apply_scalar_binary(a, b, [](const auto& c, const auto& d) {
-    using std::pow;
-    return pow(c, d);
-  });
+  return apply_scalar_binary(
+      a, b, [](const auto& c, const auto& d) { return stan::math::pow(c, d); });
 }
 }  // namespace math
 }  // namespace stan
diff --git a/stan/math/prim/fun/sum.hpp b/stan/math/prim/fun/sum.hpp
index f1256c375e4..0440997fda0 100644
--- a/stan/math/prim/fun/sum.hpp
+++ b/stan/math/prim/fun/sum.hpp
@@ -29,7 +29,7 @@ inline T sum(T&& m) {
  * @param m Standard vector to sum.
  * @return Sum of elements.
  */
-template <typename T, require_not_var_t<T>* = nullptr>
+template <typename T, require_not_autodiff_t<T>* = nullptr>
 inline T sum(const std::vector<T>& m) {
   return std::accumulate(m.begin(), m.end(), T{0});
 }
diff --git a/stan/math/rev/fun/pow.hpp b/stan/math/rev/fun/pow.hpp
index 8a2383880b9..669f53908b4 100644
--- a/stan/math/rev/fun/pow.hpp
+++ b/stan/math/rev/fun/pow.hpp
@@ -61,45 +61,54 @@ namespace math {
    \end{cases}
    \f]
  *
+ * @tparam Scal1 Either a `var`, `arithmetic`, or `complex` type with an inner
+ `var` or `arithmetic` type.
+ * @tparam Scal2 Either a `var`, `arithmetic`, or `complex` type with an inner
+ `var` or `arithmetic` type.
  * @param base Base variable.
  * @param exponent Exponent variable.
  * @return Base raised to the exponent.
  */
 template <typename Scal1, typename Scal2,
-          require_any_st_var<Scal1, Scal2>* = nullptr,
+          require_any_var_t<base_type_t<Scal1>, base_type_t<Scal2>>* = nullptr,
           require_all_stan_scalar_t<Scal1, Scal2>* = nullptr>
-inline var pow(const Scal1& base, const Scal2& exponent) {
-  if (is_constant<Scal2>::value) {
-    if (exponent == 0.5) {
-      return sqrt(base);
-    } else if (exponent == 1.0) {
-      return base;
-    } else if (exponent == 2.0) {
-      return square(base);
-    } else if (exponent == -2.0) {
-      return inv_square(base);
-    } else if (exponent == -1.0) {
-      return inv(base);
-    } else if (exponent == -0.5) {
-      return inv_sqrt(base);
+inline auto pow(const Scal1& base, const Scal2& exponent) {
+  if constexpr (is_complex<Scal1>::value || is_complex<Scal2>::value) {
+    return internal::complex_pow(base, exponent);
+  } else {
+    if constexpr (is_constant<Scal2>::value) {
+      if (exponent == 0.5) {
+        return sqrt(base);
+      } else if (exponent == 1.0) {
+        return base;
+      } else if (exponent == 2.0) {
+        return square(base);
+      } else if (exponent == -2.0) {
+        return inv_square(base);
+      } else if (exponent == -1.0) {
+        return inv(base);
+      } else if (exponent == -0.5) {
+        return inv_sqrt(base);
+      }
     }
-  }
-  return make_callback_var(
-      std::pow(value_of(base), value_of(exponent)),
-      [base, exponent](auto&& vi) mutable {
-        if (value_of(base) == 0.0) {
-          return;  // partials zero, avoids 0 & log(0)
-        }
-        const double vi_mul = vi.adj() * vi.val();
+    return make_callback_var(std::pow(value_of(base), value_of(exponent)),
+                             [base, exponent](auto&& vi) mutable {
+                               if (value_of(base) == 0.0) {
+                                 return;  // partials zero, avoids 0 & log(0)
+                               }
+                               const double vi_mul = vi.adj() * vi.val();
 
-        if (!is_constant<Scal1>::value) {
-          forward_as<var>(base).adj()
-              += vi_mul * value_of(exponent) / value_of(base);
-        }
-        if (!is_constant<Scal2>::value) {
-          forward_as<var>(exponent).adj() += vi_mul * std::log(value_of(base));
-        }
-      });
+                               if (!is_constant<Scal1>::value) {
+                                 forward_as<var>(base).adj()
+                                     += vi_mul * value_of(exponent)
+                                        / value_of(base);
+                               }
+                               if (!is_constant<Scal2>::value) {
+                                 forward_as<var>(exponent).adj()
+                                     += vi_mul * std::log(value_of(base));
+                               }
+                             });
+  }
 }
 
 /**
@@ -164,6 +173,7 @@ inline auto pow(const Mat1& base, const Mat2& exponent) {
  * @tparam Mat1 An Eigen type deriving from Eigen::EigenBase or
  *  a `var_value` with inner Eigen type as defined above. The `scalar_type`
  *  must be a `var` or Arithmetic.
+ * @tparam Scal1 An arithmetic type or a `var_value` with inner arithmetic type.
  * @param base Base variable.
  * @param exponent Exponent variable.
  * @return Base raised to the exponent.
@@ -225,10 +235,10 @@ inline auto pow(const Mat1& base, const Scal1& exponent) {
  * \f$\frac{d}{d y} \mbox{pow}(c, y) = c^y \log c \f$.
  *
  *
- * @tparam Mat An Eigen type deriving from Eigen::EigenBase or
+ * @tparam Mat1 An Eigen type deriving from Eigen::EigenBase or
  *  a `var_value` with inner Eigen type as defined above. The `scalar_type`
  * must be a `var`.
- *
+ * @tparam Scal1 An arithmetic type or a `var_value` with inner arithmetic type.
  * @param base Base scalar.
  * @param exponent Exponent variable.
  * @return Base raised to the exponent.
@@ -261,144 +271,23 @@ inline auto pow(Scal1 base, const Mat1& exponent) {
   return ret_type(ret);
 }
 
-// must uniquely match all pairs of { complex<var>, complex<T>, var, T }
-// with at least one var and at least one complex, where T is arithmetic:
-// 1) complex<var>, complex<var>
-// 2) complex<var>, complex<T>
-// 3) complex<var>, var
-// 4) complex<var>, T
-// 5) complex<T>, complex<var>
-// 6) complex<T>, var
-// 7) var, complex<var>
-// 8) var, complex<T>
-// 9) T, complex<var>
-
 /**
- * Return the first argument raised to the power of the second argument.
- *
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-inline std::complex<var> pow(const std::complex<var>& x,
-                             const std::complex<var>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename = require_arithmetic_t<T>>
-inline std::complex<var> pow(const std::complex<var>& x,
-                             const std::complex<T> y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-inline std::complex<var> pow(const std::complex<var>& x, const var& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename = require_arithmetic_t<T>>
-inline std::complex<var> pow(const std::complex<var>& x, T y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename = require_arithmetic_t<T>>
-inline std::complex<var> pow(std::complex<T> x, const std::complex<var>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename = require_arithmetic_t<T>>
-inline std::complex<var> pow(std::complex<T> x, const var& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-inline std::complex<var> pow(const var& x, const std::complex<var>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename = require_arithmetic_t<T>>
-inline std::complex<var> pow(const var& x, std::complex<T> y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * @tparam T arithmetic type
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
- */
-template <typename T, typename = require_arithmetic_t<T>>
-inline std::complex<var> pow(T x, const std::complex<var>& y) {
-  return internal::complex_pow(x, y);
-}
-
-/**
- * Return the first argument raised to the power of the second argument.
- *
- * Note: this overload is required because gcc still provides the
- * C++99 template function `pow(complex<T>, int)`, which introduces
- * an ambiguity.
+ * Returns the elementwise raising of the first argument to the power of the
+ * second argument.
  *
- * @param x first argument
- * @param y second argument
- * @return first argument to the power of the second argument
+ * @tparam T1 type of first argument
+ * @tparam T2 type of second argument
+ * @param a first argument
+ * @param b second argument
+ * @return the elementwise raising of the first argument to the power of the
+ * second argument.
  */
-inline std::complex<var> pow(const std::complex<var>& x, int y) {
-  return internal::complex_pow(x, y);
+template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr,
+          require_all_not_matrix_st<is_var, T1, T2>* = nullptr,
+          require_any_var_t<base_type_t<T1>, base_type_t<T2>>* = nullptr>
+inline auto pow(const T1& a, const T2& b) {
+  return apply_scalar_binary(
+      a, b, [](const auto& c, const auto& d) { return stan::math::pow(c, d); });
 }
 
 }  // namespace math
diff --git a/test/code_generator.py b/test/code_generator.py
index 2671b72429d..895267a7742 100644
--- a/test/code_generator.py
+++ b/test/code_generator.py
@@ -75,7 +75,7 @@ def build_arguments(self, signature_parser, arg_overloads, size):
 
             # The first case here is used for the array initializers in sig_utils.special_arg_values
             # Everything else uses the second case
-            if number_nested_arrays > 0 and isinstance(value, collections.Iterable):
+            if number_nested_arrays > 0 and isinstance(value, collections.abc.Iterable):
                 arg = statement_types.ArrayVariable(
                     overload,
                     "array" + suffix,
diff --git a/test/unit/math/mix/fun/pow_part1_test.cpp b/test/unit/math/mix/fun/pow_part1_test.cpp
index f97b645a8ec..72100764562 100644
--- a/test/unit/math/mix/fun/pow_part1_test.cpp
+++ b/test/unit/math/mix/fun/pow_part1_test.cpp
@@ -5,20 +5,19 @@
 
 template <typename T>
 void expect_arith_instantiate() {
-  using stan::math::pow;
-  auto a1 = pow(T(1.0), 1);
-  auto b1 = pow(T(1.0), 1.0);
-  auto c1 = pow(1, T(1.0));
-  auto d1 = pow(1.0, T(1.0));
-  auto e1 = pow(T(1.0), T(1.0));
-
-  auto a2 = pow(std::complex<T>(1.0), 1);
-  auto b2 = pow(std::complex<T>(1.0), 1.0);
-  auto c2 = pow(1, std::complex<T>(1.0));
-  auto d2 = pow(1.0, std::complex<T>(1.0));
-  auto e2 = pow(std::complex<T>(1.0), std::complex<T>(1.0));
-  auto f2 = pow(std::complex<double>(1.0), std::complex<T>(1.0));
-  auto g2 = pow(std::complex<T>(1.0), std::complex<double>(1.0));
+  auto a1 = stan::math::pow(T(1.0), 1);
+  auto b1 = stan::math::pow(T(1.0), 1.0);
+  auto c1 = stan::math::pow(1, T(1.0));
+  auto d1 = stan::math::pow(1.0, T(1.0));
+  auto e1 = stan::math::pow(T(1.0), T(1.0));
+
+  auto a2 = stan::math::pow(std::complex<T>(1.0), 1);
+  auto b2 = stan::math::pow(std::complex<T>(1.0), 1.0);
+  auto c2 = stan::math::pow(1, std::complex<T>(1.0));
+  auto d2 = stan::math::pow(1.0, std::complex<T>(1.0));
+  auto e2 = stan::math::pow(std::complex<T>(1.0), std::complex<T>(1.0));
+  auto f2 = stan::math::pow(std::complex<double>(1.0), std::complex<T>(1.0));
+  auto g2 = stan::math::pow(std::complex<T>(1.0), std::complex<double>(1.0));
 }
 
 // this one's been tricky to instantiate, so test all instances
@@ -34,10 +33,8 @@ TEST(mathMixScalFun, powInstantiations) {
 }
 
 TEST(mathMixScalFun, pow) {
-  auto f = [](const auto& x1, const auto& x2) {
-    using stan::math::pow;
-    return pow(x1, x2);
-  };
+  auto f
+      = [](const auto& x1, const auto& x2) { return stan::math::pow(x1, x2); };
 
   stan::test::expect_ad(f, -0.4, 0.5);
   stan::test::expect_ad(f, 0.5, 0.5);
diff --git a/test/unit/math/mix/fun/pow_part2_test.cpp b/test/unit/math/mix/fun/pow_part2_test.cpp
index 6fb46d08e2d..d3805f61cba 100644
--- a/test/unit/math/mix/fun/pow_part2_test.cpp
+++ b/test/unit/math/mix/fun/pow_part2_test.cpp
@@ -4,10 +4,8 @@
 #include <vector>
 
 TEST(mathMixFun, complexPow) {
-  auto f = [](const auto& x1, const auto& x2) {
-    using stan::math::pow;
-    return pow(x1, x2);
-  };
+  auto f
+      = [](const auto& x1, const auto& x2) { return stan::math::pow(x1, x2); };
   stan::test::ad_tolerances tols;
   tols.hessian_hessian_ = 5e-3;
   tols.hessian_fvar_hessian_ = 5e-3;
@@ -49,7 +47,6 @@ TEST(mathMixFun, complexPow) {
 }
 
 TEST(mathMixFun, powIntAmbiguityTest) {
-  using stan::math::pow;  // included to check ambiguities
   using stan::math::var;
   using std::complex;
   int i = 2;
@@ -58,35 +55,35 @@ TEST(mathMixFun, powIntAmbiguityTest) {
   complex<double> cd = 2.5;
   complex<var> cv = 2.5;
 
-  auto a1 = pow(i, i);
-  auto a2 = pow(i, d);
-  auto a3 = pow(i, v);
-  auto a4 = pow(i, cd);
-  auto a5 = pow(i, cv);
-
-  auto b1 = pow(d, i);
-  auto b2 = pow(d, d);
-  auto b3 = pow(d, v);
-  auto b4 = pow(d, cd);
-  auto b5 = pow(d, cv);
-
-  auto e1 = pow(v, i);
-  auto e2 = pow(v, d);
-  auto e3 = pow(v, v);
-  auto e4 = pow(v, cd);
-  auto e5 = pow(v, cv);
-
-  auto c1 = pow(cd, i);
-  auto c2 = pow(cd, d);
-  auto c3 = pow(cd, v);
-  auto c4 = pow(cd, cd);
-  auto c5 = pow(cd, cv);
-
-  auto d1 = pow(cv, i);
-  auto d2 = pow(cv, d);
-  auto d3 = pow(cv, v);
-  auto d4 = pow(cv, cd);
-  auto d5 = pow(cv, cv);
+  auto a1 = stan::math::pow(i, i);
+  auto a2 = stan::math::pow(i, d);
+  auto a3 = stan::math::pow(i, v);
+  auto a4 = stan::math::pow(i, cd);
+  auto a5 = stan::math::pow(i, cv);
+
+  auto b1 = stan::math::pow(d, i);
+  auto b2 = stan::math::pow(d, d);
+  auto b3 = stan::math::pow(d, v);
+  auto b4 = stan::math::pow(d, cd);
+  auto b5 = stan::math::pow(d, cv);
+
+  auto e1 = stan::math::pow(v, i);
+  auto e2 = stan::math::pow(v, d);
+  auto e3 = stan::math::pow(v, v);
+  auto e4 = stan::math::pow(v, cd);
+  auto e5 = stan::math::pow(v, cv);
+
+  auto c1 = stan::math::pow(cd, i);
+  auto c2 = stan::math::pow(cd, d);
+  auto c3 = stan::math::pow(cd, v);
+  auto c4 = stan::math::pow(cd, cd);
+  auto c5 = stan::math::pow(cd, cv);
+
+  auto d1 = stan::math::pow(cv, i);
+  auto d2 = stan::math::pow(cv, d);
+  auto d3 = stan::math::pow(cv, v);
+  auto d4 = stan::math::pow(cv, cd);
+  auto d5 = stan::math::pow(cv, cv);
 
   auto e = a1 + a2 + a3 + a4 + a5 + b1 + b2 + b3 + b4 + b5 + c1 + c2 + c3 + c4
            + c5 + d1 + d2 + d3 + d4 + d5 + e1 + e2 + e3 + e4 + e5;
@@ -96,7 +93,6 @@ TEST(mathMixFun, powIntAmbiguityTest) {
 
 TEST(mathMixFun, powIntAmbiguityTestFvar) {
   using stan::math::fvar;
-  using stan::math::pow;  // included to check ambiguities
   using std::complex;
   int i = 2;
   double d = 2.5;
@@ -104,29 +100,29 @@ TEST(mathMixFun, powIntAmbiguityTestFvar) {
   complex<double> cd = 2.5;
   complex<fvar<double>> cv = 2.5;
 
-  auto a1 = pow(i, i);
-  auto a2 = pow(i, d);
-  auto a3 = pow(i, v);
-  auto a4 = pow(i, cd);
-  auto a5 = pow(i, cv);
-
-  auto b1 = pow(d, i);
-  auto b2 = pow(d, d);
-  auto b3 = pow(d, v);
-  auto b4 = pow(d, cd);
-  auto b5 = pow(d, cv);
-
-  auto c1 = pow(cd, i);
-  auto c2 = pow(cd, d);
-  auto c3 = pow(cd, v);
-  auto c4 = pow(cd, cd);
-  auto c5 = pow(cd, cv);
-
-  auto d1 = pow(cv, i);
-  auto d2 = pow(cv, d);
-  auto d3 = pow(cv, v);
-  auto d4 = pow(cv, cd);
-  auto d5 = pow(cv, cv);
+  auto a1 = stan::math::pow(i, i);
+  auto a2 = stan::math::pow(i, d);
+  auto a3 = stan::math::pow(i, v);
+  auto a4 = stan::math::pow(i, cd);
+  auto a5 = stan::math::pow(i, cv);
+
+  auto b1 = stan::math::pow(d, i);
+  auto b2 = stan::math::pow(d, d);
+  auto b3 = stan::math::pow(d, v);
+  auto b4 = stan::math::pow(d, cd);
+  auto b5 = stan::math::pow(d, cv);
+
+  auto c1 = stan::math::pow(cd, i);
+  auto c2 = stan::math::pow(cd, d);
+  auto c3 = stan::math::pow(cd, v);
+  auto c4 = stan::math::pow(cd, cd);
+  auto c5 = stan::math::pow(cd, cv);
+
+  auto d1 = stan::math::pow(cv, i);
+  auto d2 = stan::math::pow(cv, d);
+  auto d3 = stan::math::pow(cv, v);
+  auto d4 = stan::math::pow(cv, cd);
+  auto d5 = stan::math::pow(cv, cv);
 
   auto e = a1 + a2 + a3 + a4 + a5 + b1 + b2 + b3 + b4 + b5 + c1 + c2 + c3 + c4
            + c5 + d1 + d2 + d3 + d4 + d5;