diff --git a/stan/math/fwd/core/operator_division.hpp b/stan/math/fwd/core/operator_division.hpp
index 6ed5bdb168e..93514e0a642 100644
--- a/stan/math/fwd/core/operator_division.hpp
+++ b/stan/math/fwd/core/operator_division.hpp
@@ -2,6 +2,9 @@
 #define STAN_MATH_FWD_CORE_OPERATOR_DIVISION_HPP
 
 #include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/prim/core/operator_division.hpp>
+#include <complex>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -28,8 +31,8 @@ inline fvar<T> operator/(const fvar<T>& x1, const fvar<T>& x2) {
  * @param x2 second argument
  * @return first argument divided by second argument
  */
-template <typename T>
-inline fvar<T> operator/(const fvar<T>& x1, double x2) {
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline fvar<T> operator/(const fvar<T>& x1, U x2) {
   return fvar<T>(x1.val_ / x2, x1.d_ / x2);
 }
 
@@ -41,11 +44,59 @@ inline fvar<T> operator/(const fvar<T>& x1, double x2) {
  * @param x2 second argument
  * @return first argument divided by second argument
  */
-template <typename T>
-inline fvar<T> operator/(double x1, const fvar<T>& x2) {
-  // TODO(carpenter): store x1 / x2.val_ and reuse
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline fvar<T> operator/(U x1, const fvar<T>& x2) {
   return fvar<T>(x1 / x2.val_, -x1 * x2.d_ / (x2.val_ * x2.val_));
 }
+
+template <typename T>
+inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1,
+                                       const std::complex<fvar<T>>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1,
+                                       const std::complex<U>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+template <typename T>
+inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1,
+                                       const fvar<T>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline std::complex<fvar<T>> operator/(const std::complex<fvar<T>>& x1, U x2) {
+  return internal::complex_divide(x1, x2);
+}
+
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline std::complex<fvar<T>> operator/(const std::complex<U>& x1,
+                                       const std::complex<fvar<T>>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline std::complex<fvar<T>> operator/(const std::complex<U>& x1,
+                                       const fvar<T>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+
+template <typename T>
+inline std::complex<fvar<T>> operator/(const fvar<T>& x1,
+                                       const std::complex<fvar<T>>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+template <typename T, typename U,
+          typename = std::enable_if_t<std::is_arithmetic<U>::value>>
+inline std::complex<fvar<T>> operator/(const fvar<T>& x1,
+                                       const std::complex<U>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+
+template <typename T, typename U, require_arithmetic_t<U>...>
+inline std::complex<fvar<T>> operator/(U x1, const std::complex<fvar<T>>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/core/std_complex.hpp b/stan/math/fwd/core/std_complex.hpp
index 6492a11f7e5..49c8cda0448 100644
--- a/stan/math/fwd/core/std_complex.hpp
+++ b/stan/math/fwd/core/std_complex.hpp
@@ -32,8 +32,12 @@ class complex<stan::math::fvar<T>>
    * @tparam Scalar real type (must be assignable to `value_type`)
    * @param[in] re real part
    */
-  template <typename U, typename = stan::require_stan_scalar_t<U>>
-  complex(U&& re) : base_t(re) {}  // NOLINT(runtime/explicit)
+  template <typename U>  // , typename = stan::require_stan_scalar_t<U>>
+  complex(const U& re) : base_t(re) {}  // NOLINT(runtime/explicit)
+
+  template <typename U>
+  complex(const std::complex<U>& z)  // NOLINT(runtime/explicit)
+      : base_t(z.real(), z.imag()) {}
 
   /**
    * Construct a complex number from the specified real and imaginary
diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp
index 7da34e5e62d..2e3c9c6c2aa 100644
--- a/stan/math/fwd/fun.hpp
+++ b/stan/math/fwd/fun.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/fwd/fun/acos.hpp>
 #include <stan/math/fwd/fun/acosh.hpp>
 #include <stan/math/fwd/fun/asin.hpp>
+#include <stan/math/fwd/fun/arg.hpp>
 #include <stan/math/fwd/fun/asinh.hpp>
 #include <stan/math/fwd/fun/atan.hpp>
 #include <stan/math/fwd/fun/atan2.hpp>
@@ -15,6 +16,7 @@
 #include <stan/math/fwd/fun/binary_log_loss.hpp>
 #include <stan/math/fwd/fun/cbrt.hpp>
 #include <stan/math/fwd/fun/ceil.hpp>
+#include <stan/math/fwd/fun/conj.hpp>
 #include <stan/math/fwd/fun/cos.hpp>
 #include <stan/math/fwd/fun/cosh.hpp>
 #include <stan/math/fwd/fun/determinant.hpp>
@@ -79,11 +81,14 @@
 #include <stan/math/fwd/fun/multiply.hpp>
 #include <stan/math/fwd/fun/multiply_log.hpp>
 #include <stan/math/fwd/fun/multiply_lower_tri_self_transpose.hpp>
+#include <stan/math/fwd/fun/norm.hpp>
 #include <stan/math/fwd/fun/owens_t.hpp>
 #include <stan/math/fwd/fun/Phi.hpp>
 #include <stan/math/fwd/fun/Phi_approx.hpp>
+#include <stan/math/fwd/fun/polar.hpp>
 #include <stan/math/fwd/fun/pow.hpp>
 #include <stan/math/fwd/fun/primitive_value.hpp>
+#include <stan/math/fwd/fun/proj.hpp>
 #include <stan/math/fwd/fun/quad_form.hpp>
 #include <stan/math/fwd/fun/quad_form_sym.hpp>
 #include <stan/math/fwd/fun/rising_factorial.hpp>
diff --git a/stan/math/fwd/fun/abs.hpp b/stan/math/fwd/fun/abs.hpp
index f52c49482da..ffd21791ff0 100644
--- a/stan/math/fwd/fun/abs.hpp
+++ b/stan/math/fwd/fun/abs.hpp
@@ -3,9 +3,10 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
-#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/abs.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -23,6 +24,18 @@ inline fvar<T> abs(const fvar<T>& x) {
   }
 }
 
+/**
+ * Return the absolute value of the complex argument.
+ *
+ * @tparam T value type of argument
+ * @param[in] z argument
+ * @return absolute value of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> abs(const std::complex<fvar<T>>& z) {
+  return internal::complex_abs(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/acos.hpp b/stan/math/fwd/fun/acos.hpp
index b6f9bfbd683..a3cf5315686 100644
--- a/stan/math/fwd/fun/acos.hpp
+++ b/stan/math/fwd/fun/acos.hpp
@@ -1,10 +1,12 @@
 #ifndef STAN_MATH_FWD_FUN_ACOS_HPP
 #define STAN_MATH_FWD_FUN_ACOS_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
 #include <stan/math/prim/fun/square.hpp>
+#include <stan/math/prim/fun/acos.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -15,6 +17,19 @@ inline fvar<T> acos(const fvar<T>& x) {
   using std::sqrt;
   return fvar<T>(acos(x.val_), x.d_ / -sqrt(1 - square(x.val_)));
 }
+
+/**
+ * Return the arc cosine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return arc cosine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> acos(const std::complex<fvar<T>>& z) {
+  return internal::complex_acos(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/acosh.hpp b/stan/math/fwd/fun/acosh.hpp
index 948e7fe226c..376daa568aa 100644
--- a/stan/math/fwd/fun/acosh.hpp
+++ b/stan/math/fwd/fun/acosh.hpp
@@ -5,7 +5,9 @@
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/prim/fun/acosh.hpp>
 #include <stan/math/prim/fun/square.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -16,6 +18,18 @@ inline fvar<T> acosh(const fvar<T>& x) {
   return fvar<T>(acosh(x.val_), x.d_ / sqrt(square(x.val_) - 1));
 }
 
+/**
+ * Return the hyperbolic arc cosine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return hyperbolic arc cosine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> acosh(const std::complex<fvar<T>>& z) {
+  return internal::complex_acosh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/arg.hpp b/stan/math/fwd/fun/arg.hpp
new file mode 100644
index 00000000000..89e3ae192ba
--- /dev/null
+++ b/stan/math/fwd/fun/arg.hpp
@@ -0,0 +1,26 @@
+#ifndef STAN_MATH_FWD_FUN_ARG_HPP
+#define STAN_MATH_FWD_FUN_ARG_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/arg.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the phase angle of the complex argument.
+ *
+ * @tparam T value type of autodiff variable
+ * @param[in] z argument
+ * @return phase angle of the argument
+ */
+template <typename T>
+inline fvar<T> arg(const std::complex<fvar<T>>& z) {
+  return internal::complex_arg(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/fwd/fun/asin.hpp b/stan/math/fwd/fun/asin.hpp
index a4d4e4de0c8..0ce1b593342 100644
--- a/stan/math/fwd/fun/asin.hpp
+++ b/stan/math/fwd/fun/asin.hpp
@@ -3,8 +3,10 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/asin.hpp>
 #include <stan/math/prim/fun/square.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -16,6 +18,18 @@ inline fvar<T> asin(const fvar<T>& x) {
   return fvar<T>(asin(x.val_), x.d_ / sqrt(1 - square(x.val_)));
 }
 
+/**
+ * Return the arc sine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return arc sine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> asin(const std::complex<fvar<T>>& z) {
+  return internal::complex_asin(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/asinh.hpp b/stan/math/fwd/fun/asinh.hpp
index 7a4b0df5693..4e430514ec4 100644
--- a/stan/math/fwd/fun/asinh.hpp
+++ b/stan/math/fwd/fun/asinh.hpp
@@ -1,11 +1,12 @@
 #ifndef STAN_MATH_FWD_FUN_ASINH_HPP
 #define STAN_MATH_FWD_FUN_ASINH_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
 #include <stan/math/prim/fun/asinh.hpp>
 #include <stan/math/prim/fun/square.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -16,6 +17,18 @@ inline fvar<T> asinh(const fvar<T>& x) {
   return fvar<T>(asinh(x.val_), x.d_ / sqrt(square(x.val_) + 1));
 }
 
+/**
+ * Return the hyperbolic arcsine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return hyperbolic arcsine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> asinh(const std::complex<fvar<T>>& z) {
+  return internal::complex_asinh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/atan.hpp b/stan/math/fwd/fun/atan.hpp
index 5ded2b132f4..a16f6ef8792 100644
--- a/stan/math/fwd/fun/atan.hpp
+++ b/stan/math/fwd/fun/atan.hpp
@@ -1,10 +1,12 @@
 #ifndef STAN_MATH_FWD_FUN_ATAN_HPP
 #define STAN_MATH_FWD_FUN_ATAN_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/fun/atan.hpp>
 #include <stan/math/prim/fun/square.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -15,6 +17,18 @@ inline fvar<T> atan(const fvar<T>& x) {
   return fvar<T>(atan(x.val_), x.d_ / (1 + square(x.val_)));
 }
 
+/**
+ * Return the arc tangent of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return arc tanget of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> atan(const std::complex<fvar<T>>& z) {
+  return internal::complex_atan(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/atanh.hpp b/stan/math/fwd/fun/atanh.hpp
index bd1154d281c..ded211fd2d4 100644
--- a/stan/math/fwd/fun/atanh.hpp
+++ b/stan/math/fwd/fun/atanh.hpp
@@ -1,10 +1,12 @@
 #ifndef STAN_MATH_FWD_FUN_ATANH_HPP
 #define STAN_MATH_FWD_FUN_ATANH_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
 #include <stan/math/prim/fun/atanh.hpp>
 #include <stan/math/prim/fun/square.hpp>
+#include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -23,6 +25,18 @@ inline fvar<T> atanh(const fvar<T>& x) {
   return fvar<T>(atanh(x.val_), x.d_ / (1 - square(x.val_)));
 }
 
+/**
+ * Return the hyperbolic arc tangent of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return hyperbolic arc tangent of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> atanh(const std::complex<fvar<T>>& z) {
+  return internal::complex_atanh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/conj.hpp b/stan/math/fwd/fun/conj.hpp
new file mode 100644
index 00000000000..706d97d2678
--- /dev/null
+++ b/stan/math/fwd/fun/conj.hpp
@@ -0,0 +1,26 @@
+#ifndef STAN_MATH_FWD_FUN_CONJ_HPP
+#define STAN_MATH_FWD_FUN_CONJ_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/conj.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the phase angle of the complex argument.
+ *
+ * @tparam T value type of autodiff variable
+ * @param[in] z argument
+ * @return phase angle of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> conj(const std::complex<fvar<T>>& z) {
+  return internal::complex_conj(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/fwd/fun/cos.hpp b/stan/math/fwd/fun/cos.hpp
index de589753db3..0cfe38166b8 100644
--- a/stan/math/fwd/fun/cos.hpp
+++ b/stan/math/fwd/fun/cos.hpp
@@ -1,9 +1,11 @@
 #ifndef STAN_MATH_FWD_FUN_COS_HPP
 #define STAN_MATH_FWD_FUN_COS_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/fun/cos.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -15,6 +17,18 @@ inline fvar<T> cos(const fvar<T>& x) {
   return fvar<T>(cos(x.val_), x.d_ * -sin(x.val_));
 }
 
+/**
+ * Return the cosine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return cosine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> cos(const std::complex<fvar<T>>& z) {
+  return internal::complex_cos(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/cosh.hpp b/stan/math/fwd/fun/cosh.hpp
index 124063b767e..d1194cf54ed 100644
--- a/stan/math/fwd/fun/cosh.hpp
+++ b/stan/math/fwd/fun/cosh.hpp
@@ -3,7 +3,10 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/fun/exp.hpp>
+#include <stan/math/prim/fun/cosh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -15,6 +18,18 @@ inline fvar<T> cosh(const fvar<T>& x) {
   return fvar<T>(cosh(x.val_), x.d_ * sinh(x.val_));
 }
 
+/**
+ * Return the hyperbolic cosine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return hyperbolic cosine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> cosh(const std::complex<fvar<T>>& z) {
+  return stan::math::internal::complex_cosh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/exp.hpp b/stan/math/fwd/fun/exp.hpp
index c15b8a2637b..d0690fd803d 100644
--- a/stan/math/fwd/fun/exp.hpp
+++ b/stan/math/fwd/fun/exp.hpp
@@ -3,7 +3,9 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/exp.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -13,6 +15,18 @@ inline fvar<T> exp(const fvar<T>& x) {
   return fvar<T>(exp(x.val_), x.d_ * exp(x.val_));
 }
 
+/**
+ * Return the natural exponentiation (base e) of the specified complex number.
+ *
+ * @tparam T value type of autodiff variable
+ * @param z complex argument
+ * @return exponentiation of argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> exp(const std::complex<fvar<T>>& z) {
+  return internal::complex_exp(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/log.hpp b/stan/math/fwd/fun/log.hpp
index f1d5d1d96c7..d7ce96f0bfa 100644
--- a/stan/math/fwd/fun/log.hpp
+++ b/stan/math/fwd/fun/log.hpp
@@ -4,7 +4,9 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/log.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -18,6 +20,19 @@ inline fvar<T> log(const fvar<T>& x) {
     return fvar<T>(log(x.val_), x.d_ / x.val_);
   }
 }
+
+/**
+ * Return the natural logarithm (base e) of the specified complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param z complex argument
+ * @return natural logarithm of argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> log(const std::complex<fvar<T>>& z) {
+  return internal::complex_log(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/log10.hpp b/stan/math/fwd/fun/log10.hpp
index 189f55e9a29..f418eeb4b37 100644
--- a/stan/math/fwd/fun/log10.hpp
+++ b/stan/math/fwd/fun/log10.hpp
@@ -4,7 +4,9 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/log10.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -20,6 +22,18 @@ inline fvar<T> log10(const fvar<T>& x) {
   }
 }
 
+/**
+ * Return the base 10 logarithm of the specified complex number.
+ *
+ * @tparam T autodiff value type
+ * @param z complex argument
+ * @return base 10 log of argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> log10(const std::complex<fvar<T>>& z) {
+  return internal::complex_log10(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/norm.hpp b/stan/math/fwd/fun/norm.hpp
new file mode 100644
index 00000000000..1de3e640223
--- /dev/null
+++ b/stan/math/fwd/fun/norm.hpp
@@ -0,0 +1,26 @@
+#ifndef STAN_MATH_FWD_FUN_NORM_HPP
+#define STAN_MATH_FWD_FUN_NORM_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/norm.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the squared magnitude of the complex argument.
+ *
+ * @tparam T value type of autodiff variable
+ * @param[in] z argument
+ * @return phase squared magnitude of the argument
+ */
+template <typename T>
+inline fvar<T> norm(const std::complex<fvar<T>>& z) {
+  return internal::complex_norm(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/fwd/fun/polar.hpp b/stan/math/fwd/fun/polar.hpp
new file mode 100644
index 00000000000..891c2f0d500
--- /dev/null
+++ b/stan/math/fwd/fun/polar.hpp
@@ -0,0 +1,56 @@
+#ifndef STAN_MATH_FWD_FUN_POLAR_HPP
+#define STAN_MATH_FWD_FUN_POLAR_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/fun/polar.hpp>
+#include <complex>
+#include <type_traits>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @tparam T autodiff value type
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+template <typename T>
+inline std::complex<fvar<T>> polar(const fvar<T>& r, const fvar<T>& theta) {
+  return internal::complex_polar(r, theta);
+}
+
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @tparam T autodiff value type for magnitude
+ * @tparam U arithmetic type for phase angle
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+template <typename T, typename U>
+inline std::complex<fvar<T>> polar(const fvar<T>& r, U theta) {
+  return internal::complex_polar(r, theta);
+}
+
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @tparam T autodiff value type for phase angle
++* * @tparam U arithmetic type for magnitude
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+template <typename T, typename U>
+inline std::complex<fvar<T>> polar(U r, const fvar<T>& theta) {
+  return internal::complex_polar(r, theta);
+}
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/fwd/fun/pow.hpp b/stan/math/fwd/fun/pow.hpp
index 71ecb438d0e..c181df79b62 100644
--- a/stan/math/fwd/fun/pow.hpp
+++ b/stan/math/fwd/fun/pow.hpp
@@ -7,7 +7,10 @@
 #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/prim/fun/pow.hpp>
 #include <cmath>
+#include <complex>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -21,19 +24,18 @@ inline fvar<T> pow(const fvar<T>& x1, const fvar<T>& x2) {
                                 * pow_x1_x2);
 }
 
-template <typename T>
-inline fvar<T> pow(double x1, const fvar<T>& x2) {
+template <typename T, typename U, typename = require_arithmetic_t<U>>
+inline fvar<T> pow(U x1, const fvar<T>& 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>
-inline fvar<T> pow(const fvar<T>& x1, double x2) {
+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);
   }
@@ -52,9 +54,166 @@ inline fvar<T> pow(const fvar<T>& x1, double x2) {
   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.
+ *
+ * @tparam T autodiff value type
+ * @param x first argument
+ * @param y second argument
+ * @return 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);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/proj.hpp b/stan/math/fwd/fun/proj.hpp
new file mode 100644
index 00000000000..931d2dc4b72
--- /dev/null
+++ b/stan/math/fwd/fun/proj.hpp
@@ -0,0 +1,27 @@
+#ifndef STAN_MATH_FWD_FUN_PROJ_HPP
+#define STAN_MATH_FWD_FUN_PROJ_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/fun/proj.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the projection of the complex argument onto the Riemann
+ * sphere.
+ *
+ * @tparam T value type of autodiff variable
+ * @param[in] z argument
+ * @return projection of the argument onto the Riemann sphere
+ */
+template <typename T>
+inline std::complex<fvar<T>> proj(const std::complex<fvar<T>>& z) {
+  return internal::complex_proj(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/fwd/fun/sin.hpp b/stan/math/fwd/fun/sin.hpp
index 2adc1ac6473..c938cf4948d 100644
--- a/stan/math/fwd/fun/sin.hpp
+++ b/stan/math/fwd/fun/sin.hpp
@@ -1,8 +1,11 @@
 #ifndef STAN_MATH_FWD_FUN_SIN_HPP
 #define STAN_MATH_FWD_FUN_SIN_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/fun/sin.hpp>
+#include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -13,6 +16,19 @@ inline fvar<T> sin(const fvar<T>& x) {
   using std::sin;
   return fvar<T>(sin(x.val_), x.d_ * cos(x.val_));
 }
+
+/**
+ * Return the sine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return sine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> sin(const std::complex<fvar<T>>& z) {
+  return internal::complex_sin(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/sinh.hpp b/stan/math/fwd/fun/sinh.hpp
index 2260cf7f8bf..0a7e6d4343a 100644
--- a/stan/math/fwd/fun/sinh.hpp
+++ b/stan/math/fwd/fun/sinh.hpp
@@ -1,8 +1,10 @@
 #ifndef STAN_MATH_FWD_FUN_SINH_HPP
 #define STAN_MATH_FWD_FUN_SINH_HPP
 
+#include <stan/math/prim/fun/sinh.hpp>
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -13,6 +15,19 @@ inline fvar<T> sinh(const fvar<T>& x) {
   using std::sinh;
   return fvar<T>(sinh(x.val_), x.d_ * cosh(x.val_));
 }
+
+/**
+ * Return the hyperbolic sine of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return hyperbolic sine of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> sinh(const std::complex<fvar<T>>& z) {
+  return internal::complex_sinh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/sqrt.hpp b/stan/math/fwd/fun/sqrt.hpp
index 92b712cf6f6..e72bcea4a2f 100644
--- a/stan/math/fwd/fun/sqrt.hpp
+++ b/stan/math/fwd/fun/sqrt.hpp
@@ -4,7 +4,11 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/prim/fun/inv_sqrt.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
+#include <stan/math/fwd/fun/inv_sqrt.hpp>
+#include <stan/math/fwd/fun/sqrt.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -14,6 +18,19 @@ inline fvar<T> sqrt(const fvar<T>& x) {
   using std::sqrt;
   return fvar<T>(sqrt(x.val_), 0.5 * x.d_ * inv_sqrt(x.val_));
 }
+
+/**
+ * Return the square root of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return square root of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> sqrt(const std::complex<fvar<T>>& z) {
+  return internal::complex_sqrt(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/tan.hpp b/stan/math/fwd/fun/tan.hpp
index 7fe5f92de5e..503ca10e2ef 100644
--- a/stan/math/fwd/fun/tan.hpp
+++ b/stan/math/fwd/fun/tan.hpp
@@ -1,8 +1,11 @@
 #ifndef STAN_MATH_FWD_FUN_TAN_HPP
 #define STAN_MATH_FWD_FUN_TAN_HPP
 
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/fun/tan.hpp>
+#include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -13,6 +16,18 @@ inline fvar<T> tan(const fvar<T>& x) {
   using std::tan;
   return fvar<T>(tan(x.val_), x.d_ / (cos(x.val_) * cos(x.val_)));
 }
+
+/**
+ * Return the tangent of the complex argument.
+ *
+ * @param[in] z argument
+ * @return tangent of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> tan(const std::complex<fvar<T>>& z) {
+  return stan::math::internal::complex_tan(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/fun/tanh.hpp b/stan/math/fwd/fun/tanh.hpp
index 2dff5e3ee58..9bb120ba4e3 100644
--- a/stan/math/fwd/fun/tanh.hpp
+++ b/stan/math/fwd/fun/tanh.hpp
@@ -3,7 +3,10 @@
 
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/fun/exp.hpp>
+#include <stan/math/prim/fun/tanh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -14,6 +17,19 @@ inline fvar<T> tanh(const fvar<T>& x) {
   T u = tanh(x.val_);
   return fvar<T>(u, x.d_ * (1 - u * u));
 }
+
+/**
+ * Return the hyperbolic tangent of the complex argument.
+ *
+ * @tparam T autodiff value type
+ * @param[in] z argument
+ * @return hyperbolic tangent of the argument
+ */
+template <typename T>
+inline std::complex<fvar<T>> tanh(const std::complex<fvar<T>>& z) {
+  return stan::math::internal::complex_tanh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/core.hpp b/stan/math/prim/core.hpp
index 680c7b0139f..3a96e419fd6 100644
--- a/stan/math/prim/core.hpp
+++ b/stan/math/prim/core.hpp
@@ -2,5 +2,13 @@
 #define STAN_MATH_PRIM_CORE_HPP
 
 #include <stan/math/prim/core/init_threadpool_tbb.hpp>
+#include <stan/math/prim/core/operator_addition.hpp>
+#include <stan/math/prim/core/operator_division.hpp>
+#include <stan/math/prim/core/operator_equal_equal.hpp>
+#include <stan/math/prim/core/operator_minus.hpp>
+#include <stan/math/prim/core/operator_multiplication.hpp>
+#include <stan/math/prim/core/operator_not_equal.hpp>
+#include <stan/math/prim/core/operator_plus.hpp>
+#include <stan/math/prim/core/operator_subtraction.hpp>
 
 #endif
diff --git a/stan/math/prim/core/complex_base.hpp b/stan/math/prim/core/complex_base.hpp
index 4d21c88b585..222a93009ef 100644
--- a/stan/math/prim/core/complex_base.hpp
+++ b/stan/math/prim/core/complex_base.hpp
@@ -39,7 +39,7 @@ class complex_base {
    * @tparam U real type (assignable to `value_type`)
    * @param[in] re real part
    */
-  template <typename U, typename = require_stan_scalar_t<U>>
+  template <typename U>  // , typename = require_stan_scalar_t<U>>
   complex_base(const U& re) : re_(re) {}  // NOLINT(runtime/explicit)
 
   /**
diff --git a/stan/math/prim/core/operator_addition.hpp b/stan/math/prim/core/operator_addition.hpp
new file mode 100644
index 00000000000..030e3e7e69b
--- /dev/null
+++ b/stan/math/prim/core/operator_addition.hpp
@@ -0,0 +1,75 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_ADDITION_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_ADDITION_HPP
+
+#include <stan/math/prim/meta/return_type.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+/**
+ * Return the sum of the specified arguments.  At least one of the
+ * arguments must be a complex number.
+ *
+ * @tparam U type of first argument
+ * @tparam V type of second argument
+ * @param[in] lhs first argument
+ * @param[in] rhs second argument
+ * @return sum of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> complex_add(const U& lhs, const V& rhs) {
+  complex_return_t<U, V> y(lhs);
+  y += rhs;
+  return y;
+}
+}  // namespace internal
+
+/**
+ * Return the sum of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return sum of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator+(const std::complex<U>& x,
+                                        const std::complex<V>& y) {
+  return internal::complex_add(x, y);
+}
+
+/**
+ * Return the sum of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return sum of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator+(const std::complex<U>& x, const V& y) {
+  return internal::complex_add(x, y);
+}
+
+/**
+ * Return the sum of the arguments.
+ *
+ * @tparam U type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return sum of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator+(const U& x, const std::complex<V>& y) {
+  return internal::complex_add(x, y);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/core/operator_division.hpp b/stan/math/prim/core/operator_division.hpp
new file mode 100644
index 00000000000..090970b616c
--- /dev/null
+++ b/stan/math/prim/core/operator_division.hpp
@@ -0,0 +1,75 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_DIVISION_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_DIVISION_HPP
+
+#include <stan/math/prim/meta/return_type.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+/**
+ * Return the quotient of the specified arguments.  At least one of the
+ * arguments must be a complex number.
+ *
+ * @tparam U type of first argumentx
+ * @tparam V type of second argument
+ * @param[in] lhs first argument
+ * @param[in] rhs second argument
+ * @return quotient of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> complex_divide(const U& lhs, const V& rhs) {
+  complex_return_t<U, V> y(lhs);
+  y /= rhs;
+  return y;
+}
+}  // namespace internal
+
+/**
+ * Return the quotient of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return quotient of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator/(const std::complex<U>& x,
+                                        const std::complex<V>& y) {
+  return internal::complex_divide(x, y);
+}
+
+/**
+ * Return the quotient of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return quotient of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator/(const std::complex<U>& x, const V& y) {
+  return internal::complex_divide(x, y);
+}
+
+/**
+ * Return the quotient of the arguments.
+ *
+ * @tparam U type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return quotient of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator/(const U& x, const std::complex<V>& y) {
+  return internal::complex_divide(x, y);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/core/operator_equal_equal.hpp b/stan/math/prim/core/operator_equal_equal.hpp
new file mode 100644
index 00000000000..4c66cd90b5a
--- /dev/null
+++ b/stan/math/prim/core/operator_equal_equal.hpp
@@ -0,0 +1,58 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_EQUAL_EQUAL_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_EQUAL_EQUAL_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return `true` if the complex numbers have equal imaginary and
+ * complex parts.
+ *
+ * @tparam U value type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return `true` if the arguments are equal
+ */
+template <typename U, typename V, typename = require_any_autodiff_t<U, V>>
+inline bool operator==(const std::complex<U>& x, const std::complex<V>& y) {
+  return x.real() == y.real() && x.imag() == y.imag();
+}
+
+/**
+ * Return `true` if the first argument's real part is equal to the
+ * second argument and the first argument's imaginary part is zero.
+ *
+ * @tparam U value type of first argument
+ * @tparam V type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return `true` if the arguments are equal
+ */
+template <typename U, typename V, typename = require_any_autodiff_t<U, V>>
+inline bool operator==(const std::complex<U>& x, const V& y) {
+  return x.real() == y && x.imag() == 0;
+}
+
+/**
+ * Return `true` if the first argument is equal to the real part of
+ * the second argument and the imaginary part of the second argument
+ * is zero.
+ *
+ * @tparam U type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return `true` if the arguments are equal
+ */
+template <typename U, typename V, typename = require_any_autodiff_t<U, V>>
+inline bool operator==(const U& x, const std::complex<V>& y) {
+  return x == y.real() && 0 == y.imag();
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/core/operator_minus.hpp b/stan/math/prim/core/operator_minus.hpp
new file mode 100644
index 00000000000..e1b238aff5c
--- /dev/null
+++ b/stan/math/prim/core/operator_minus.hpp
@@ -0,0 +1,25 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_MINUS_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_MINUS_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the negation of the argument.
+ *
+ * @tparam U value type argument
+ * @param x argument
+ * @return negation of the argument
+ */
+template <typename U, require_autodiff_t<U>>
+inline std::complex<U> operator-(const std::complex<U>& x) {
+  return {-x.real(), -x.imag()};
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/core/operator_multiplication.hpp b/stan/math/prim/core/operator_multiplication.hpp
new file mode 100644
index 00000000000..bb70715033b
--- /dev/null
+++ b/stan/math/prim/core/operator_multiplication.hpp
@@ -0,0 +1,75 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_MULTIPLICATION_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_MULTIPLICATION_HPP
+
+#include <stan/math/prim/meta/return_type.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+/**
+ * Return the product of the specified arguments.  At least one of the
+ * arguments must be a complex number.
+ *
+ * @tparam U type of first argumentx
+ * @tparam V type of second argument
+ * @param[in] lhs first argument
+ * @param[in] rhs second argument
+ * @return product of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> complex_multiply(const U& lhs, const V& rhs) {
+  complex_return_t<U, V> y(lhs);
+  y *= rhs;
+  return y;
+}
+}  // namespace internal
+
+/**
+ * Return the product of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return product of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator*(const std::complex<U>& x,
+                                        const std::complex<V>& y) {
+  return internal::complex_multiply(x, y);
+}
+
+/**
+ * Return the product of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return product of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator*(const std::complex<U>& x, const V& y) {
+  return internal::complex_multiply(x, y);
+}
+
+/**
+ * Return the product of the arguments.
+ *
+ * @tparam U type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return product of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator*(const U& x, const std::complex<V>& y) {
+  return internal::complex_multiply(x, y);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/core/operator_not_equal.hpp b/stan/math/prim/core/operator_not_equal.hpp
new file mode 100644
index 00000000000..a3dc733ed51
--- /dev/null
+++ b/stan/math/prim/core/operator_not_equal.hpp
@@ -0,0 +1,60 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_NOT_EQUAL_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_NOT_EQUAL_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return `true` if the complex numbers have unequal imaginary or
+ * complex parts.
+ *
+ * @tparam U value type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return `true` if the arguments are equal
+ */
+template <typename U, typename V, typename = require_any_autodiff_t<U, V>>
+inline bool operator!=(const std::complex<U>& x, const std::complex<V>& y) {
+  return !(x.real() == y.real() && x.imag() == y.imag());
+}
+
+/**
+ * Return `true` if the first argument's real part is unequal to the
+ * second argument or the first argument's imaginary part is unequal
+ * to zero.
+ *
+ * @tparam U value type of first argument
+ * @tparam V type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return `true` if the arguments are equal
+ */
+template <typename U, typename V, typename = require_any_autodiff_t<U, V>>
+inline bool operator!=(const std::complex<U>& x, const V& y) {
+  return !(x.real() == y && x.imag() == 0);
+}
+
+/**
+ * Return `true` if the first argument is unequal to the real part of
+ * the second argument or the imaginary part of the second argument
+ * is nonzero.
+ *
+ * @tparam U type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return `true` if the arguments are not equal
+ */
+template <typename U, typename V, typename = require_any_autodiff_t<U, V>>
+inline bool operator!=(const U& x, const std::complex<V>& y) {
+  return !(x == y.real() && 0 == y.imag());
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/core/operator_plus.hpp b/stan/math/prim/core/operator_plus.hpp
new file mode 100644
index 00000000000..cee306e56b5
--- /dev/null
+++ b/stan/math/prim/core/operator_plus.hpp
@@ -0,0 +1,25 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_PLUS_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_PLUS_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the argument.
+ *
+ * @tparam U value type argument
+ * @param x argument
+ * @return argument
+ */
+template <typename U, require_autodiff_t<U>>
+inline std::complex<U> operator+(const std::complex<U>& x) {
+  return x;
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/core/operator_subtraction.hpp b/stan/math/prim/core/operator_subtraction.hpp
new file mode 100644
index 00000000000..cc1c35bc2b3
--- /dev/null
+++ b/stan/math/prim/core/operator_subtraction.hpp
@@ -0,0 +1,75 @@
+#ifndef STAN_MATH_PRIM_CORE_OPERATOR_SUBTRACTION_HPP
+#define STAN_MATH_PRIM_CORE_OPERATOR_SUBTRACTION_HPP
+
+#include <stan/math/prim/meta/return_type.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+namespace internal {
+/**
+ * Return the difference between specified arguments.  At least one of
+ * the arguments must be a complex number.
+ *
+ * @tparam U type of first argument
+ * @tparam V type of second argument
+ * @param[in] lhs first argument
+ * @param[in] rhs second argument
+ * @return sum of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> complex_subtract(const U& lhs, const V& rhs) {
+  complex_return_t<U, V> y(lhs);
+  y -= rhs;
+  return y;
+}
+}  // namespace internal
+
+/**
+ * Return the difference of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return difference of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator-(const std::complex<U>& x,
+                                        const std::complex<V>& y) {
+  return internal::complex_subtract(x, y);
+}
+
+/**
+ * Return the difference of the arguments.
+ *
+ * @tparam U value type of first argument
+ * @tparam V type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return difference of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator-(const std::complex<U>& x, const V& y) {
+  return internal::complex_subtract(x, y);
+}
+
+/**
+ * Return the difference of the arguments.
+ *
+ * @tparam U type of first argument
+ * @tparam V value type of second argument
+ * @param x first argument
+ * @param y second argument
+ * @return difference of the arguments
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> operator-(const U& x, const std::complex<V>& y) {
+  return internal::complex_subtract(x, y);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp
index 72d6a26d557..f214acd63ea 100644
--- a/stan/math/prim/fun.hpp
+++ b/stan/math/prim/fun.hpp
@@ -10,6 +10,7 @@
 #include <stan/math/prim/fun/append_array.hpp>
 #include <stan/math/prim/fun/append_col.hpp>
 #include <stan/math/prim/fun/append_row.hpp>
+#include <stan/math/prim/fun/arg.hpp>
 #include <stan/math/prim/fun/array_builder.hpp>
 #include <stan/math/prim/fun/as_bool.hpp>
 #include <stan/math/prim/fun/asin.hpp>
@@ -39,6 +40,7 @@
 #include <stan/math/prim/fun/cols.hpp>
 #include <stan/math/prim/fun/columns_dot_product.hpp>
 #include <stan/math/prim/fun/columns_dot_self.hpp>
+#include <stan/math/prim/fun/conj.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/copysign.hpp>
 #include <stan/math/prim/fun/corr_constrain.hpp>
@@ -115,7 +117,9 @@
 #include <stan/math/prim/fun/identity_constrain.hpp>
 #include <stan/math/prim/fun/identity_free.hpp>
 #include <stan/math/prim/fun/identity_matrix.hpp>
+#include <stan/math/prim/fun/i_times.hpp>
 #include <stan/math/prim/fun/if_else.hpp>
+#include <stan/math/prim/fun/imag.hpp>
 #include <stan/math/prim/fun/inc_beta.hpp>
 #include <stan/math/prim/fun/initialize.hpp>
 #include <stan/math/prim/fun/int_step.hpp>
@@ -147,6 +151,7 @@
 #include <stan/math/prim/fun/linspaced_vector.hpp>
 #include <stan/math/prim/fun/lmgamma.hpp>
 #include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/logb.hpp>
 #include <stan/math/prim/fun/log10.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/log1m_exp.hpp>
@@ -204,6 +209,7 @@
 #include <stan/math/prim/fun/multiply.hpp>
 #include <stan/math/prim/fun/multiply_log.hpp>
 #include <stan/math/prim/fun/multiply_lower_tri_self_transpose.hpp>
+#include <stan/math/prim/fun/norm.hpp>
 #include <stan/math/prim/fun/num_elements.hpp>
 #include <stan/math/prim/fun/offset_multiplier_constrain.hpp>
 #include <stan/math/prim/fun/offset_multiplier_free.hpp>
@@ -220,14 +226,17 @@
 #include <stan/math/prim/fun/owens_t.hpp>
 #include <stan/math/prim/fun/Phi.hpp>
 #include <stan/math/prim/fun/Phi_approx.hpp>
+#include <stan/math/prim/fun/polar.hpp>
 #include <stan/math/prim/fun/positive_constrain.hpp>
 #include <stan/math/prim/fun/positive_free.hpp>
 #include <stan/math/prim/fun/positive_ordered_constrain.hpp>
 #include <stan/math/prim/fun/positive_ordered_free.hpp>
+#include <stan/math/prim/fun/pow.hpp>
 #include <stan/math/prim/fun/primitive_value.hpp>
 #include <stan/math/prim/fun/prob_constrain.hpp>
 #include <stan/math/prim/fun/prob_free.hpp>
 #include <stan/math/prim/fun/prod.hpp>
+#include <stan/math/prim/fun/proj.hpp>
 #include <stan/math/prim/fun/promote_elements.hpp>
 #include <stan/math/prim/fun/promote_scalar.hpp>
 #include <stan/math/prim/fun/qr_Q.hpp>
@@ -242,6 +251,7 @@
 #include <stan/math/prim/fun/read_corr_matrix.hpp>
 #include <stan/math/prim/fun/read_cov_L.hpp>
 #include <stan/math/prim/fun/read_cov_matrix.hpp>
+#include <stan/math/prim/fun/real.hpp>
 #include <stan/math/prim/fun/rep_array.hpp>
 #include <stan/math/prim/fun/rep_matrix.hpp>
 #include <stan/math/prim/fun/rep_row_vector.hpp>
@@ -254,6 +264,7 @@
 #include <stan/math/prim/fun/rows.hpp>
 #include <stan/math/prim/fun/rows_dot_product.hpp>
 #include <stan/math/prim/fun/rows_dot_self.hpp>
+#include <stan/math/prim/fun/scalbn.hpp>
 #include <stan/math/prim/fun/scale_matrix_exp_multiply.hpp>
 #include <stan/math/prim/fun/scaled_add.hpp>
 #include <stan/math/prim/fun/sd.hpp>
diff --git a/stan/math/prim/fun/abs.hpp b/stan/math/prim/fun/abs.hpp
index cdf0a52dcd6..4a4d76350c4 100644
--- a/stan/math/prim/fun/abs.hpp
+++ b/stan/math/prim/fun/abs.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_PRIM_FUN_ABS_HPP
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/hypot.hpp>
 #include <cmath>
 
 namespace stan {
@@ -18,6 +19,20 @@ namespace math {
  */
 inline double abs(double x) { return std::fabs(x); }
 
+namespace internal {
+/**
+ * Return the absolute value of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return absolute value of the argument
+ */
+template <typename V>
+inline V complex_abs(const std::complex<V>& z) {
+  return hypot(z.real(), z.imag());
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/acos.hpp b/stan/math/prim/fun/acos.hpp
index 65b80388ad3..88317641b0e 100644
--- a/stan/math/prim/fun/acos.hpp
+++ b/stan/math/prim/fun/acos.hpp
@@ -1,9 +1,17 @@
 #ifndef STAN_MATH_PRIM_FUN_ACOS_HPP
 #define STAN_MATH_PRIM_FUN_ACOS_HPP
 
+#include <stan/math/prim/core.hpp>
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/asin.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
+#include <stan/math/prim/fun/isnan.hpp>
+#include <stan/math/prim/fun/polar.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -53,6 +61,20 @@ inline auto acos(const Container& x) {
       x, [](const auto& v) { return v.array().acos(); });
 }
 
+namespace internal {
+/**
+ * Return the arc cosine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return arc cosine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_acos(const std::complex<V>& z) {
+  return 0.5 * pi() - asin(z);
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/acosh.hpp b/stan/math/prim/fun/acosh.hpp
index e0374ce8031..015bc014112 100644
--- a/stan/math/prim/fun/acosh.hpp
+++ b/stan/math/prim/fun/acosh.hpp
@@ -3,9 +3,13 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/copysign.hpp>
 #include <stan/math/prim/fun/is_inf.hpp>
 #include <stan/math/prim/fun/is_nan.hpp>
+#include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -83,6 +87,24 @@ inline auto acosh(const T& x) {
   return apply_scalar_unary<acosh_fun, T>::apply(x);
 }
 
+namespace internal {
+
+/**
+ * Return the hyperbolic arc cosine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return hyperbolic arc cosine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_acosh(const std::complex<V>& z) {
+  std::complex<double> y_d = acosh(value_of_rec(z));
+  auto y = log(z + sqrt(z * z - 1));
+  return copysign(y, y_d);
+}
+
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/arg.hpp b/stan/math/prim/fun/arg.hpp
new file mode 100644
index 00000000000..e8d5029c904
--- /dev/null
+++ b/stan/math/prim/fun/arg.hpp
@@ -0,0 +1,26 @@
+#ifndef STAN_MATH_PRIM_FUN_ARG_HPP
+#define STAN_MATH_PRIM_FUN_ARG_HPP
+
+#include <cmath>
+#include <complex>
+
+namespace stan {
+namespace math {
+namespace internal {
+/**
+ * Return the phase angle of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return phase angle of the argument
+ */
+template <typename V>
+inline V complex_arg(const std::complex<V>& z) {
+  using std::atan2;
+  return atan2(z.imag(), z.real());
+}
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/asin.hpp b/stan/math/prim/fun/asin.hpp
index 8ddd0272d05..c5db2b43a15 100644
--- a/stan/math/prim/fun/asin.hpp
+++ b/stan/math/prim/fun/asin.hpp
@@ -1,9 +1,15 @@
 #ifndef STAN_MATH_PRIM_FUN_ASIN_HPP
 #define STAN_MATH_PRIM_FUN_ASIN_HPP
 
+#include <stan/math/prim/core.hpp>
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/asinh.hpp>
+#include <stan/math/prim/fun/copysign.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/i_times.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -53,6 +59,22 @@ inline auto asin(const Container& x) {
       x, [](const auto& v) { return v.array().asin(); });
 }
 
+namespace internal {
+/**
+ * Return the arc sine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return arc sine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_asin(const std::complex<V>& z) {
+  auto y_d = asin(value_of_rec(z));
+  auto y = neg_i_times(asinh(i_times(z)));
+  return copysign(y, y_d);
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/asinh.hpp b/stan/math/prim/fun/asinh.hpp
index 7f3fe63476a..06eb3726b8b 100644
--- a/stan/math/prim/fun/asinh.hpp
+++ b/stan/math/prim/fun/asinh.hpp
@@ -1,8 +1,20 @@
 #ifndef STAN_MATH_PRIM_FUN_ASINH_HPP
 #define STAN_MATH_PRIM_FUN_ASINH_HPP
 
+#include <stan/math/prim/core.hpp>
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/abs.hpp>
+#include <stan/math/prim/fun/arg.hpp>
+#include <stan/math/prim/fun/copysign.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/isnan.hpp>
+#include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/polar.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -35,6 +47,22 @@ inline auto asinh(const T& x) {
   return apply_scalar_unary<asinh_fun, T>::apply(x);
 }
 
+namespace internal {
+/**
+ * Return the hyperbolic arc sine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return hyperbolic arc sine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_asinh(const std::complex<V>& z) {
+  std::complex<double> y_d = asinh(value_of_rec(z));
+  auto y = log(z + sqrt(1 + z * z));
+  return copysign(y, y_d);
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/atan.hpp b/stan/math/prim/fun/atan.hpp
index 90e95a04c5d..d40ada7ff6f 100644
--- a/stan/math/prim/fun/atan.hpp
+++ b/stan/math/prim/fun/atan.hpp
@@ -1,9 +1,13 @@
 #ifndef STAN_MATH_PRIM_FUN_ATAN_HPP
 #define STAN_MATH_PRIM_FUN_ATAN_HPP
 
+#include <stan/math/prim/core.hpp>
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/atanh.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/i_times.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -53,6 +57,20 @@ inline auto atan(const Container& x) {
       x, [](const auto& v) { return v.array().atan(); });
 }
 
+namespace internal {
+/**
+ * Return the arc tangent of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return arc tangent of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_atan(const std::complex<V>& z) {
+  return neg_i_times(atanh(i_times(z)));
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/atanh.hpp b/stan/math/prim/fun/atanh.hpp
index c04730015ea..d6aeccd5cbb 100644
--- a/stan/math/prim/fun/atanh.hpp
+++ b/stan/math/prim/fun/atanh.hpp
@@ -1,10 +1,14 @@
 #ifndef STAN_MATH_PRIM_FUN_ATANH_HPP
 #define STAN_MATH_PRIM_FUN_ATANH_HPP
 
-#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/core.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/copysign.hpp>
 #include <stan/math/prim/fun/is_nan.hpp>
+#include <stan/math/prim/fun/log.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -72,6 +76,23 @@ inline auto atanh(const T& x) {
   return apply_scalar_unary<atanh_fun, T>::apply(x);
 }
 
+namespace internal {
+/**
+ * Return the hyperbolic arc tangent of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return hyperbolic arc tangent of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_atanh(const std::complex<V>& z) {
+  std::complex<double> y_d = atanh(value_of_rec(z));
+  V one(1);
+  auto y = 0.5 * (log(one + z) - log(one - z));
+  return copysign(y, y_d);
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/conj.hpp b/stan/math/prim/fun/conj.hpp
new file mode 100644
index 00000000000..a356cd0a9ee
--- /dev/null
+++ b/stan/math/prim/fun/conj.hpp
@@ -0,0 +1,24 @@
+#ifndef STAN_MATH_PRIM_FUN_CONJ_HPP
+#define STAN_MATH_PRIM_FUN_CONJ_HPP
+
+#include <complex>
+
+namespace stan {
+namespace math {
+namespace internal {
+/**
+ * Return the complex conjugate the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return complex conjugate of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_conj(const std::complex<V>& z) {
+  return {z.real(), -z.imag()};
+}
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/copysign.hpp b/stan/math/prim/fun/copysign.hpp
index bccf23bd1cc..3273ab2db1f 100644
--- a/stan/math/prim/fun/copysign.hpp
+++ b/stan/math/prim/fun/copysign.hpp
@@ -2,6 +2,9 @@
 #define STAN_MATH_PRIM_SCAL_FUN_COPYSIGN_HPP
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -16,18 +19,62 @@ namespace math {
  * Overload of `std::copysign` from `cmath` for argument-dependent
  * lookup.
  *
- * @tparam ADType type of first argument
+ * @tparam T type of first argument
  * @tparam U type of second argument
  * @param[in] x first complex argument
  * @param[in] y second complex argument
- * @return copy of second argument, negated if necessary to match sign
- * of first argument
+ * @return copy of the first argument, negated if necessary to match
+ * the sign of the second argument
  */
-template <typename ADType, typename U>
-inline ADType copysign(const ADType& x, const U& y) {
-  // +0 is positive; second condition handles -0
-  return (x < 0 && y >= 0) || (x >= 0 && y < 0) ? -x : x;
+template <typename T, typename U>
+inline T copysign(const T& x, const U& y) {
+  return (std::signbit(value_of_rec(x)) != std::signbit(value_of_rec(y))) ? -x
+                                                                          : x;
 }
+
+/**
+ * Return the negation of the first argument if the first and second
+ * arguments have different signs and the first argument is not zero,
+ * otherwise return a copy of the first argument.
+ *
+ * @tparam T type of first argument
+ * @tparam U type of second argument
+ * @param[in] x first complex argument
+ * @param[in] y second complex argument
+ * @return copy of the first argument, negated if necessary to match
+ * the sign of the second argument
+ * @see copysign
+ */
+template <typename T, typename U>
+inline T copysign_non_zero(const T& x, const U& y) {
+  return x != 0 ? stan::math::copysign(x, y) : x;
+}
+
+/**
+ * Return the complex number composed of the real and complex parts
+ * with signs copied from the real and complex parts of the first
+ * arguments to the real and complex parts of the second.
+ *
+ * This is an overload of the standard libary `copysign` for complex
+ * numbers that will be used with argument-dependent lookup.  Rather
+ * than using the standard library `copysign`, it uses
+ * `copysign_non_zero`, which does not change sign if the reference
+ * value is zero (`-0.0` or `0.0`).
+ *
+ * @tparam T value type of first argument
+ * @tparam U value type of second argument
+ * @param[in] x first complex argument
+ * @param[in] y second complex argument
+ * @return copy of second argument, with components negated if
+ * necessary to match sign of first argument
+ */
+template <typename T, typename U>
+inline std::complex<T> copysign(const std::complex<T>& x,
+                                const std::complex<U>& y) {
+  return {copysign_non_zero(x.real(), y.real()),
+          copysign_non_zero(x.imag(), y.imag())};
+}
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/cos.hpp b/stan/math/prim/fun/cos.hpp
index 48672fd12cb..44fc7636c46 100644
--- a/stan/math/prim/fun/cos.hpp
+++ b/stan/math/prim/fun/cos.hpp
@@ -2,8 +2,11 @@
 #define STAN_MATH_PRIM_FUN_COS_HPP
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/cosh.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/i_times.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -53,6 +56,20 @@ inline auto cos(const Container& x) {
       x, [](const auto& v) { return v.array().cos(); });
 }
 
+namespace internal {
+/**
+ * Return the cosine of the complex argument.
+ *
+ * @tparam T value type of argument
+ * @param[in] z argument
+ * @return cosine of the argument
+ */
+template <typename T>
+inline std::complex<T> complex_cos(const std::complex<T>& z) {
+  return cosh(i_times(z));
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/cosh.hpp b/stan/math/prim/fun/cosh.hpp
index 219a5a0847d..c482d19987a 100644
--- a/stan/math/prim/fun/cosh.hpp
+++ b/stan/math/prim/fun/cosh.hpp
@@ -1,8 +1,10 @@
 #ifndef STAN_MATH_PRIM_FUN_COSH_HPP
 #define STAN_MATH_PRIM_FUN_COSH_HPP
 
+#include <stan/math/prim/core.hpp>
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/exp.hpp>
 #include <cmath>
 
 namespace stan {
@@ -52,6 +54,21 @@ inline auto cosh(const Container& x) {
   return apply_vector_unary<Container>::apply(
       x, [](const auto& v) { return v.array().cosh(); });
 }
+
+namespace internal {
+/**
+ * Return the hyperbolic cosine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return hyperbolic cosine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_cosh(const std::complex<V>& z) {
+  return 0.5 * (exp(z) + exp(-z));
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/exp.hpp b/stan/math/prim/fun/exp.hpp
index e3c3ea7ace7..283c13ed0a1 100644
--- a/stan/math/prim/fun/exp.hpp
+++ b/stan/math/prim/fun/exp.hpp
@@ -4,6 +4,8 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <cmath>
+#include <complex>
+#include <limits>
 
 namespace stan {
 namespace math {
@@ -57,6 +59,45 @@ inline auto exp(const Container& x) {
       x, [](const auto& v) { return v.array().exp(); });
 }
 
+namespace internal {
+/**
+ * Return the natural (base e) complex exponentiation of the specified
+ * complex argument.
+ *
+ * @tparam V value type (must be Stan autodiff type)
+ * @param z complex number
+ * @return natural exponentiation of specified complex number
+ * @see documentation for `std::complex` for boundary condition and
+ * branch cut details
+ */
+template <typename V>
+inline std::complex<V> complex_exp(const std::complex<V>& z) {
+  if (is_inf(z.real()) && z.real() > 0) {
+    if (is_nan(z.imag()) || z.imag() == 0) {
+      // (+inf, nan), (+inf, 0)
+      return z;
+    } else if (is_inf(z.imag()) && z.imag() > 0) {
+      // (+inf, +inf)
+      return {z.real(), std::numeric_limits<double>::quiet_NaN()};
+    } else if (is_inf(z.imag()) && z.imag() < 0) {
+      // (+inf, -inf)
+      return {std::numeric_limits<double>::quiet_NaN(),
+              std::numeric_limits<double>::quiet_NaN()};
+    }
+  }
+  if (is_inf(z.real()) && z.real() < 0
+      && (is_nan(z.imag()) || is_inf(z.imag()))) {
+    // (-inf, nan), (-inf, -inf), (-inf, inf)
+    return {0, 0};
+  }
+  if (is_nan(z.real()) && z.imag() == -0.0) {
+    // (nan, -0)
+    return z;
+  }
+  V exp_re = exp(z.real());
+  return {exp_re * cos(z.imag()), exp_re * sin(z.imag())};
+}
+}  // namespace internal
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/i_times.hpp b/stan/math/prim/fun/i_times.hpp
new file mode 100644
index 00000000000..002f8420edf
--- /dev/null
+++ b/stan/math/prim/fun/i_times.hpp
@@ -0,0 +1,55 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_I_TIMES_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_I_TIMES_HPP
+
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the specified complex number multiplied by `i`.
+ *
+ * This compound function is more efficient than mulitplying by a
+ * constant `i` because it involves only a single arithmetic negation.
+ *
+ * @tparam value type of complex argument
+ * @param[in] z complex argument
+ * @return argument multipled by `i`
+ */
+template <typename T>
+inline std::complex<T> i_times(const std::complex<T>& z) {
+  return {-z.imag(), z.real()};
+}
+
+/**
+ * Return the specified complex number multiplied by `-i`.
+ *
+ * This compound function is more efficient than mulitplying by the
+ * constant `-i` because it involves only a single arithmetic
+ * negation.
+ *
+ * @tparam value type of complex argument
+ * @param[in] z complex argument
+ * @return argument multipled by `-i`
+ */
+template <typename T>
+inline std::complex<T> neg_i_times(const std::complex<T>& z) {
+  return {z.imag(), -z.real()};
+}
+
+/**
+ * Return the complex negation of the specified complex argument.
+ *
+ * @tparam V value type of complex argument
+ * @param[in] z argument
+ * @return negation of argument
+ */
+template <typename V>
+inline std::complex<V> complex_negate(const std::complex<V>& z) {
+  return {-z.real(), -z.imag()};
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/imag.hpp b/stan/math/prim/fun/imag.hpp
new file mode 100644
index 00000000000..f6aed7b3778
--- /dev/null
+++ b/stan/math/prim/fun/imag.hpp
@@ -0,0 +1,25 @@
+#ifndef STAN_MATH_PRIM_FUN_IMAG_HPP
+#define STAN_MATH_PRIM_FUN_IMAG_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the imaginary part of the complex argument.
+ *
+ * @tparam T value type of argument
+ * @param[in] z argument
+ * @return imaginary part of argument
+ */
+template <typename T, require_autodiff_t<T>>
+T imag(const std::complex<T>& z) {
+  return z.imag();
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/log.hpp b/stan/math/prim/fun/log.hpp
index 0dfc1344013..ff50a31c96a 100644
--- a/stan/math/prim/fun/log.hpp
+++ b/stan/math/prim/fun/log.hpp
@@ -3,7 +3,11 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/is_inf.hpp>
+#include <stan/math/prim/fun/is_nan.hpp>
 #include <cmath>
+#include <complex>
+#include <limits>
 
 namespace stan {
 namespace math {
@@ -57,6 +61,28 @@ inline auto log(const Container& x) {
       x, [](const auto& v) { return v.array().log(); });
 }
 
+namespace internal {
+/**
+ * Return the natural logarithm of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return natural logarithm of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_log(const std::complex<V>& z) {
+  static const double inf = std::numeric_limits<double>::infinity();
+  static const double nan = std::numeric_limits<double>::quiet_NaN();
+  if ((is_nan(z.real()) && is_inf(z.imag()))
+      || (is_inf(z.real()) && is_nan(z.imag()))) {
+    return {inf, nan};
+  }
+  V r = sqrt(norm(z));
+  V theta = arg(z);
+  return {log(r), theta};
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/fun/log10.hpp b/stan/math/prim/fun/log10.hpp
index 0476af967ce..7353e83bccd 100644
--- a/stan/math/prim/fun/log10.hpp
+++ b/stan/math/prim/fun/log10.hpp
@@ -3,7 +3,9 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/log.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -52,6 +54,21 @@ inline auto log10(const Container& x) {
       x, [](const auto& v) { return v.array().log10(); });
 }
 
+namespace internal {
+/**
+ * Return the base 10 logarithm of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return base 10 logarithm of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_log10(const std::complex<V>& z) {
+  static const double inv_log_10 = 1 / std::log(10);
+  return log(z) * inv_log_10;
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/logb.hpp b/stan/math/prim/fun/logb.hpp
new file mode 100644
index 00000000000..f0df32a4076
--- /dev/null
+++ b/stan/math/prim/fun/logb.hpp
@@ -0,0 +1,31 @@
+#ifndef STAN_MATH_PRIM_FUN_LOGB_HPP
+#define STAN_MATH_PRIM_FUN_LOGB_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the value of the unbiased radix-independent exponent from
+ * the floating-point argument.
+ *
+ * Implementation note: This implementation works for all autodiff
+ * types because it is a step function, so can return a `double` which
+ * will produce the correct zero derivative if promoted to an autodiff
+ * variable.
+ *
+ * @tparam T floating-point type
+ * @param[in] x floating-point argument
+ * @return unbiased radix-independent exponent of the argument
+ */
+template <typename T, typename = require_autodiff_t<T>>
+double logb(const T& x) {
+  return std::logb(value_of_rec(x));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/norm.hpp b/stan/math/prim/fun/norm.hpp
new file mode 100644
index 00000000000..89fe76552f6
--- /dev/null
+++ b/stan/math/prim/fun/norm.hpp
@@ -0,0 +1,25 @@
+#ifndef STAN_MATH_PRIM_FUN_NORM_HPP
+#define STAN_MATH_PRIM_FUN_NORM_HPP
+
+#include <cmath>
+#include <complex>
+
+namespace stan {
+namespace math {
+namespace internal {
+/**
+ * Return the squared magnitude of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return squared magnitude of the argument
+ */
+template <typename V>
+inline V complex_norm(const std::complex<V>& z) {
+  return square(z.real()) + square(z.imag());
+}
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/polar.hpp b/stan/math/prim/fun/polar.hpp
new file mode 100644
index 00000000000..9ba15ee3aae
--- /dev/null
+++ b/stan/math/prim/fun/polar.hpp
@@ -0,0 +1,56 @@
+#ifndef STAN_MATH_PRIM_FUN_POLAR_HPP
+#define STAN_MATH_PRIM_FUN_POLAR_HPP
+
+#include <stan/math/prim/fun/cos.hpp>
+#include <stan/math/prim/fun/is_inf.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <cmath>
+#include <complex>
+#include <limits>
+
+namespace stan {
+namespace math {
+namespace internal {
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> complex_polar(const U& r, const V& theta) {
+  using std::cos;
+  using std::sin;
+  if (!(r >= 0) || is_inf(theta)) {
+    return {std::numeric_limits<double>::quiet_NaN()};
+  }
+  return {r * cos(theta), r * sin(theta)};
+}
+}  // namespace internal
+
+/**
+ * Returns the complex number with specified magnitude and phase angle.
+ *
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+inline std::complex<double> polar(double r, double theta) {
+  return internal::complex_polar(r, theta);
+}
+inline std::complex<double> polar(double r, int theta) {
+  return internal::complex_polar(r, static_cast<double>(theta));
+}
+inline std::complex<double> polar(int r, double theta) {
+  return internal::complex_polar(static_cast<double>(r), theta);
+}
+inline std::complex<double> polar(int r, int theta) {
+  return internal::complex_polar(static_cast<double>(r),
+                                 static_cast<double>(theta));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/pow.hpp b/stan/math/prim/fun/pow.hpp
new file mode 100644
index 00000000000..da61000e00f
--- /dev/null
+++ b/stan/math/prim/fun/pow.hpp
@@ -0,0 +1,49 @@
+#ifndef STAN_MATH_PRIM_FUN_POW_HPP
+#define STAN_MATH_PRIM_FUN_POW_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <cmath>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+// This overload is required becuase MINGW32 does not include the
+// specialization std::pow(double, int) as required by the C++11
+// standard library interface specification.  Use with other compilers
+// than MINGW32 would introduce an ambiguity with the standard library.
+#if __MINGW32__
+/**
+ * Return the first argument raised to the power of the second
+ * argument.
+ *
+ * @param x base
+ * @param y exponent
+ * @return base raised to the power of the exponent
+ */
+inline double pow(double x, int y) {
+  return std::pow(x, static_cast<double>(y));
+}
+#endif
+
+namespace internal {
+
+/**
+ * Return the first argument raised to the power of the second
+ * argument.  At least one of the arguments must be a complex number.
+ *
+ * @tparam U type of base
+ * @tparam V type of exponent
+ * @param[in] x base
+ * @param[in] y exponent
+ * @return base raised to the power of the exponent
+ */
+template <typename U, typename V>
+inline complex_return_t<U, V> complex_pow(const U& x, const V& y) {
+  return exp(y * log(x));
+}
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/fun/proj.hpp b/stan/math/prim/fun/proj.hpp
new file mode 100644
index 00000000000..a8cb7439a95
--- /dev/null
+++ b/stan/math/prim/fun/proj.hpp
@@ -0,0 +1,30 @@
+#ifndef STAN_MATH_PRIM_FUN_PROJ_HPP
+#define STAN_MATH_PRIM_FUN_PROJ_HPP
+
+#include <stan/math/prim/fun/is_inf.hpp>
+#include <complex>
+#include <limits>
+
+namespace stan {
+namespace math {
+namespace internal {
+/**
+ * Return the projection of the complex argument onto the Riemann
+ * sphere.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return projection of the argument onto the Riemann sphere
+ */
+template <typename V>
+inline std::complex<V> complex_proj(const std::complex<V>& z) {
+  if (is_inf(z.real()) || is_inf(z.imag())) {
+    return {std::numeric_limits<V>::infinity(), z.imag() < 0 ? -0.0 : 0.0};
+  }
+  return z;
+}
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/real.hpp b/stan/math/prim/fun/real.hpp
new file mode 100644
index 00000000000..a14c0d404d0
--- /dev/null
+++ b/stan/math/prim/fun/real.hpp
@@ -0,0 +1,25 @@
+#ifndef STAN_MATH_PRIM_FUN_REAL_HPP
+#define STAN_MATH_PRIM_FUN_REAL_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the real part of the complex argument.
+ *
+ * @tparam T value type of argument
+ * @param[in] z argument
+ * @return real part of argument
+ */
+template <typename T, require_autodiff_t<T>>
+T real(const std::complex<T>& z) {
+  return z.real();
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/scalbn.hpp b/stan/math/prim/fun/scalbn.hpp
new file mode 100644
index 00000000000..6c3e54cdc24
--- /dev/null
+++ b/stan/math/prim/fun/scalbn.hpp
@@ -0,0 +1,19 @@
+#ifndef STAN_MATH_PRIM_FUN_SCALBN_HPP
+#define STAN_MATH_PRIM_FUN_SCALBN_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+template <typename T, typename = require_autodiff_t<T>>
+double scalbn(const T& x, int n) {
+  return std::scalbn(value_of_rec(x), n);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/fun/sin.hpp b/stan/math/prim/fun/sin.hpp
index 9bcd0da5052..a12f85aade4 100644
--- a/stan/math/prim/fun/sin.hpp
+++ b/stan/math/prim/fun/sin.hpp
@@ -3,7 +3,10 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/i_times.hpp>
+#include <stan/math/prim/fun/sinh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -51,6 +54,20 @@ inline auto sin(const Container& x) {
       x, [&](const auto& v) { return v.array().sin(); });
 }
 
+namespace internal {
+/**
+ * Return the sine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return sine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_sin(const std::complex<V>& z) {
+  return neg_i_times(sinh(i_times(z)));
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/sinh.hpp b/stan/math/prim/fun/sinh.hpp
index 7a9486ed838..272c1970016 100644
--- a/stan/math/prim/fun/sinh.hpp
+++ b/stan/math/prim/fun/sinh.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -52,6 +53,20 @@ inline auto sinh(const Container& x) {
       x, [](const auto& v) { return v.array().sinh(); });
 }
 
+namespace internal {
+/**
+ * Return the hyperbolic sine of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return hyperbolic sine of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_sinh(const std::complex<V>& z) {
+  return 0.5 * (exp(z) - exp(-z));
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/sqrt.hpp b/stan/math/prim/fun/sqrt.hpp
index 2ff7cac8d9f..91b32d04932 100644
--- a/stan/math/prim/fun/sqrt.hpp
+++ b/stan/math/prim/fun/sqrt.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -52,6 +53,22 @@ inline auto sqrt(const Container& x) {
       x, [](const auto& v) { return v.array().sqrt(); });
 }
 
+namespace internal {
+/**
+ * Return the square root of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return square root of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_sqrt(const std::complex<V>& z) {
+  auto m = sqrt(hypot(z.real(), z.imag()));
+  auto at = 0.5 * atan2(z.imag(), z.real());
+  return {m * cos(at), m * sin(at)};
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/tan.hpp b/stan/math/prim/fun/tan.hpp
index 72cf84481e6..aa1aa10f872 100644
--- a/stan/math/prim/fun/tan.hpp
+++ b/stan/math/prim/fun/tan.hpp
@@ -3,7 +3,10 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/i_times.hpp>
+#include <stan/math/prim/fun/tanh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -52,6 +55,20 @@ inline auto tan(const Container& x) {
       x, [](const auto& v) { return v.array().tan(); });
 }
 
+namespace internal {
+/**
+ * Return the tangent of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return tangent of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_tan(const std::complex<V>& z) {
+  return neg_i_times(tanh(i_times(z)));
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/tanh.hpp b/stan/math/prim/fun/tanh.hpp
index 108b034fdbd..71100f8047b 100644
--- a/stan/math/prim/fun/tanh.hpp
+++ b/stan/math/prim/fun/tanh.hpp
@@ -3,7 +3,10 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/core/operator_division.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -52,6 +55,25 @@ inline auto tanh(const Container& x) {
       x, [](const auto& v) { return v.array().tanh(); });
 }
 
+namespace internal {
+/**
+ * Return the hyperbolic tangent of the complex argument.
+ *
+ * @tparam V value type of argument
+ * @param[in] z argument
+ * @return hyperbolic tangent of the argument
+ */
+template <typename V>
+inline std::complex<V> complex_tanh(const std::complex<V>& z) {
+  using std::exp;
+  using stan::math::operator/;
+  auto exp_z = exp(z);
+  auto exp_neg_z = exp(-z);
+  return stan::math::internal::complex_divide(exp_z - exp_neg_z,
+                                              exp_z + exp_neg_z);
+}
+}  // namespace internal
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/value_of_rec.hpp b/stan/math/prim/fun/value_of_rec.hpp
index c6f4fd1ab9c..8f6ee17deb6 100644
--- a/stan/math/prim/fun/value_of_rec.hpp
+++ b/stan/math/prim/fun/value_of_rec.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <complex>
 #include <cstddef>
 #include <vector>
 
@@ -42,6 +43,18 @@ inline double value_of_rec(const T x) {
  */
 inline double value_of_rec(double x) { return x; }
 
+/**
+ * Recursively apply value-of to the parts of the argument.
+ *
+ * @tparam T value type of argument
+ * @param[in] x argument
+ * @return real complex value of argument
+ */
+template <typename T>
+inline std::complex<double> value_of_rec(const std::complex<T>& x) {
+  return {value_of_rec(x.real()), value_of_rec(x.imag())};
+}
+
 /**
  * Convert a std::vector of type T to a std::vector of doubles.
  *
diff --git a/stan/math/rev/core/operator_division.hpp b/stan/math/rev/core/operator_division.hpp
index 5d5dc1e86c5..420ec04d07d 100644
--- a/stan/math/rev/core/operator_division.hpp
+++ b/stan/math/rev/core/operator_division.hpp
@@ -2,12 +2,20 @@
 #define STAN_MATH_REV_CORE_OPERATOR_DIVISION_HPP
 
 #include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/core/operator_division.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/is_any_nan.hpp>
 #include <stan/math/rev/core/var.hpp>
+#include <stan/math/rev/core/std_complex.hpp>
 #include <stan/math/rev/core/vv_vari.hpp>
 #include <stan/math/rev/core/vd_vari.hpp>
 #include <stan/math/rev/core/dv_vari.hpp>
-#include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/is_any_nan.hpp>
+#include <stan/math/rev/core/operator_addition.hpp>
+#include <stan/math/rev/core/operator_multiplication.hpp>
+#include <stan/math/rev/core/operator_subtraction.hpp>
+#include <stan/math/rev/meta/apply_scalar_unary.hpp>
+#include <complex>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -132,6 +140,11 @@ inline var operator/(Arith dividend, var divisor) {
   return {new internal::divide_dv_vari(dividend, divisor.vi_)};
 }
 
+inline std::complex<var> operator/(const std::complex<var>& x1,
+                                   const std::complex<var>& x2) {
+  return internal::complex_divide(x1, x2);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/core/operator_equal.hpp b/stan/math/rev/core/operator_equal.hpp
index a6f33b12d41..1f123ac1da4 100644
--- a/stan/math/rev/core/operator_equal.hpp
+++ b/stan/math/rev/core/operator_equal.hpp
@@ -55,6 +55,32 @@ inline bool operator==(Arith a, var b) {
   return a == b.val();
 }
 
+/**
+ * Return `true` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ *
+ * @param x real number
+ * @param z complex number
+ * @return `true` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ */
+inline bool operator==(const var& x, const std::complex<var>& z) {
+  return x == z.real() && 0 == z.imag();
+}
+
+/**
+ * Return `true` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ *
+ * @param z complex number
+ * @param y real number
+ * @return `true` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ */
+inline bool operator==(const std::complex<var>& z, const var& y) {
+  return z.real() == y && z.imag() == 0;
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/core/operator_not_equal.hpp b/stan/math/rev/core/operator_not_equal.hpp
index c5e6d71a7a7..9a5148294fe 100644
--- a/stan/math/rev/core/operator_not_equal.hpp
+++ b/stan/math/rev/core/operator_not_equal.hpp
@@ -1,8 +1,11 @@
 #ifndef STAN_MATH_REV_CORE_OPERATOR_NOT_EQUAL_HPP
 #define STAN_MATH_REV_CORE_OPERATOR_NOT_EQUAL_HPP
 
+#include <stan/math/rev/core/std_complex.hpp>
+#include <stan/math/rev/core/operator_equal.hpp>
 #include <stan/math/rev/core/var.hpp>
 #include <stan/math/prim/meta.hpp>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -56,6 +59,32 @@ inline bool operator!=(Arith a, var b) {
   return a != b.val();
 }
 
+/**
+ * Return `false` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ *
+ * @param x real number
+ * @param z complex number
+ * @return `false` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ */
+inline bool operator!=(const var& x, const std::complex<var>& z) {
+  return !(x == z.real() && 0 == z.imag());
+}
+
+/**
+ * Return `false` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ *
+ * @param z complex number
+ * @param y real number
+ * @return `false` if the real number is equal to the real part of the
+ * complex number, and the imaginary part of the complex number is zero
+ */
+inline bool operator!=(const std::complex<var>& z, const var& y) {
+  return !(z.real() == y && z.imag() == 0);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/core/std_complex.hpp b/stan/math/rev/core/std_complex.hpp
index be02b304c73..921be9a1146 100644
--- a/stan/math/rev/core/std_complex.hpp
+++ b/stan/math/rev/core/std_complex.hpp
@@ -42,8 +42,12 @@ class complex<stan::math::var>
    * @tparam U type of real part (assignable to `value_type`)
    * @param[in] re real part
    */
-  template <typename U, typename = stan::require_stan_scalar_t<U>>
-  complex(U&& re) : base_t(re) {}  // NOLINT(runtime/explicit)
+  template <typename U>  // , typename = stan::require_stan_scalar_t<U>>
+  complex(const U& re) : base_t(re) {}  // NOLINT(runtime/explicit)
+
+  template <typename U>
+  complex(const std::complex<U>& z)  // NOLINT(runtime/explicit)
+      : base_t(z.real(), z.imag()) {}
 
   /**
    * Set the real and imaginary components of this complex number to
diff --git a/stan/math/rev/core/var.hpp b/stan/math/rev/core/var.hpp
index 1eb6587e227..1f9cb114000 100644
--- a/stan/math/rev/core/var.hpp
+++ b/stan/math/rev/core/var.hpp
@@ -186,71 +186,7 @@ class var {
   var(unsigned long x)
       : vi_(new vari(static_cast<double>(x), false)) {}  // NOLINT
 
-  /**
-   * Construct a variable from the specified arithmetic argument
-   * by constructing a new <code>vari</code> with the argument
-   * cast to <code>double</code>, and a zero adjoint. Only works
-   * if the imaginary part is zero.
-   *
-   * @param x Value of the variable.
-   */
-  explicit var(const std::complex<double>& x) {
-    if (imag(x) == 0) {
-      vi_ = new vari(real(x), false);
-    } else {
-      std::stringstream ss;
-      ss << "Imaginary part of std::complex used to construct var"
-         << " must be zero. Found real part = " << real(x) << " and "
-         << " found imaginary part = " << imag(x) << std::endl;
-      std::string msg = ss.str();
-      throw std::invalid_argument(msg);
-    }
-  }
-
-  /**
-   * Construct a variable from the specified arithmetic argument
-   * by constructing a new <code>vari</code> with the argument
-   * cast to <code>double</code>, and a zero adjoint. Only works
-   * if the imaginary part is zero.
-   *
-   * @param x Value of the variable.
-   */
-  explicit var(const std::complex<float>& x) {
-    if (imag(x) == 0) {
-      vi_ = new vari(static_cast<double>(real(x)), false);
-    } else {
-      std::stringstream ss;
-      ss << "Imaginary part of std::complex used to construct var"
-         << " must be zero. Found real part = " << real(x) << " and "
-         << " found imaginary part = " << imag(x) << std::endl;
-      std::string msg = ss.str();
-      throw std::invalid_argument(msg);
-    }
-  }
-
-  /**
-   * Construct a variable from the specified arithmetic argument
-   * by constructing a new <code>vari</code> with the argument
-   * cast to <code>double</code>, and a zero adjoint. Only works
-   * if the imaginary part is zero.
-   *
-   * @param x Value of the variable.
-   */
-  explicit var(const std::complex<long double>& x) {
-    if (imag(x) == 0) {
-      vi_ = new vari(static_cast<double>(real(x)), false);
-    } else {
-      std::stringstream ss;
-      ss << "Imaginary part of std::complex used to construct var"
-         << " must be zero. Found real part = " << real(x) << " and "
-         << " found imaginary part = " << imag(x) << std::endl;
-      std::string msg = ss.str();
-      throw std::invalid_argument(msg);
-    }
-  }
-
 #ifdef _WIN64
-
   // these two ctors are for Win64 to enable 64-bit signed
   // and unsigned integers, because long and unsigned long
   // are still 32-bit
diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp
index 1dee72acd23..a6025090f21 100644
--- a/stan/math/rev/fun.hpp
+++ b/stan/math/rev/fun.hpp
@@ -13,6 +13,7 @@
 #include <stan/math/rev/fun/acos.hpp>
 #include <stan/math/rev/fun/acosh.hpp>
 #include <stan/math/rev/fun/as_bool.hpp>
+#include <stan/math/rev/fun/arg.hpp>
 #include <stan/math/rev/fun/asin.hpp>
 #include <stan/math/rev/fun/asinh.hpp>
 #include <stan/math/rev/fun/atan.hpp>
@@ -27,6 +28,7 @@
 #include <stan/math/rev/fun/cholesky_decompose.hpp>
 #include <stan/math/rev/fun/columns_dot_product.hpp>
 #include <stan/math/rev/fun/columns_dot_self.hpp>
+#include <stan/math/rev/fun/conj.hpp>
 #include <stan/math/rev/fun/cos.hpp>
 #include <stan/math/rev/fun/cosh.hpp>
 #include <stan/math/rev/fun/cov_exp_quad.hpp>
@@ -102,11 +104,14 @@
 #include <stan/math/rev/fun/multiply.hpp>
 #include <stan/math/rev/fun/multiply_log.hpp>
 #include <stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp>
+#include <stan/math/rev/fun/norm.hpp>
 #include <stan/math/rev/fun/ordered_constrain.hpp>
 #include <stan/math/rev/fun/owens_t.hpp>
+#include <stan/math/rev/fun/polar.hpp>
 #include <stan/math/rev/fun/positive_ordered_constrain.hpp>
 #include <stan/math/rev/fun/pow.hpp>
 #include <stan/math/rev/fun/primitive_value.hpp>
+#include <stan/math/rev/fun/proj.hpp>
 #include <stan/math/rev/fun/quad_form.hpp>
 #include <stan/math/rev/fun/quad_form_sym.hpp>
 #include <stan/math/rev/fun/rising_factorial.hpp>
diff --git a/stan/math/rev/fun/abs.hpp b/stan/math/rev/fun/abs.hpp
index 9aad50f2865..d77abd3e09e 100644
--- a/stan/math/rev/fun/abs.hpp
+++ b/stan/math/rev/fun/abs.hpp
@@ -1,8 +1,11 @@
 #ifndef STAN_MATH_REV_FUN_ABS_HPP
 #define STAN_MATH_REV_FUN_ABS_HPP
 
+#include <stan/math/prim/fun/abs.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/fun/fabs.hpp>
+#include <stan/math/rev/fun/hypot.hpp>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -35,6 +38,16 @@ namespace math {
  */
 inline var abs(const var& a) { return fabs(a); }
 
+/**
+ * Return the absolute value of the complex argument.
+ *
+ * @param[in] z argument
+ * @return absolute value of the argument
+ */
+inline std::complex<var> abs(const std::complex<var>& z) {
+  return internal::complex_abs(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/acos.hpp b/stan/math/rev/fun/acos.hpp
index f6510e46191..909d1a73ab6 100644
--- a/stan/math/rev/fun/acos.hpp
+++ b/stan/math/rev/fun/acos.hpp
@@ -1,9 +1,17 @@
 #ifndef STAN_MATH_REV_FUN_ACOS_HPP
 #define STAN_MATH_REV_FUN_ACOS_HPP
 
-#include <stan/math/rev/meta.hpp>
+#include <stan/math/prim/fun/acos.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/arg.hpp>
+#include <stan/math/rev/fun/asin.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
+#include <stan/math/rev/fun/polar.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -56,6 +64,16 @@ class acos_vari : public op_v_vari {
  */
 inline var acos(const var& a) { return var(new internal::acos_vari(a.vi_)); }
 
+/**
+ * Return the arc cosine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return arc cosine of the argument
+ */
+inline std::complex<var> acos(const std::complex<var>& z) {
+  return stan::math::internal::complex_acos(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/acosh.hpp b/stan/math/rev/fun/acosh.hpp
index 3d4833445e5..dfc3b27bdab 100644
--- a/stan/math/rev/fun/acosh.hpp
+++ b/stan/math/rev/fun/acosh.hpp
@@ -1,10 +1,22 @@
 #ifndef STAN_MATH_REV_FUN_ACOSH_HPP
 #define STAN_MATH_REV_FUN_ACOSH_HPP
 
-#include <stan/math/rev/meta.hpp>
-#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/core.hpp>
+#include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/acosh.hpp>
+#include <stan/math/prim/fun/isnan.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/value_of_rec.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/arg.hpp>
+#include <stan/math/rev/fun/cosh.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
+#include <stan/math/rev/fun/log.hpp>
+#include <stan/math/rev/fun/polar.hpp>
+#include <stan/math/rev/fun/sqrt.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -62,6 +74,16 @@ inline var acosh(const var& a) {
   return var(new internal::acosh_vari(acosh(a.val()), a.vi_));
 }
 
+/**
+ * Return the hyperbolic arc cosine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return hyperbolic arc cosine of the argument
+ */
+inline std::complex<var> acosh(const std::complex<var>& z) {
+  return stan::math::internal::complex_acosh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/arg.hpp b/stan/math/rev/fun/arg.hpp
new file mode 100644
index 00000000000..293dfd6b059
--- /dev/null
+++ b/stan/math/rev/fun/arg.hpp
@@ -0,0 +1,23 @@
+#ifndef STAN_MATH_REV_FUN_ARG_HPP
+#define STAN_MATH_REV_FUN_ARG_HPP
+
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/atan2.hpp>
+#include <stan/math/prim/fun/arg.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the phase angle of the complex argument.
+ *
+ * @param[in] z argument
+ * @return phase angle of the argument
+ */
+inline var arg(const std::complex<var>& z) { return internal::complex_arg(z); }
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/fun/asin.hpp b/stan/math/rev/fun/asin.hpp
index 64e63d0a5fb..ed03e449ec8 100644
--- a/stan/math/rev/fun/asin.hpp
+++ b/stan/math/rev/fun/asin.hpp
@@ -1,9 +1,15 @@
 #ifndef STAN_MATH_REV_FUN_ASIN_HPP
 #define STAN_MATH_REV_FUN_ASIN_HPP
 
-#include <stan/math/rev/meta.hpp>
+#include <stan/math/prim/fun/asin.hpp>
+#include <stan/math/prim/fun/abs.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/asinh.hpp>
+#include <stan/math/rev/fun/value_of_rec.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -56,6 +62,16 @@ class asin_vari : public op_v_vari {
  */
 inline var asin(const var& a) { return var(new internal::asin_vari(a.vi_)); }
 
+/**
+ * Return the arc sine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return arc sine of the argument
+ */
+inline std::complex<var> asin(const std::complex<var>& z) {
+  return stan::math::internal::complex_asin(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/asinh.hpp b/stan/math/rev/fun/asinh.hpp
index bbcfbb9f78a..5ff296a19de 100644
--- a/stan/math/rev/fun/asinh.hpp
+++ b/stan/math/rev/fun/asinh.hpp
@@ -1,10 +1,21 @@
 #ifndef STAN_MATH_REV_FUN_ASINH_HPP
 #define STAN_MATH_REV_FUN_ASINH_HPP
 
-#include <stan/math/rev/meta.hpp>
-#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/core.hpp>
+#include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/asinh.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/value_of_rec.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/arg.hpp>
+#include <stan/math/rev/fun/cosh.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
+#include <stan/math/rev/fun/log.hpp>
+#include <stan/math/rev/fun/polar.hpp>
+#include <stan/math/rev/fun/sqrt.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -56,6 +67,16 @@ inline var asinh(const var& a) {
   return var(new internal::asinh_vari(asinh(a.val()), a.vi_));
 }
 
+/**
+ * Return the hyperbolic arcsine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return hyperbolic arcsine of the argument
+ */
+inline std::complex<var> asinh(const std::complex<var>& z) {
+  return stan::math::internal::complex_asinh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/atan.hpp b/stan/math/rev/fun/atan.hpp
index 8a62d347fa3..9d4ead191fc 100644
--- a/stan/math/rev/fun/atan.hpp
+++ b/stan/math/rev/fun/atan.hpp
@@ -1,9 +1,22 @@
 #ifndef STAN_MATH_REV_FUN_ATAN_HPP
 #define STAN_MATH_REV_FUN_ATAN_HPP
 
-#include <stan/math/rev/meta.hpp>
+#include <stan/math/prim/fun/abs.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/isnan.hpp>
+#include <stan/math/prim/fun/atan.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/hypot.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/cos.hpp>
+#include <stan/math/rev/fun/atanh.hpp>
+#include <stan/math/rev/fun/log.hpp>
+#include <stan/math/rev/fun/sqrt.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -48,6 +61,16 @@ class atan_vari : public op_v_vari {
  */
 inline var atan(const var& a) { return var(new internal::atan_vari(a.vi_)); }
 
+/**
+ * Return the arc tangent of the complex argument.
+ *
+ * @param[in] z argument
+ * @return arc tangent of the argument
+ */
+inline std::complex<var> atan(const std::complex<var>& z) {
+  return stan::math::internal::complex_atan(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/atanh.hpp b/stan/math/rev/fun/atanh.hpp
index e1c0383f68d..b7a5b66c63f 100644
--- a/stan/math/rev/fun/atanh.hpp
+++ b/stan/math/rev/fun/atanh.hpp
@@ -1,9 +1,17 @@
 #ifndef STAN_MATH_REV_FUN_ATANH_HPP
 #define STAN_MATH_REV_FUN_ATANH_HPP
 
+#include <stan/math/prim/fun/atanh.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
-#include <stan/math/prim/fun/atanh.hpp>
+#include <stan/math/rev/fun/value_of_rec.hpp>
+#include <stan/math/rev/fun/atan2.hpp>
+#include <stan/math/rev/fun/cosh.hpp>
+#include <stan/math/rev/fun/log.hpp>
+#include <stan/math/rev/fun/sinh.hpp>
+#include <stan/math/rev/fun/hypot.hpp>
+#include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -59,6 +67,16 @@ inline var atanh(const var& a) {
   return var(new internal::atanh_vari(atanh(a.val()), a.vi_));
 }
 
+/**
+ * Return the hyperbolic arc tangent of the complex argument.
+ *
+ * @param[in] z argument
+ * @return hyperbolic arc tangent of the argument
+ */
+inline std::complex<var> atanh(const std::complex<var>& z) {
+  return stan::math::internal::complex_atanh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/conj.hpp b/stan/math/rev/fun/conj.hpp
new file mode 100644
index 00000000000..2537f17c8e8
--- /dev/null
+++ b/stan/math/rev/fun/conj.hpp
@@ -0,0 +1,24 @@
+#ifndef STAN_MATH_REV_FUN_CONJ_HPP
+#define STAN_MATH_REV_FUN_CONJ_HPP
+
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/fun/conj.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the complex conjugate of the complex argument.
+ *
+ * @param[in] z argument
+ * @return complex conjugate of the argument
+ */
+inline std::complex<var> conj(const std::complex<var>& z) {
+  return internal::complex_conj(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/fun/cos.hpp b/stan/math/rev/fun/cos.hpp
index 7256f8ab23c..da4c8067e32 100644
--- a/stan/math/rev/fun/cos.hpp
+++ b/stan/math/rev/fun/cos.hpp
@@ -1,9 +1,17 @@
 #ifndef STAN_MATH_REV_FUN_COS_HPP
 #define STAN_MATH_REV_FUN_COS_HPP
 
-#include <stan/math/rev/meta.hpp>
+#include <stan/math/prim/fun/abs.hpp>
+#include <stan/math/prim/fun/cos.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/cosh.hpp>
+#include <stan/math/rev/fun/sinh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -45,6 +53,16 @@ class cos_vari : public op_v_vari {
  */
 inline var cos(const var& a) { return var(new internal::cos_vari(a.vi_)); }
 
+/**
+ * Return the cosine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return cosine of the argument
+ */
+inline std::complex<var> cos(const std::complex<var>& z) {
+  return stan::math::internal::complex_cos(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/cosh.hpp b/stan/math/rev/fun/cosh.hpp
index 4f1c377fa67..0b3605e8c9b 100644
--- a/stan/math/rev/fun/cosh.hpp
+++ b/stan/math/rev/fun/cosh.hpp
@@ -1,9 +1,15 @@
 #ifndef STAN_MATH_REV_FUN_COSH_HPP
 #define STAN_MATH_REV_FUN_COSH_HPP
 
+#include <stan/math/prim/fun/cosh.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/core/operator_multiplication.hpp>
+#include <stan/math/rev/fun/exp.hpp>
+#include <stan/math/rev/fun/sin.hpp>
+#include <stan/math/rev/fun/sinh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -45,6 +51,16 @@ class cosh_vari : public op_v_vari {
  */
 inline var cosh(const var& a) { return var(new internal::cosh_vari(a.vi_)); }
 
+/**
+ * Return the hyperbolic cosine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return hyperbolic cosine of the argument
+ */
+inline std::complex<var> cosh(const std::complex<var>& z) {
+  return stan::math::internal::complex_cosh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/exp.hpp b/stan/math/rev/fun/exp.hpp
index 0136194f894..4a14a4811b1 100644
--- a/stan/math/rev/fun/exp.hpp
+++ b/stan/math/rev/fun/exp.hpp
@@ -3,7 +3,15 @@
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/prim/fun/cos.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
+#include <stan/math/rev/fun/sin.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -40,6 +48,15 @@ class exp_vari : public op_v_vari {
  */
 inline var exp(const var& a) { return var(new internal::exp_vari(a.vi_)); }
 
+/**
+ * Return the exponentiation (base e) of the specified complex number.
+ * @param z argument
+ * @return exponentiation of argument
+ */
+inline std::complex<var> exp(const std::complex<var>& z) {
+  return internal::complex_exp(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/log.hpp b/stan/math/rev/fun/log.hpp
index 2ab780ab289..52b891da2ec 100644
--- a/stan/math/rev/fun/log.hpp
+++ b/stan/math/rev/fun/log.hpp
@@ -1,8 +1,18 @@
 #ifndef STAN_MATH_REV_FUN_LOG_HPP
 #define STAN_MATH_REV_FUN_LOG_HPP
 
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/is_inf.hpp>
+#include <stan/math/prim/fun/log.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/abs.hpp>
+#include <stan/math/rev/fun/arg.hpp>
+#include <stan/math/rev/fun/atan2.hpp>
+#include <stan/math/rev/fun/cos.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
+#include <stan/math/rev/fun/sqrt.hpp>
 #include <cmath>
 
 namespace stan {
@@ -46,6 +56,16 @@ class log_vari : public op_v_vari {
  */
 inline var log(const var& a) { return var(new internal::log_vari(a.vi_)); }
 
+/**
+ * Return the natural logarithm (base e) of the specified complex argument.
+ *
+ * @param z complex argument
+ * @return natural logarithm of argument
+ */
+inline std::complex<var> log(const std::complex<var>& z) {
+  return internal::complex_log(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/log10.hpp b/stan/math/rev/fun/log10.hpp
index e5d246f6908..460e45f89b9 100644
--- a/stan/math/rev/fun/log10.hpp
+++ b/stan/math/rev/fun/log10.hpp
@@ -1,10 +1,16 @@
 #ifndef STAN_MATH_REV_FUN_LOG10_HPP
 #define STAN_MATH_REV_FUN_LOG10_HPP
 
+#include <stan/math/prim/core.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/log10.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
-#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/rev/fun/atan2.hpp>
+#include <stan/math/rev/fun/hypot.hpp>
+#include <stan/math/rev/fun/log.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -50,6 +56,16 @@ class log10_vari : public op_v_vari {
  */
 inline var log10(const var& a) { return var(new internal::log10_vari(a.vi_)); }
 
+/**
+ * Return the base 10 logarithm of the specified complex number.
+ *
+ * @param z complex argument
+ * @return base 10 log of argument
+ */
+inline std::complex<var> log10(const std::complex<var>& z) {
+  return internal::complex_log10(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/norm.hpp b/stan/math/rev/fun/norm.hpp
new file mode 100644
index 00000000000..41a6b6af10c
--- /dev/null
+++ b/stan/math/rev/fun/norm.hpp
@@ -0,0 +1,24 @@
+#ifndef STAN_MATH_REV_FUN_NORM_HPP
+#define STAN_MATH_REV_FUN_NORM_HPP
+
+#include <stan/math/prim/fun/norm.hpp>
+#include <stan/math/rev/core.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the squared magnitude of the complex argument.
+ *
+ * @param[in] z argument
+ * @return squared magnitude of the argument
+ */
+inline var norm(const std::complex<var>& z) {
+  return internal::complex_norm(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/fun/polar.hpp b/stan/math/rev/fun/polar.hpp
new file mode 100644
index 00000000000..d2c73bf835b
--- /dev/null
+++ b/stan/math/rev/fun/polar.hpp
@@ -0,0 +1,59 @@
+#ifndef STAN_MATH_REV_FUN_POLAR_HPP
+#define STAN_MATH_REV_FUN_POLAR_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/polar.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
+#include <stan/math/prim/fun/is_inf.hpp>
+#include <stan/math/prim/fun/sin.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/rev/fun/sin.hpp>
+#include <complex>
+#include <type_traits>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+inline std::complex<var> polar(const var& r, const var& theta) {
+  return internal::complex_polar(r, theta);
+}
+
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @tparam T arithmetic type of magnitude
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+template <typename T>
+inline std::complex<var> polar(T r, const var& theta) {
+  return internal::complex_polar(r, theta);
+}
+
+/**
+ * Returns complex number with specified magnitude and phase angle.
+ *
+ * @tparam T arithmetic type of phase angle
+ * @param[in] r magnitude
+ * @param[in] theta phase angle
+ * @return complex number with magnitude and phase angle
+ */
+template <typename T>
+inline std::complex<var> polar(const var& r, T theta) {
+  return internal::complex_polar(r, theta);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/fun/pow.hpp b/stan/math/rev/fun/pow.hpp
index 1b401abda32..577c23381fd 100644
--- a/stan/math/rev/fun/pow.hpp
+++ b/stan/math/rev/fun/pow.hpp
@@ -1,16 +1,27 @@
 #ifndef STAN_MATH_REV_FUN_POW_HPP
 #define STAN_MATH_REV_FUN_POW_HPP
 
+#include <stan/math/prim/core.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/copysign.hpp>
+#include <stan/math/prim/fun/is_any_nan.hpp>
+#include <stan/math/prim/fun/isnan.hpp>
+#include <stan/math/prim/fun/is_nan.hpp>
+#include <stan/math/prim/fun/pow.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/inv.hpp>
 #include <stan/math/rev/fun/inv_sqrt.hpp>
 #include <stan/math/rev/fun/inv_square.hpp>
+#include <stan/math/rev/fun/is_nan.hpp>
+#include <stan/math/rev/fun/log.hpp>
 #include <stan/math/rev/fun/sqrt.hpp>
 #include <stan/math/rev/fun/square.hpp>
-#include <stan/math/prim/fun/constants.hpp>
-#include <stan/math/prim/fun/is_any_nan.hpp>
+#include <stan/math/rev/fun/value_of_rec.hpp>
 #include <cmath>
+#include <complex>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -120,15 +131,13 @@ inline var pow(const var& base, const var& exponent) {
  * The template parameters are coded as they are so that arithmetic
  * types will not be promoted into the `var` slots.
  *
- * @tparam Var var type
- * @tparam Arith arithmetic type
+ * @tparam T arithmetic type
  * @param base Base variable.
  * @param exponent Exponent scalar.
  * @return Base raised to the exponent.
  */
-template <typename Var, typename Arith, require_var_t<Var>...,
-          require_arithmetic_t<Arith>...>
-inline var pow(const Var& base, Arith exponent) {
+template <typename T, typename = require_arithmetic_t<T>>
+inline var pow(const var& base, T exponent) {
   if (exponent == 0.5) {
     return sqrt(base);
   }
@@ -161,19 +170,157 @@ inline var pow(const Var& base, Arith exponent) {
  * The template parameters are coded as they are so that arithmetic
  * types will not be promoted into the `var` slots.
  *
- * @tparam Var var type
- * @tparam Arith arithmetic type
+ * @tparam T arithmetic type
  *
  * @param base Base scalar.
  * @param exponent Exponent variable.
  * @return Base raised to the exponent.
  */
-template <typename Arith, typename Var, require_arithmetic_t<Arith>...,
-          require_var_t<Var>>
-inline var pow(Arith base, const Var& exponent) {
+template <typename T, typename = require_arithmetic_t<T>>
+inline var pow(T base, const var& exponent) {
   return {new internal::pow_dv_vari(base, exponent.vi_)};
 }
 
+// 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.
+ *
+ * @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, int y) {
+  return internal::complex_pow(x, y);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/proj.hpp b/stan/math/rev/fun/proj.hpp
new file mode 100644
index 00000000000..ceb83db8bfa
--- /dev/null
+++ b/stan/math/rev/fun/proj.hpp
@@ -0,0 +1,27 @@
+#ifndef STAN_MATH_REV_FUN_PROJ_HPP
+#define STAN_MATH_REV_FUN_PROJ_HPP
+
+#include <stan/math/prim/fun/proj.hpp>
+#include <stan/math/prim/fun/is_inf.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the projection of the complex argument onto the Riemann
+ * sphere.
+ *
+ * @param[in] z argument
+ * @return projection of the argument onto the Riemann sphere
+ */
+inline std::complex<var> proj(const std::complex<var>& z) {
+  return internal::complex_proj(z);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/fun/sin.hpp b/stan/math/rev/fun/sin.hpp
index 9a478e7b731..61dfc54d9af 100644
--- a/stan/math/rev/fun/sin.hpp
+++ b/stan/math/rev/fun/sin.hpp
@@ -1,9 +1,17 @@
 #ifndef STAN_MATH_REV_FUN_SIN_HPP
 #define STAN_MATH_REV_FUN_SIN_HPP
 
-#include <stan/math/rev/meta.hpp>
+#include <stan/math/prim/fun/cos.hpp>
+#include <stan/math/prim/fun/isfinite.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/sin.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/rev/fun/cosh.hpp>
+#include <stan/math/rev/fun/sinh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -45,6 +53,16 @@ class sin_vari : public op_v_vari {
  */
 inline var sin(const var& a) { return var(new internal::sin_vari(a.vi_)); }
 
+/**
+ * Return the sine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return sine of the argument
+ */
+inline std::complex<var> sin(const std::complex<var>& z) {
+  return stan::math::internal::complex_sin(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/sinh.hpp b/stan/math/rev/fun/sinh.hpp
index 822b70940c2..fecc1f431fc 100644
--- a/stan/math/rev/fun/sinh.hpp
+++ b/stan/math/rev/fun/sinh.hpp
@@ -1,9 +1,15 @@
 #ifndef STAN_MATH_REV_FUN_SINH_HPP
 #define STAN_MATH_REV_FUN_SINH_HPP
 
+#include <stan/math/prim/core.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/sinh.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/exp.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -45,6 +51,16 @@ class sinh_vari : public op_v_vari {
  */
 inline var sinh(const var& a) { return var(new internal::sinh_vari(a.vi_)); }
 
+/**
+ * Return the hyperbolic sine of the complex argument.
+ *
+ * @param[in] z argument
+ * @return hyperbolic sine of the argument
+ */
+inline std::complex<var> sinh(const std::complex<var>& z) {
+  return stan::math::internal::complex_sinh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/sqrt.hpp b/stan/math/rev/fun/sqrt.hpp
index 7b8bddc627e..5dafb4dac4d 100644
--- a/stan/math/rev/fun/sqrt.hpp
+++ b/stan/math/rev/fun/sqrt.hpp
@@ -1,9 +1,14 @@
 #ifndef STAN_MATH_REV_FUN_SQRT_HPP
 #define STAN_MATH_REV_FUN_SQRT_HPP
 
+#include <stan/math/prim/fun/sqrt.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/atan2.hpp>
+#include <stan/math/rev/fun/cos.hpp>
+#include <stan/math/rev/fun/hypot.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -46,6 +51,16 @@ class sqrt_vari : public op_v_vari {
  */
 inline var sqrt(const var& a) { return var(new internal::sqrt_vari(a.vi_)); }
 
+/**
+ * Return the square root of the complex argument.
+ *
+ * @param[in] z argument
+ * @return square root of the argument
+ */
+inline std::complex<var> sqrt(const std::complex<var>& z) {
+  return internal::complex_sqrt(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/tan.hpp b/stan/math/rev/fun/tan.hpp
index 28f03b6eb5e..e7f231b28ae 100644
--- a/stan/math/rev/fun/tan.hpp
+++ b/stan/math/rev/fun/tan.hpp
@@ -1,9 +1,15 @@
 #ifndef STAN_MATH_REV_FUN_TAN_HPP
 #define STAN_MATH_REV_FUN_TAN_HPP
 
+#include <stan/math/prim/fun/isinf.hpp>
+#include <stan/math/prim/fun/tan.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/is_inf.hpp>
+#include <stan/math/rev/fun/sinh.hpp>
+#include <stan/math/rev/fun/tanh.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -45,6 +51,16 @@ class tan_vari : public op_v_vari {
  */
 inline var tan(const var& a) { return var(new internal::tan_vari(a.vi_)); }
 
+/**
+ * Return the tangent of the complex argument.
+ *
+ * @param[in] z argument
+ * @return tangent of the argument
+ */
+inline std::complex<var> tan(const std::complex<var>& z) {
+  return stan::math::internal::complex_tan(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/fun/tanh.hpp b/stan/math/rev/fun/tanh.hpp
index dde1ec7f9a8..ef9763989e2 100644
--- a/stan/math/rev/fun/tanh.hpp
+++ b/stan/math/rev/fun/tanh.hpp
@@ -1,9 +1,12 @@
 #ifndef STAN_MATH_REV_FUN_TANH_HPP
 #define STAN_MATH_REV_FUN_TANH_HPP
 
+#include <stan/math/prim/fun/tanh.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/exp.hpp>
 #include <cmath>
+#include <complex>
 
 namespace stan {
 namespace math {
@@ -48,6 +51,16 @@ class tanh_vari : public op_v_vari {
  */
 inline var tanh(const var& a) { return var(new internal::tanh_vari(a.vi_)); }
 
+/**
+ * Return the hyperbolic tangent of the complex argument.
+ *
+ * @param[in] z argument
+ * @return hyperbolic tangent of the argument
+ */
+inline std::complex<var> tanh(const std::complex<var>& z) {
+  return stan::math::internal::complex_tanh(z);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/test/unit/math/expect_near_rel.hpp b/test/unit/math/expect_near_rel.hpp
index 9f7f7ae0696..f3f7e645331 100644
--- a/test/unit/math/expect_near_rel.hpp
+++ b/test/unit/math/expect_near_rel.hpp
@@ -122,6 +122,17 @@ void expect_near_rel(const std::string& msg, EigMat1&& x1, EigMat2&& x2,
     expect_near_rel(msg2, x1_eval(i), x2_eval(i), tol);
 }
 
+/**
+ * Tewsts that the elements of the specified standard vectors are
+ * relatively near one another by calling `expect_near_rel`
+ * recursively.
+ *
+ * @tparam T1 value type for first vector
+ * @tparam T2 value type for second vector
+ * @param[in] x1 first vector
+ * @param[in] x2 second vector
+ * @param[in] tol relative tolerance
+ */
 template <typename T1, typename T2>
 void expect_near_rel(const std::string& msg, const std::vector<T1>& x1,
                      const std::vector<T2>& x2,
@@ -136,6 +147,68 @@ void expect_near_rel(const std::string& msg, const std::vector<T1>& x1,
     expect_near_rel(msg2, x1[i], x2[i], tol);
 }
 
+/**
+ * Tests that the real and complex parts of the specified complex numbers
+ * by calling `expect_near_rel` recursively with the specified message
+ * and tolerance.
+ *
+ * @tparam T1 value type of first complex number
+ * @tparam T2 value type of second complex number
+ * @param msg[in] message to print under failure
+ * @param z1[in] first complex number
+ * @param z2[in] second complex number
+ * @param tol[in] tolerance for comparison
+ */
+template <typename T1, typename T2>
+void expect_near_rel(const std::string& msg, const std::complex<T1>& z1,
+                     const std::complex<T2>& z2,
+                     relative_tolerance tol = relative_tolerance()) {
+  expect_near_rel(msg, z1.real(), z2.real(), tol);
+  expect_near_rel(msg, z1.imag(), z2.imag(), tol);
+}
+
+/**
+ * Tests that the specified real number is near the real part of the
+ * specified complex number and the complex number's imaginary part is
+ * near zero by calling `expect_near_rel` recursively with the
+ * specified message and tolerance.
+ *
+ * @tparam T1 value type of first number
+ * @tparam T2 value type of second complex number
+ * @param msg[in] message to print under failure
+ * @param x1[in] real number
+ * @param z2[in] complex number
+ * @param tol[in] tolerance for comparison
+ */
+template <typename T1, typename T2>
+void expect_near_rel(const std::string& msg, const T1& x1,
+                     const std::complex<T2>& z2,
+                     relative_tolerance tol = relative_tolerance()) {
+  expect_near_rel(msg, x1, z2.real(), tol);
+  expect_near_rel(msg, 0, z2.imag(), tol);
+}
+
+/**
+ * Tests that the specified real number is near the real part of the
+ * specified complex number and the complex number's imaginary part is
+ * near zero by calling `expect_near_rel` recursively with the
+ * specified message and tolerance.
+ *
+ * @tparam T1 value type of first number
+ * @tparam T2 value type of second complex number
+ * @param msg[in] message to print under failure
+ * @param z1[in] complex number
+ * @param x2[in] real number
+ * @param tol[in] tolerance for comparison
+ */
+template <typename T1, typename T2>
+void expect_near_rel(const std::string& msg, const std::complex<T1>& z1,
+                     const T2& x2,
+                     relative_tolerance tol = relative_tolerance()) {
+  expect_near_rel(msg, z1.real(), x2, tol);
+  expect_near_rel(msg, z1.imag(), 0, tol);
+}
+
 }  // namespace test
 }  // namespace stan
 #endif
diff --git a/test/unit/math/mix/core/operator_addition_test.cpp b/test/unit/math/mix/core/operator_addition_test.cpp
index 5fce58a166d..62bb8d0f189 100644
--- a/test/unit/math/mix/core/operator_addition_test.cpp
+++ b/test/unit/math/mix/core/operator_addition_test.cpp
@@ -1,7 +1,10 @@
+#include <stan/math/mix.hpp>
+#include <stan/math/prim/core/operator_addition.hpp>
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixCore, operatorAddition) {
   auto f = [](const auto& x1, const auto& x2) { return x1 + x2; };
   bool disable_lhs_int = true;
   stan::test::expect_common_binary(f, disable_lhs_int);
+  stan::test::expect_complex_common_binary(f);
 }
diff --git a/test/unit/math/mix/core/operator_division_test.cpp b/test/unit/math/mix/core/operator_division_test.cpp
index d61ec16afba..06bf2e9b291 100644
--- a/test/unit/math/mix/core/operator_division_test.cpp
+++ b/test/unit/math/mix/core/operator_division_test.cpp
@@ -1,7 +1,35 @@
 #include <test/unit/math/test_ad.hpp>
+#include <vector>
 
 TEST(mathMixCore, operatorDivision) {
   auto f = [](const auto& x1, const auto& x2) { return x1 / x2; };
   bool disable_lhs_int = true;
   stan::test::expect_common_binary(f, disable_lhs_int);
+
+  std::vector<double> common_finite = {-2.9, -1, -0.0, 0.0, 1, 1.39};
+  std::vector<double> common_finite_nz = {-3.1, -1, 1, 2.7};
+  for (auto re1 : common_finite) {
+    for (auto im1 : common_finite_nz) {
+      for (auto re2 : common_finite) {
+        for (auto im2 : common_finite_nz) {
+          stan::test::expect_ad(f, std::complex<double>(re1, im1),
+                                std::complex<double>(re2, im2));
+        }
+      }
+    }
+  }
+  for (auto re1 : common_finite) {
+    for (auto re2 : common_finite) {
+      for (auto im2 : common_finite_nz) {
+        stan::test::expect_ad(f, re1, std::complex<double>(re2, im2));
+      }
+    }
+  }
+  for (auto re1 : common_finite) {
+    for (auto im1 : common_finite_nz) {
+      for (auto re2 : common_finite) {
+        stan::test::expect_ad(f, std::complex<double>(re1, im1), re2);
+      }
+    }
+  }
 }
diff --git a/test/unit/math/mix/core/operator_equal_test.cpp b/test/unit/math/mix/core/operator_equal_test.cpp
index 9c425a468e9..d3ebdace1b7 100644
--- a/test/unit/math/mix/core/operator_equal_test.cpp
+++ b/test/unit/math/mix/core/operator_equal_test.cpp
@@ -3,4 +3,5 @@
 TEST(mathMixCore, operatorEqual) {
   auto f = [](const auto& x1, const auto& x2) { return x1 == x2; };
   stan::test::expect_common_comparison(f);
+  stan::test::expect_complex_common_comparison(f);
 }
diff --git a/test/unit/math/mix/core/operator_multiplication_test.cpp b/test/unit/math/mix/core/operator_multiplication_test.cpp
index c752957c5fc..32b553f30d2 100644
--- a/test/unit/math/mix/core/operator_multiplication_test.cpp
+++ b/test/unit/math/mix/core/operator_multiplication_test.cpp
@@ -4,4 +4,5 @@ TEST(mathMixCore, operatorMultiplication) {
   auto f = [](const auto& x1, const auto& x2) { return x1 * x2; };
   bool disable_lhs_int = true;
   stan::test::expect_common_binary(f, disable_lhs_int);
+  stan::test::expect_complex_common_binary(f);
 }
diff --git a/test/unit/math/mix/core/operator_not_equal_test.cpp b/test/unit/math/mix/core/operator_not_equal_test.cpp
index 37b998fba20..0289a031f7c 100644
--- a/test/unit/math/mix/core/operator_not_equal_test.cpp
+++ b/test/unit/math/mix/core/operator_not_equal_test.cpp
@@ -3,4 +3,5 @@
 TEST(mathMixCore, operatorNotEqual) {
   auto f = [](const auto& x1, const auto& x2) { return x1 != x2; };
   stan::test::expect_common_comparison(f);
+  stan::test::expect_complex_common_comparison(f);
 }
diff --git a/test/unit/math/mix/core/operator_subtraction_test.cpp b/test/unit/math/mix/core/operator_subtraction_test.cpp
index c657571210b..0cb0efd11cb 100644
--- a/test/unit/math/mix/core/operator_subtraction_test.cpp
+++ b/test/unit/math/mix/core/operator_subtraction_test.cpp
@@ -4,4 +4,5 @@ TEST(mathMixCore, operatorSubtraction) {
   auto f = [](const auto& x1, const auto& x2) { return x1 - x2; };
   bool disable_lhs_int = true;
   stan::test::expect_common_binary(f, disable_lhs_int);
+  stan::test::expect_complex_common_binary(f);
 }
diff --git a/test/unit/math/mix/core/operator_unary_minus_test.cpp b/test/unit/math/mix/core/operator_unary_minus_test.cpp
index 80e5b242d74..6e650963f74 100644
--- a/test/unit/math/mix/core/operator_unary_minus_test.cpp
+++ b/test/unit/math/mix/core/operator_unary_minus_test.cpp
@@ -3,4 +3,5 @@
 TEST(mathMixCore, operatorUnaryMinus) {
   auto f = [](const auto& x1) { return -x1; };
   stan::test::expect_common_unary(f);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/core/operator_unary_plus_test.cpp b/test/unit/math/mix/core/operator_unary_plus_test.cpp
index a61674f3f5b..450fa549da7 100644
--- a/test/unit/math/mix/core/operator_unary_plus_test.cpp
+++ b/test/unit/math/mix/core/operator_unary_plus_test.cpp
@@ -3,4 +3,5 @@
 TEST(mathMixCore, operatorUnaryPlus) {
   auto f = [](const auto& x1) { return +x1; };
   stan::test::expect_common_unary(f);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/abs_test.cpp b/test/unit/math/mix/fun/abs_test.cpp
index 19d43f75b00..11f5da342b2 100644
--- a/test/unit/math/mix/fun/abs_test.cpp
+++ b/test/unit/math/mix/fun/abs_test.cpp
@@ -1,7 +1,13 @@
 #include <test/unit/math/test_ad.hpp>
+#include <cmath>
+#include <complex>
+#include <vector>
 
 TEST(mixScalFun, abs) {
-  auto f = [](const auto& x) { return stan::math::abs(x); };
+  auto f = [](const auto& x) {
+    using std::abs;
+    return abs(x);
+  };
   stan::test::expect_common_nonzero_unary(f);
   stan::test::expect_value(f, 0);
   stan::test::expect_value(f, 0.0);
@@ -15,4 +21,11 @@ TEST(mixScalFun, abs) {
   stan::test::expect_ad(f, 0.68);
   stan::test::expect_ad(f, 2.0);
   stan::test::expect_ad(f, 4.0);
+
+  // not differentiable at zero
+  for (double re : std::vector<double>{-4, -2.5, -1.5, -0.3, 1.3, 2.1, 3.9}) {
+    for (double im : std::vector<double>{-4, -2.5, -1.5, -0.3, 1.3, 2.1, 3.9}) {
+      stan::test::expect_ad(f, std::complex<double>(re, im));
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/acos_test.cpp b/test/unit/math/mix/fun/acos_test.cpp
index 7d274d7d19c..350fc5ef687 100644
--- a/test/unit/math/mix/fun/acos_test.cpp
+++ b/test/unit/math/mix/fun/acos_test.cpp
@@ -1,13 +1,23 @@
 #include <test/unit/math/test_ad.hpp>
+#include <complex>
 #include <limits>
+#include <vector>
 
 TEST(mathMixMatFun, acos) {
   using stan::test::expect_unary_vectorized;
-  auto f = [](const auto& x1) { return stan::math::acos(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::acos;
+    return acos(x1);
+  };
   // can't autodiff acos through integers
   for (auto x : stan::test::internal::common_nonzero_args())
     stan::test::expect_unary_vectorized(f, x);
   expect_unary_vectorized(f, -2.2, -0.8, 0.5,
                           1 + std::numeric_limits<double>::epsilon(), 1.5, 3,
                           3.4, 4);
+  for (double re : std::vector<double>{-0.2, 0, 0.3}) {
+    for (double im : std::vector<double>{-0.3, 0, 0.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/acosh_test.cpp b/test/unit/math/mix/fun/acosh_test.cpp
index 3784354e961..3f3bd50ee4c 100644
--- a/test/unit/math/mix/fun/acosh_test.cpp
+++ b/test/unit/math/mix/fun/acosh_test.cpp
@@ -1,8 +1,18 @@
 #include <test/unit/math/test_ad.hpp>
+#include <vector>
 
 TEST(mathMixMatFun, acosh) {
-  auto f = [](const auto& x1) { return stan::math::acosh(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::acosh;
+    return acosh(x1);
+  };
   for (double x : stan::test::internal::common_args())
     stan::test::expect_unary_vectorized(x);
   stan::test::expect_unary_vectorized(f, 1.5, 3.2, 5, 10, 12.9);
+  // avoid pole at complex zero that can't be autodiffed
+  for (double re : std::vector<double>{-0.2, 0, 0.3}) {
+    for (double im : std::vector<double>{-0.3, 0.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/arg_test.cpp b/test/unit/math/mix/fun/arg_test.cpp
new file mode 100644
index 00000000000..93c62aa1f09
--- /dev/null
+++ b/test/unit/math/mix/fun/arg_test.cpp
@@ -0,0 +1,15 @@
+#include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <vector>
+
+TEST(mixScalFun, arg) {
+  auto f = [](const auto& x) { return arg(x); };
+
+  // undefined with 0 in denominator
+  stan::test::expect_ad(f, std::complex<double>(0.9, 0.8));
+  for (double re : std::vector<double>{-3.6, -0.0, 0.0, 0.5}) {
+    for (double im : std::vector<double>{-1.3, 2.3}) {
+      stan::test::expect_ad(f, std::complex<double>(re, im));
+    }
+  }
+}
diff --git a/test/unit/math/mix/fun/asin_test.cpp b/test/unit/math/mix/fun/asin_test.cpp
index badc174d3f2..6873a587510 100644
--- a/test/unit/math/mix/fun/asin_test.cpp
+++ b/test/unit/math/mix/fun/asin_test.cpp
@@ -1,10 +1,18 @@
 #include <test/unit/math/test_ad.hpp>
-#include <limits>
+#include <vector>
 
 TEST(mathMixMatFun, asin) {
-  auto f = [](const auto& x1) { return stan::math::asin(x1); };
+  auto f = [](const auto& x) {
+    using stan::math::asin;
+    return asin(x);
+  };
   // can't autodiff asin through integers
   for (auto x : stan::test::internal::common_nonzero_args())
     stan::test::expect_unary_vectorized(f, x);
   stan::test::expect_unary_vectorized(f, -2.6, -2, -0.2, 0.7, 1.3, 3.4, 5);
+  for (double re : std::vector<double>{-0.2, 0, 0.3}) {
+    for (double im : std::vector<double>{-0.3, 0, 0.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/asinh_test.cpp b/test/unit/math/mix/fun/asinh_test.cpp
index 9e4c518bdff..f218bccdac1 100644
--- a/test/unit/math/mix/fun/asinh_test.cpp
+++ b/test/unit/math/mix/fun/asinh_test.cpp
@@ -1,7 +1,19 @@
 #include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <vector>
 
-TEST(mathMixMatFun, asinh) {
-  auto f = [](const auto& x1) { return stan::math::asinh(x1); };
+TEST(mathMixFun, asinh) {
+  auto f = [](const auto& x1) {
+    using stan::math::asinh;
+    return asinh(x1);
+  };
   stan::test::expect_common_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2.6, -1.2, -0.2, 0.5, 2, -1.2);
+  stan::test::expect_ad(f, std::complex<double>{0.2});
+  // avoid pole at real zero that can't be autodiffed
+  for (double re : std::vector<double>{-0.2, 0.3}) {
+    for (double im : std::vector<double>{-0.3, 0, 0.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/atan_test.cpp b/test/unit/math/mix/fun/atan_test.cpp
index 506bcc16099..3c668654fc4 100644
--- a/test/unit/math/mix/fun/atan_test.cpp
+++ b/test/unit/math/mix/fun/atan_test.cpp
@@ -1,7 +1,18 @@
 #include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <vector>
 
 TEST(mathMixMatFun, atan) {
-  auto f = [](const auto& x1) { return stan::math::atan(x1); };
+  auto f = [](const auto& x) {
+    using stan::math::atan;
+    return atan(x);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2.6, -2, -0.2, 0.5, 1, 1.3, 1.5, 3);
+  // avoid 0 imaginary component where autodiff doesn't work
+  for (double re : std::vector<double>{-0.2, 0, 0.3}) {
+    for (double im : std::vector<double>{-0.3, 0.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/atanh_test.cpp b/test/unit/math/mix/fun/atanh_test.cpp
index 831ce1c697e..fbf9dcf6da4 100644
--- a/test/unit/math/mix/fun/atanh_test.cpp
+++ b/test/unit/math/mix/fun/atanh_test.cpp
@@ -1,7 +1,16 @@
 #include <test/unit/math/test_ad.hpp>
+#include <vector>
 
 TEST(mathMixMatFun, atanh) {
-  auto f = [](const auto& x1) { return stan::math::atanh(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::atanh;
+    return atanh(x1);
+  };
   stan::test::expect_common_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -0.9, 0.5);
+  for (double re : std::vector<double>{-0.2, 0, 0.3}) {
+    for (double im : std::vector<double>{-0.3, 0, 0.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/conj_test.cpp b/test/unit/math/mix/fun/conj_test.cpp
new file mode 100644
index 00000000000..877ced049c2
--- /dev/null
+++ b/test/unit/math/mix/fun/conj_test.cpp
@@ -0,0 +1,8 @@
+#include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <vector>
+
+TEST(mixScalFun, conj) {
+  auto f = [](const auto& x) { return conj(x); };
+  stan::test::expect_complex_common(f);
+}
diff --git a/test/unit/math/mix/fun/copysign_test.cpp b/test/unit/math/mix/fun/copysign_test.cpp
index 67cd1685393..e7b111c994f 100644
--- a/test/unit/math/mix/fun/copysign_test.cpp
+++ b/test/unit/math/mix/fun/copysign_test.cpp
@@ -1,24 +1,72 @@
 #include <stan/math/mix.hpp>
 #include <gtest/gtest.h>
+#include <cmath>
+#include <complex>
 #include <limits>
 #include <vector>
 
+template <typename U, typename V>
+void expect_eq_signbit(const U& u, const V& v) {
+  using stan::math::signbit;
+  using std::signbit;
+  EXPECT_EQ(signbit(u), signbit(v));
+}
+
 template <typename T>
 void expect_copysign() {
   using stan::math::copysign;
-  using stan::math::signbit;
+  using stan::math::copysign_non_zero;
+  using stan::math::is_nan;
+  using stan::math::value_of_rec;
   using std::copysign;
   using std::numeric_limits;
-  using std::signbit;
-  T inf = numeric_limits<double>::infinity();
-  T nan = std::numeric_limits<double>::quiet_NaN();
-  std::vector<T> ys{inf, -inf, -1, 0, 1};  // inf, -inf, -1.0, 0.0, 1.0 };
 
-  for (const T& x : ys) {
-    for (const T& y : ys) {
-      EXPECT_EQ(signbit(y), signbit(copysign(T(x), y)));
-      EXPECT_EQ(signbit(y), signbit(copysign(x, T(y))));
-      EXPECT_EQ(signbit(y), signbit(copysign(T(x), T(y))));
+  double nan = numeric_limits<double>::quiet_NaN();
+  double inf = numeric_limits<double>::infinity();
+  std::vector<double> vs{-inf, -1, -0.0, 0.0, 1, inf, nan};
+
+  // real
+  for (double x : vs) {
+    for (double y : vs) {
+      expect_eq_signbit(y, copysign(T(x), y));
+      expect_eq_signbit(y, copysign(x, T(y)));
+      expect_eq_signbit(y, copysign(T(x), T(y)));
+    }
+  }
+
+  // complex
+  for (double re1 : vs) {
+    for (double im1 : vs) {
+      auto x_d = std::complex<double>(re1, im1);
+      auto x_t = std::complex<T>(re1, im1);
+      for (double re2 : vs) {
+        for (double im2 : vs) {
+          auto y_d = std::complex<double>(re2, im2);
+          auto y_t = std::complex<T>(re2, im2);
+          if (re1 == 0) {
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_d, y_d).real()));
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_d, y_t).real()));
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_t, y_d).real()));
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_t, y_t).real()));
+          } else {
+            expect_eq_signbit(re2, copysign(x_d, y_d).real());
+            expect_eq_signbit(re2, copysign(x_d, y_t).real());
+            expect_eq_signbit(re2, copysign(x_t, y_d).real());
+            expect_eq_signbit(re2, copysign(x_t, y_t).real());
+          }
+          if (im1 == 0) {
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_d, y_d).imag()));
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_d, y_t).imag()));
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_t, y_d).imag()));
+            EXPECT_FLOAT_EQ(0.0, value_of_rec(copysign(x_t, y_t).imag()));
+          } else {
+            expect_eq_signbit(im2, copysign(x_d, y_d).imag());
+            expect_eq_signbit(im2, copysign(x_d, y_t).imag());
+            expect_eq_signbit(im2, copysign(x_t, y_d).imag());
+            expect_eq_signbit(im2, copysign(x_t, y_t).imag());
+          }
+        }
+      }
     }
   }
 }
diff --git a/test/unit/math/mix/fun/cos_test.cpp b/test/unit/math/mix/fun/cos_test.cpp
index ccd3a79f5c7..bda64a2f8b6 100644
--- a/test/unit/math/mix/fun/cos_test.cpp
+++ b/test/unit/math/mix/fun/cos_test.cpp
@@ -1,8 +1,12 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, cos) {
-  auto f = [](const auto& x1) { return stan::math::cos(x1); };
+  auto f = [](const auto& x) {
+    using stan::math::cos;
+    return cos(x);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2.6, -2, -0.2, -0.5, 0, 1.5, 3, 5,
                                       5.3);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/cosh_test.cpp b/test/unit/math/mix/fun/cosh_test.cpp
index 93922de9b05..4f9c90c5a51 100644
--- a/test/unit/math/mix/fun/cosh_test.cpp
+++ b/test/unit/math/mix/fun/cosh_test.cpp
@@ -1,8 +1,12 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, cosh) {
-  auto f = [](const auto& x1) { return stan::math::cosh(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::cosh;
+    return cosh(x1);
+  };
   stan::test::expect_common_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2.6, -2, -1.2, -0.2, 0.5, 1, 1.3,
                                       1.5);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/exp_test.cpp b/test/unit/math/mix/fun/exp_test.cpp
index ce8ea90d6c1..1a75ed232e8 100644
--- a/test/unit/math/mix/fun/exp_test.cpp
+++ b/test/unit/math/mix/fun/exp_test.cpp
@@ -1,8 +1,12 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, exp) {
-  auto f = [](const auto& x1) { return stan::math::exp(x1); };
+  auto f = [](const auto& x) {
+    using stan::math::exp;
+    return exp(x);
+  };
   stan::test::expect_common_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -15.2, -10, -0.5, 0.5, 1, 1.0, 1.3, 5,
                                       10);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/i_times_test.cpp b/test/unit/math/mix/fun/i_times_test.cpp
new file mode 100644
index 00000000000..5234906194d
--- /dev/null
+++ b/test/unit/math/mix/fun/i_times_test.cpp
@@ -0,0 +1,16 @@
+#include <test/unit/math/test_ad.hpp>
+
+TEST(matPrimFun, iTimes) {
+  auto f = [](const auto& x) { return stan::math::i_times(x); };
+  stan::test::expect_complex_common(f);
+}
+
+TEST(matPrimFun, negITimes) {
+  auto f = [](const auto& x) { return stan::math::neg_i_times(x); };
+  stan::test::expect_complex_common(f);
+}
+
+TEST(matPrimFun, complexNegate) {
+  auto f = [](const auto& x) { return stan::math::complex_negate(x); };
+  stan::test::expect_complex_common(f);
+}
diff --git a/test/unit/math/mix/fun/imag_test.cpp b/test/unit/math/mix/fun/imag_test.cpp
new file mode 100644
index 00000000000..195cb6f5c65
--- /dev/null
+++ b/test/unit/math/mix/fun/imag_test.cpp
@@ -0,0 +1,6 @@
+#include <test/unit/math/test_ad.hpp>
+
+TEST(mathMixMatFun, imag) {
+  auto f = [](const auto& z) { return imag(z); };
+  stan::test::expect_complex_common(f);
+}
diff --git a/test/unit/math/mix/fun/log10_test.cpp b/test/unit/math/mix/fun/log10_test.cpp
index 8b66bddfde4..b572db0dcd5 100644
--- a/test/unit/math/mix/fun/log10_test.cpp
+++ b/test/unit/math/mix/fun/log10_test.cpp
@@ -1,8 +1,27 @@
 #include <test/unit/math/test_ad.hpp>
+#include <vector>
 
 TEST(mathMixMatFun, log10) {
-  auto f = [](const auto& x1) { return stan::math::log10(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::log10;
+    return log10(x1);
+  };
   stan::test::expect_common_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -0.2, 1e-3, 1, 1.3, 3, 3.7, 10, 10.2,
                                       1e6);
+
+  // non-zero real and imaginary components
+  for (auto re : std::vector<double>{-2.7, 1, 2.3}) {
+    for (auto im : std::vector<double>{-1.5, 1.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
+  // zero tests which lead to finite vals
+  stan::test::expect_ad(f, std::complex<double>{0, 2.1});
+  stan::test::expect_ad(f, std::complex<double>{0, -2.1});
+  stan::test::expect_ad(f, std::complex<double>{2.1, 0});
+  stan::test::expect_ad(f, std::complex<double>{-0.0, 2.1});
+  stan::test::expect_ad(f, std::complex<double>{-0.0, -2.1});
+  stan::test::expect_ad(f, std::complex<double>{2.1, -0.0});
+  // (negative real and zero imaginary illegal)
 }
diff --git a/test/unit/math/mix/fun/log_test.cpp b/test/unit/math/mix/fun/log_test.cpp
index 6dcd7314132..21ff2c8fa03 100644
--- a/test/unit/math/mix/fun/log_test.cpp
+++ b/test/unit/math/mix/fun/log_test.cpp
@@ -1,8 +1,27 @@
 #include <test/unit/math/test_ad.hpp>
+#include <vector>
 
 TEST(mathMixMatFun, log) {
-  auto f = [](const auto& x1) { return stan::math::log(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::log;
+    return log(x1);
+  };
   stan::test::expect_common_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -0.2, 1e-3, 1, 1.3, 3, 3.7, 10, 10.2,
                                       1e6);
+
+  // non-zero real and imaginary components
+  for (auto re : std::vector<double>{-2.7, 1, 2.3}) {
+    for (auto im : std::vector<double>{-1.5, 1.2}) {
+      stan::test::expect_ad(f, std::complex<double>{re, im});
+    }
+  }
+  // zero tests which lead to finite vals
+  stan::test::expect_ad(f, std::complex<double>{0, 2.1});
+  stan::test::expect_ad(f, std::complex<double>{0, -2.1});
+  stan::test::expect_ad(f, std::complex<double>{2.1, 0});
+  stan::test::expect_ad(f, std::complex<double>{-0.0, 2.1});
+  stan::test::expect_ad(f, std::complex<double>{-0.0, -2.1});
+  stan::test::expect_ad(f, std::complex<double>{2.1, -0.0});
+  // (negative real and zero imaginary illegal)
 }
diff --git a/test/unit/math/mix/fun/logb_test.cpp b/test/unit/math/mix/fun/logb_test.cpp
new file mode 100644
index 00000000000..840ee6f3d47
--- /dev/null
+++ b/test/unit/math/mix/fun/logb_test.cpp
@@ -0,0 +1,30 @@
+#include <test/unit/math/test_ad.hpp>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+template <typename T>
+void expect_logb(double x) {
+  using std::logb;
+  T xt(x);
+  stan::test::expect_near_rel("logb", logb(x), logb(xt));
+}
+
+void expect_all_logb(double x) {
+  using stan::math::fvar;
+  using stan::math::var;
+  expect_logb<double>(x);
+  expect_logb<var>(x);
+  expect_logb<fvar<double>>(x);
+  expect_logb<fvar<fvar<double>>>(x);
+  expect_logb<fvar<var>>(x);
+  expect_logb<fvar<fvar<var>>>(x);
+}
+
+TEST(mathMixMatFun, logb) {
+  double inf = std::numeric_limits<double>::infinity();
+  double nan = std::numeric_limits<double>::quiet_NaN();
+  for (double x : std::vector<double>{-inf, -2.3, 0, 1, 2, 3.9, inf, nan}) {
+    expect_all_logb(x);
+  }
+}
diff --git a/test/unit/math/mix/fun/norm_test.cpp b/test/unit/math/mix/fun/norm_test.cpp
new file mode 100644
index 00000000000..9fec91afb6d
--- /dev/null
+++ b/test/unit/math/mix/fun/norm_test.cpp
@@ -0,0 +1,8 @@
+#include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <vector>
+
+TEST(mixScalFun, norm) {
+  auto f = [](const auto& x) { return norm(x); };
+  stan::test::expect_complex_common(f);
+}
diff --git a/test/unit/math/mix/fun/polar_test.cpp b/test/unit/math/mix/fun/polar_test.cpp
new file mode 100644
index 00000000000..f1a1a8805b9
--- /dev/null
+++ b/test/unit/math/mix/fun/polar_test.cpp
@@ -0,0 +1,50 @@
+#include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <iostream>
+#include <vector>
+
+TEST(mathMixFun, polar) {
+  auto f = [](const auto& r, const auto& theta) {
+    using stan::math::polar;
+    auto y = polar(r, theta);
+    return y;
+  };
+
+  stan::test::expect_ad(f, 0.2, 0.5);
+  stan::test::expect_ad(f, 1, 0.5);
+  stan::test::expect_ad(f, 0.2, 1);
+  stan::test::expect_ad(f, 1, 1);
+
+  auto f1 = [](const auto& r, const auto& theta) {
+    using stan::math::polar;
+    auto y = polar(r, theta);
+    return y.real();
+  };
+
+  auto f2 = [](const auto& r, const auto& theta) {
+    using stan::math::polar;
+    auto y = polar(r, theta);
+    return y.imag();
+  };
+  stan::test::expect_ad(f1, 1, 1);
+  stan::test::expect_ad(f2, 1, 1);
+
+  std::vector<double> rs{-2.3, -0.3, 0.4, 1.2};
+  std::vector<double> thetas{-0.2, 0, 0.1};
+
+  for (double r : rs) {
+    for (double theta : thetas) {
+      stan::test::expect_ad(f1, r, theta);
+      stan::test::expect_ad(f2, r, theta);
+    }
+  }
+
+  for (double r : rs) {
+    stan::test::expect_ad(f1, r, 1);
+    stan::test::expect_ad(f2, r, 1);
+  }
+  for (double theta : thetas) {
+    stan::test::expect_ad(f1, 1, theta);
+    stan::test::expect_ad(f2, 1, theta);
+  }
+}
diff --git a/test/unit/math/mix/fun/pow_test.cpp b/test/unit/math/mix/fun/pow_test.cpp
index 199001770a6..4a48fa235b7 100644
--- a/test/unit/math/mix/fun/pow_test.cpp
+++ b/test/unit/math/mix/fun/pow_test.cpp
@@ -1,22 +1,31 @@
 #include <test/unit/math/test_ad.hpp>
 #include <cmath>
 #include <limits>
+#include <vector>
 
 template <typename T>
 void expect_arith_instantiate() {
   using stan::math::pow;
   using std::pow;
-  auto a = pow(T(1.0), 1);
-  auto b = pow(T(1.0), 1.0);
-  auto c = pow(1, T(1.0));
-  auto d = pow(1.0, T(1.0));
+  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));
 }
 
 // this one's been tricky to instantiate, so test all instances
 TEST(mathMixScalFun, powInstantiations) {
   using stan::math::fvar;
   using stan::math::var;
-  expect_arith_instantiate<int>();
   expect_arith_instantiate<double>();
   expect_arith_instantiate<var>();
   expect_arith_instantiate<fvar<double>>();
@@ -46,3 +55,123 @@ TEST(mathMixScalFun, pow) {
   stan::test::expect_ad(f, nan, 1.0);
   stan::test::expect_ad(f, nan, nan);
 }
+TEST(mathMixFun, complexPow) {
+  auto f = [](const auto& x1, const auto& x2) {
+    using std::pow;
+    using stan::math::pow;
+    return pow(x1, x2);
+  };
+  stan::test::ad_tolerances tols;
+  tols.hessian_hessian_ = 5e-3;
+  tols.hessian_fvar_hessian_ = 5e-3;
+  // complex, complex
+  for (auto re1 : std::vector<double>{-1.8, 3.4}) {
+    for (auto im1 : std::vector<double>{-1.8, 3.4}) {
+      for (auto re2 : std::vector<double>{-2.7, 1, 2.3}) {
+        for (auto im2 : std::vector<double>{-1.5, 1.2}) {
+          stan::test::expect_ad(tols, f, std::complex<double>{re1, im1},
+                                std::complex<double>{re2, im2});
+        }
+      }
+    }
+  }
+  // if real, first arg must be positive
+  for (auto re1 : std::vector<double>{3.4}) {
+    for (auto re2 : std::vector<double>{-2.7, 1, 2.3}) {
+      for (auto im2 : std::vector<double>{-1.5, 1.2}) {
+        stan::test::expect_ad(tols, f, re1, std::complex<double>{re2, im2});
+      }
+    }
+  }
+  for (auto re1 : std::vector<double>{-1.8, 3.4}) {
+    for (auto im1 : std::vector<double>{-1.8, 3.4}) {
+      for (auto re2 : std::vector<double>{-2.7, 1, 2.3}) {
+        stan::test::expect_ad(tols, f, std::complex<double>{re1, im1}, re2);
+      }
+    }
+  }
+}
+
+TEST(mathMixFun, powIntAmbiguityTest) {
+  using stan::math::pow;  // included to check ambiguities
+  using stan::math::var;
+  using std::complex;
+  using std::pow;
+  int i = 2;
+  double d = 2.5;
+  var v = 2.5;
+  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 e = a1 + a2 + a3 + a4 + a5 + b1 + b2 + b3 + b4 + b5 + c1 + c2 + c3 + c4
+           + c5 + d1 + d2 + d3 + d4 + d5;
+}
+
+TEST(mathMixFun, powIntAmbiguityTestFvar) {
+  using stan::math::fvar;
+  using stan::math::pow;  // included to check ambiguities
+  using std::complex;
+  using std::pow;
+  int i = 2;
+  double d = 2.5;
+  fvar<double> v = 2.5;
+  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 e = a1 + a2 + a3 + a4 + a5 + b1 + b2 + b3 + b4 + b5 + c1 + c2 + c3 + c4
+           + c5 + d1 + d2 + d3 + d4 + d5;
+}
diff --git a/test/unit/math/mix/fun/proj_test.cpp b/test/unit/math/mix/fun/proj_test.cpp
new file mode 100644
index 00000000000..ea0006cbc7c
--- /dev/null
+++ b/test/unit/math/mix/fun/proj_test.cpp
@@ -0,0 +1,8 @@
+#include <test/unit/math/test_ad.hpp>
+#include <complex>
+#include <vector>
+
+TEST(mixScalFun, proj) {
+  auto f = [](const auto& x) { return proj(x); };
+  stan::test::expect_complex_common(f);
+}
diff --git a/test/unit/math/mix/fun/real_test.cpp b/test/unit/math/mix/fun/real_test.cpp
new file mode 100644
index 00000000000..94be5b39e52
--- /dev/null
+++ b/test/unit/math/mix/fun/real_test.cpp
@@ -0,0 +1,6 @@
+#include <test/unit/math/test_ad.hpp>
+
+TEST(mathMixMatFun, real) {
+  auto f = [](const auto& z) { return real(z); };
+  stan::test::expect_complex_common(f);
+}
diff --git a/test/unit/math/mix/fun/scalbn_test.cpp b/test/unit/math/mix/fun/scalbn_test.cpp
new file mode 100644
index 00000000000..d68618fc38e
--- /dev/null
+++ b/test/unit/math/mix/fun/scalbn_test.cpp
@@ -0,0 +1,32 @@
+#include <test/unit/math/test_ad.hpp>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+template <typename T>
+void expect_scalbn(double x) {
+  using std::scalbn;
+  T xt(x);
+  for (int n = 2; n < 5; ++n) {
+    stan::test::expect_near_rel("scalbn", scalbn(x, n), scalbn(xt, n));
+  }
+}
+
+void expect_all_scalbn(double x) {
+  using stan::math::fvar;
+  using stan::math::var;
+  expect_scalbn<double>(x);
+  expect_scalbn<var>(x);
+  expect_scalbn<fvar<double>>(x);
+  expect_scalbn<fvar<fvar<double>>>(x);
+  expect_scalbn<fvar<var>>(x);
+  expect_scalbn<fvar<fvar<var>>>(x);
+}
+
+TEST(mathMixMatFun, scalbn) {
+  double inf = std::numeric_limits<double>::infinity();
+  double nan = std::numeric_limits<double>::quiet_NaN();
+  for (double x : std::vector<double>{2.3}) {
+    expect_all_scalbn(x);
+  }
+}
diff --git a/test/unit/math/mix/fun/sin_test.cpp b/test/unit/math/mix/fun/sin_test.cpp
index 88ce5db4630..896351d76ee 100644
--- a/test/unit/math/mix/fun/sin_test.cpp
+++ b/test/unit/math/mix/fun/sin_test.cpp
@@ -1,7 +1,11 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, sin) {
-  auto f = [](const auto& x1) { return stan::math::sin(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::sin;
+    return sin(x1);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2.6, -2, -0.2, 3, 5, 5.3);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/sinh_test.cpp b/test/unit/math/mix/fun/sinh_test.cpp
index b1282ee79b9..de285a71651 100644
--- a/test/unit/math/mix/fun/sinh_test.cpp
+++ b/test/unit/math/mix/fun/sinh_test.cpp
@@ -1,8 +1,12 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, sinh) {
-  auto f = [](const auto& x1) { return stan::math::sinh(x1); };
+  auto f = [](const auto& x) {
+    using stan::math::sinh;
+    return sinh(x);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2, -1.2, -0.5, -0.2, 0.5, 1.3, 1.5,
                                       3);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/sqrt_test.cpp b/test/unit/math/mix/fun/sqrt_test.cpp
index 251042259f8..16d2c5c0be5 100644
--- a/test/unit/math/mix/fun/sqrt_test.cpp
+++ b/test/unit/math/mix/fun/sqrt_test.cpp
@@ -1,7 +1,19 @@
 #include <test/unit/math/test_ad.hpp>
+#include <vector>
 
 TEST(mathMixMatFun, sqrt) {
-  auto f = [](const auto& x1) { return stan::math::sqrt(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::sqrt;
+    return sqrt(x1);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -6, -5.2, 1.3, 7, 10.7, 36, 1e6);
+
+  // undefined with 0 in denominator
+  stan::test::expect_ad(f, std::complex<double>(0.9, 0.8));
+  for (double im : std::vector<double>{-1.3, 2.3}) {
+    for (double re : std::vector<double>{-3.6, -0.0, 0.0, 0.5}) {
+      stan::test::expect_ad(f, std::complex<double>(re, im));
+    }
+  }
 }
diff --git a/test/unit/math/mix/fun/tan_test.cpp b/test/unit/math/mix/fun/tan_test.cpp
index 36c994d35e1..6570cd3b48b 100644
--- a/test/unit/math/mix/fun/tan_test.cpp
+++ b/test/unit/math/mix/fun/tan_test.cpp
@@ -1,7 +1,11 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, tan) {
-  auto f = [](const auto& x1) { return stan::math::tan(x1); };
+  auto f = [](const auto& x) {
+    using stan::math::tan;
+    return tan(x);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2, -0.5, 0.5, 1.5, 3, 4.4);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/mix/fun/tanh_test.cpp b/test/unit/math/mix/fun/tanh_test.cpp
index b98a28298dc..78bca8cef01 100644
--- a/test/unit/math/mix/fun/tanh_test.cpp
+++ b/test/unit/math/mix/fun/tanh_test.cpp
@@ -1,7 +1,11 @@
 #include <test/unit/math/test_ad.hpp>
 
 TEST(mathMixMatFun, tanh) {
-  auto f = [](const auto& x1) { return stan::math::tanh(x1); };
+  auto f = [](const auto& x1) {
+    using stan::math::tanh;
+    return tanh(x1);
+  };
   stan::test::expect_common_nonzero_unary_vectorized(f);
   stan::test::expect_unary_vectorized(f, -2.6, -2, -1.2, -0.5, 0.5, 1.5);
+  stan::test::expect_complex_common(f);
 }
diff --git a/test/unit/math/rev/core/var_test.cpp b/test/unit/math/rev/core/var_test.cpp
index 72b455af121..54d66f75db0 100644
--- a/test/unit/math/rev/core/var_test.cpp
+++ b/test/unit/math/rev/core/var_test.cpp
@@ -69,16 +69,6 @@ TEST_F(AgradRev, ctorOverloads) {
   // ptrdiff_t
   EXPECT_FLOAT_EQ(37, var(static_cast<ptrdiff_t>(37)).val());
   EXPECT_FLOAT_EQ(0, var(static_cast<ptrdiff_t>(0)).val());
-
-  // complex but with zero imaginary part
-  EXPECT_FLOAT_EQ(37, var(std::complex<double>(37, 0)).val());
-  EXPECT_FLOAT_EQ(37, var(std::complex<float>(37, 0)).val());
-  EXPECT_FLOAT_EQ(37, var(std::complex<long double>(37, 0)).val());
-
-  // complex but with non-zero imaginary part
-  EXPECT_THROW(var(std::complex<double>(37, 10)), std::invalid_argument);
-  EXPECT_THROW(var(std::complex<float>(37, 10)), std::invalid_argument);
-  EXPECT_THROW(var(std::complex<long double>(37, 10)), std::invalid_argument);
 }
 
 TEST_F(AgradRev, a_eq_x) {
diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp
index 031ed8b449e..d9c9e5dd94c 100644
--- a/test/unit/math/test_ad.hpp
+++ b/test/unit/math/test_ad.hpp
@@ -1721,6 +1721,100 @@ auto ldlt_factor(const Eigen::Matrix<T, -1, -1>& x) {
   return ldlt_x;
 }
 
+std::vector<double> common_complex_parts() {
+  return {-4, -2.5, -1.5, -0.3, -0.0, 0.0, 1.3, 2.1, 3.9};
+}
+
+std::vector<std::complex<double>> common_complex() {
+  std::vector<std::complex<double>> zs;
+  for (double re : common_complex_parts())
+    for (double im : common_complex_parts())
+      zs.emplace_back(re, im);
+  return zs;
+}
+
+template <typename F>
+void expect_complex_common(const F& f) {
+  auto zs = common_complex();
+  for (auto z : zs) {
+    expect_ad(f, z);
+  }
+}
+
+template <typename F>
+void expect_complex_common_binary(const F& f) {
+  auto xs = common_complex_parts();
+  auto zs = common_complex();
+  // complex, complex
+  for (auto z1 : zs) {
+    for (auto z2 : zs) {
+      expect_ad(f, z1, z2);
+    }
+  }
+  // complex, real
+  for (auto z1 : zs) {
+    for (auto x2 : xs) {
+      expect_ad(f, z1, x2);
+    }
+  }
+  // real, complex
+  for (auto x1 : xs) {
+    for (auto z2 : zs) {
+      expect_ad(f, x1, z2);
+    }
+  }
+}
+
+template <typename T, typename F>
+void expect_complex_compare(const F& f, const std::complex<double>& z1,
+                            const std::complex<double>& z2) {
+  using c_t = std::complex<T>;
+  c_t cz1{z1};
+  c_t cz2{z2};
+  T z1r{z1.real()};
+  T z2r{z2.real()};
+
+  EXPECT_EQ(f(z1, z2), f(cz1, cz2));
+  EXPECT_EQ(f(z1, z2), f(cz1, z2));
+  EXPECT_EQ(f(z1, z2), f(z1, cz2));
+
+  EXPECT_EQ(f(z1.real(), z2), f(z1r, cz2));
+  EXPECT_EQ(f(z1.real(), z2), f(z1r, z2));
+
+  EXPECT_EQ(f(z1, z2.real()), f(cz1, z2r));
+  EXPECT_EQ(f(z1, z2.real()), f(z1, z2r));
+}
+
+template <typename F>
+void expect_complex_comparison(const F& f, const std::complex<double>& z1,
+                               const std::complex<double>& z2) {
+  using stan::math::fvar;
+  using stan::math::var;
+  using std::complex;
+  expect_complex_compare<double>(f, z1, z2);              // PASS
+  expect_complex_compare<var>(f, z1, z2);                 // FAIL
+  expect_complex_compare<fvar<double>>(f, z1, z2);        // PASS
+  expect_complex_compare<fvar<fvar<double>>>(f, z1, z2);  // PASS
+  expect_complex_compare<fvar<var>>(f, z1, z2);           // PASS
+  expect_complex_compare<fvar<fvar<var>>>(f, z1, z2);     // PASS
+}
+
+/**
+ * Test the specified comparison operation provides results matching
+ * those for the double version for all the common complex numbers.
+ *
+ * @tparam F type of function to test
+ * @param f function to test
+ */
+template <typename F>
+void expect_complex_common_comparison(const F& f) {
+  for (auto z1 : common_complex()) {
+    for (auto z2 : common_complex()) {
+      expect_complex_comparison(f, z1, z2);
+    }
+  }
+}
+
 }  // namespace test
 }  // namespace stan
 #endif