diff --git a/make/libraries b/make/libraries
index 76ef9bfe823..0af959afb78 100644
--- a/make/libraries
+++ b/make/libraries
@@ -78,8 +78,8 @@ $(BOOST_LIB)/mpi.so: $(BOOST)/user-config.jam
 $(BOOST_LIB)/libboost_serialization.so: $(BOOST_LIB)/mpi.so
 
 $(BOOST_LIB)/libboost_serialization.dylib: $(BOOST_LIB)/mpi.so
-	install_name_tool -add_rpath "$(BOOST_LIB_ABS)" "$(BOOST_LIB)/libboost_serialization.dylib"
-	install_name_tool -id @rpath/libboost_serialization.dylib "$(BOOST_LIB)/libboost_serialization.dylib"
+	install_name_tool -add_rpath "$(BOOST_LIB_ABS)" $(BOOST_LIB)/libboost_serialization.dylib
+	install_name_tool -id @rpath/libboost_serialization.dylib $(BOOST_LIB)/libboost_serialization.dylib
 
 $(BOOST_LIB)/libboost_mpi.so: $(BOOST_LIB)/mpi.so
 
diff --git a/stan/math/fwd/core/fvar.hpp b/stan/math/fwd/core/fvar.hpp
index 3d4defb41a9..a2689a511c1 100644
--- a/stan/math/fwd/core/fvar.hpp
+++ b/stan/math/fwd/core/fvar.hpp
@@ -4,8 +4,8 @@
 #include <stan/math/prim/scal/meta/likely.hpp>
 #include <stan/math/prim/scal/fun/is_nan.hpp>
 #include <stan/math/fwd/scal/meta/ad_promotable.hpp>
-#include <boost/utility/enable_if.hpp>
 #include <ostream>
+#include <type_traits>
 
 namespace stan {
 namespace math {
@@ -97,13 +97,10 @@ struct fvar {
    * @tparam V type of value (must be assignable to the value and
    *   tangent type T)
    * @param[in] v value
-   * @param[in] dummy value given by default with enable-if
-   *   metaprogramming
    */
-  template <typename V>
-  fvar(const V& v,
-       typename boost::enable_if_c<ad_promotable<V, T>::value>::type* dummy = 0)
-      : val_(v), d_(0.0) {
+  template <typename V,
+            typename std::enable_if_t<ad_promotable<V, T>::value>* = nullptr>
+  fvar(const V& v) : val_(v), d_(0.0) {  // NOLINT(runtime/explicit)
     if (unlikely(is_nan(v)))
       d_ = v;
   }
diff --git a/stan/math/fwd/core/operator_division.hpp b/stan/math/fwd/core/operator_division.hpp
index 6ed5bdb168e..9c5fe2308be 100644
--- a/stan/math/fwd/core/operator_division.hpp
+++ b/stan/math/fwd/core/operator_division.hpp
@@ -46,6 +46,7 @@ inline fvar<T> operator/(double x1, const fvar<T>& x2) {
   // TODO(carpenter): store x1 / x2.val_ and reuse
   return fvar<T>(x1 / x2.val_, -x1 * x2.d_ / (x2.val_ * x2.val_));
 }
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/cplx.hpp b/stan/math/fwd/cplx.hpp
new file mode 100644
index 00000000000..69ba16bb0d7
--- /dev/null
+++ b/stan/math/fwd/cplx.hpp
@@ -0,0 +1,7 @@
+#ifndef STAN_MATH_FWD_CPLX_HPP
+#define STAN_MATH_FWD_CPLX_HPP
+
+#include <stan/math/fwd/cplx/fvar.hpp>
+#include <stan/math/fwd/cplx/operator_division.hpp>
+
+#endif
diff --git a/stan/math/fwd/cplx/fvar.hpp b/stan/math/fwd/cplx/fvar.hpp
new file mode 100644
index 00000000000..3f1a1dec339
--- /dev/null
+++ b/stan/math/fwd/cplx/fvar.hpp
@@ -0,0 +1,21 @@
+#ifndef STAN_MATH_FWD_CPLX_FVAR_HPP
+#define STAN_MATH_FWD_CPLX_FVAR_HPP
+
+#include <stan/math/fwd/scal/meta/ad_promotable.hpp>
+#include <stan/math/fwd/scal/meta/is_fvar.hpp>
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/prim/cplx/complex.hpp>
+
+namespace std {
+
+/**
+ *  std::complex<fvar> inherits from stan::math::complex. Details
+ *  are provided there.
+ */
+template <class T>
+struct complex<stan::math::fvar<T>> : stan::math::complex<stan::math::fvar<T>> {
+  using stan::math::complex<stan::math::fvar<T>>::complex;
+};
+
+}  // namespace std
+#endif
diff --git a/stan/math/fwd/cplx/operator_division.hpp b/stan/math/fwd/cplx/operator_division.hpp
new file mode 100644
index 00000000000..3044d63ff4b
--- /dev/null
+++ b/stan/math/fwd/cplx/operator_division.hpp
@@ -0,0 +1,29 @@
+#ifndef STAN_MATH_FWD_CPLX_OPERATOR_DIVISION_HPP
+#define STAN_MATH_FWD_CPLX_OPERATOR_DIVISION_HPP
+
+#include <stan/math/fwd/cplx/fvar.hpp>
+#include <stan/math/fwd/scal/fun/pow.hpp>
+#include <stan/math/prim/cplx/operator_division.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the result of dividing the first argument by the second.
+ * This overload exists because libc++ uses logb and scalbn for
+ * complex division, which aren't defined for fvars.
+ *
+ * @tparam T type of complex fvar value and tangent
+ * @param t first argument
+ * @param u second argument
+ * @return first argument divided by second argument
+ */
+template <class T>
+inline std::complex<fvar<T>> operator/(std::complex<fvar<T>> const& t,
+                                       std::complex<fvar<T>> const& u) {
+  return stan::math::operator_division(t, u);  // no recursion
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp b/stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp
index 3d96fe52758..6350704555e 100644
--- a/stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp
+++ b/stan/math/fwd/mat/vectorize/apply_scalar_unary.hpp
@@ -18,7 +18,7 @@ namespace math {
  * autodiff variable.
  */
 template <typename F, typename T>
-struct apply_scalar_unary<F, fvar<T> > {
+struct apply_scalar_unary<F, fvar<T>> {
   /**
    * Function return type, which is same as the argument type for
    * the function, <code>fvar&lt;T&gt;</code>.
diff --git a/stan/math/fwd/scal.hpp b/stan/math/fwd/scal.hpp
index 3cbd4cc2d0d..6fc07280971 100644
--- a/stan/math/fwd/scal.hpp
+++ b/stan/math/fwd/scal.hpp
@@ -1,6 +1,8 @@
 #ifndef STAN_MATH_FWD_SCAL_HPP
 #define STAN_MATH_FWD_SCAL_HPP
 
+#include <stan/math/fwd/cplx.hpp>
+
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/scal/meta/ad_promotable.hpp>
 #include <stan/math/fwd/scal/meta/is_fvar.hpp>
@@ -79,6 +81,7 @@
 #include <stan/math/fwd/scal/fun/Phi_approx.hpp>
 #include <stan/math/fwd/scal/fun/pow.hpp>
 #include <stan/math/fwd/scal/fun/primitive_value.hpp>
+#include <stan/math/fwd/scal/fun/proj.hpp>
 #include <stan/math/fwd/scal/fun/rising_factorial.hpp>
 #include <stan/math/fwd/scal/fun/round.hpp>
 #include <stan/math/fwd/scal/fun/sin.hpp>
diff --git a/stan/math/fwd/scal/fun/abs.hpp b/stan/math/fwd/scal/fun/abs.hpp
index 14eb7dc123f..a35b3f62f11 100644
--- a/stan/math/fwd/scal/fun/abs.hpp
+++ b/stan/math/fwd/scal/fun/abs.hpp
@@ -1,7 +1,9 @@
 #ifndef STAN_MATH_FWD_SCAL_FUN_ABS_HPP
 #define STAN_MATH_FWD_SCAL_FUN_ABS_HPP
 
+#include <stan/math/fwd/cplx.hpp>
 #include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/scal/fun/hypot.hpp>
 #include <stan/math/fwd/scal/fun/value_of.hpp>
 #include <stan/math/prim/scal/fun/constants.hpp>
 #include <stan/math/prim/scal/fun/abs.hpp>
@@ -25,4 +27,12 @@ inline fvar<T> abs(const fvar<T>& x) {
 
 }  // namespace math
 }  // namespace stan
+
+namespace std {
+// constrained complex overload to forward zeroing fvar to std::abs
+template <class T>
+inline stan::math::fvar<T> abs(std::complex<stan::math::fvar<T>> const& t) {
+  return abs(std::complex<stan::math::zeroing<stan::math::fvar<T>>>(t));
+}
+}  // namespace std
 #endif
diff --git a/stan/math/fwd/scal/fun/is_inf.hpp b/stan/math/fwd/scal/fun/is_inf.hpp
index db3d1f54c0e..c34f9da8e76 100644
--- a/stan/math/fwd/scal/fun/is_inf.hpp
+++ b/stan/math/fwd/scal/fun/is_inf.hpp
@@ -20,6 +20,12 @@ inline int is_inf(const fvar<T>& x) {
   return is_inf(x.val());
 }
 
+// forwarding for ADL
+template <class T>
+inline auto isinf(const fvar<T>& a) {
+  return is_inf(a);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/scal/fun/is_nan.hpp b/stan/math/fwd/scal/fun/is_nan.hpp
index 29492b41f1a..984942724eb 100644
--- a/stan/math/fwd/scal/fun/is_nan.hpp
+++ b/stan/math/fwd/scal/fun/is_nan.hpp
@@ -20,6 +20,12 @@ inline int is_nan(const fvar<T>& x) {
   return is_nan(x.val());
 }
 
+// forwarding for ADL
+template <class T>
+inline auto isnan(const fvar<T>& a) {
+  return is_nan(a);
+}
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/fwd/scal/fun/proj.hpp b/stan/math/fwd/scal/fun/proj.hpp
new file mode 100644
index 00000000000..e3096e6fa9b
--- /dev/null
+++ b/stan/math/fwd/scal/fun/proj.hpp
@@ -0,0 +1,28 @@
+#ifndef STAN_MATH_FWD_SCAL_FUN_PROJ_HPP
+#define STAN_MATH_FWD_SCAL_FUN_PROJ_HPP
+
+#include <stan/math/fwd/cplx.hpp>
+#include <stan/math/fwd/scal/fun/is_inf.hpp>
+#include <stan/math/prim/scal/fun/proj.hpp>
+
+namespace std {
+
+/**
+ * Return a projection of a complex number onto
+ * the Riemann sphere.
+ *
+ * This overloads an erroneous definition in
+ * libstdc++ for general types.
+ *
+ * @tparam T auto diff variable type
+ * @param[in] t Variable input.
+ * @return projection onto Riemann sphere
+ */
+template <class T>
+inline std::complex<stan::math::fvar<T>> proj(
+    std::complex<stan::math::fvar<T>> const& t) {
+  return stan::math::proj(t);
+}
+
+}  // namespace std
+#endif
diff --git a/stan/math/prim/cplx.hpp b/stan/math/prim/cplx.hpp
new file mode 100644
index 00000000000..df55310abb1
--- /dev/null
+++ b/stan/math/prim/cplx.hpp
@@ -0,0 +1,11 @@
+#ifndef STAN_MATH_PRIM_CPLX_HPP
+#define STAN_MATH_PRIM_CPLX_HPP
+
+#include <stan/math/prim/cplx/complex.hpp>
+#include <stan/math/prim/cplx/operator_addition.hpp>
+#include <stan/math/prim/cplx/operator_division.hpp>
+#include <stan/math/prim/cplx/operator_multiplication.hpp>
+#include <stan/math/prim/cplx/operator_subtraction.hpp>
+#include <stan/math/prim/cplx/zeroing.hpp>
+
+#endif
diff --git a/stan/math/prim/cplx/complex.hpp b/stan/math/prim/cplx/complex.hpp
new file mode 100644
index 00000000000..edd95624822
--- /dev/null
+++ b/stan/math/prim/cplx/complex.hpp
@@ -0,0 +1,57 @@
+#ifndef STAN_MATH_PRIM_CPLX_COMPLEX_HPP
+#define STAN_MATH_PRIM_CPLX_COMPLEX_HPP
+
+#include <stan/math/prim/cplx/zeroing.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ *  Complex class that forwards the interface of std::complex, brining
+ *  stan::math's namespace into ADL for functions called on data types
+ *  that inherit from this class, as well as allowing template
+ *  specializations of std::complex to indirectly work with different
+ *  underlying data types. For example, std::complex<var> specializes
+ *  the std::complex<T> template to inherit from
+ *  stan::math::complex<T>, which in turn inherits from
+ *  std::complex<zeroing<T>>. This causes stan::math to be pulled into
+ *  ADL on operations involving std::complex<var>, to handle, e.g.
+ *  division operations that the base std::complex doesn't have
+ *  defined usably for var, and, since stan::math::zeroing<var> zero
+ *  initializes itself, will also correctly work with the remaining
+ *  algorithms in std::complex<T> that require the expression T() to
+ *  be 0.
+ */
+template <class T>
+struct complex : std::complex<zeroing<T>> {
+  using std::complex<zeroing<T>>::complex;
+  // NOLINTNEXTLINE(runtime/explicit)
+  complex(const std::complex<zeroing<T>>& t)
+      : std::complex<zeroing<T>>(t) {}  ///< allow promotion
+  /**
+   * Without this converting ctor, copy initialization of
+   * std::complex<zeroing<(f)var>>> from an (f)var fails, because the
+   * default implementation in the std::complex template requires the
+   * argument type to match the template type. This is the same reason
+   * why std::complex<double> can't be multiplied by an int.
+   *
+   * tparam R type of real argument. Won't fire on non-arithmetic types.
+   * tparam I type of imaginary argument.
+   * param[in] real, the pure real component of the complex number.
+   * param[in] imag, the pure imaginary component of the complex number
+   */
+  template <class R = double, class I = R,
+            std::enable_if_t<is_arith_like<R>::value>* = nullptr>
+  // NOLINTNEXTLINE(runtime/explicit)
+  complex(R const& real = 0, I const& imag = 0)
+      : std::complex<zeroing<T>>(zeroing<T>(real), zeroing<T>(imag)) {}
+  // downcast
+  operator std::complex<T>() const {
+    return *static_cast<std::complex<T> const*>(this);
+  }
+};
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/cplx/operator_addition.hpp b/stan/math/prim/cplx/operator_addition.hpp
new file mode 100644
index 00000000000..f9457ea61b4
--- /dev/null
+++ b/stan/math/prim/cplx/operator_addition.hpp
@@ -0,0 +1,50 @@
+#ifndef STAN_MATH_PRIM_CPLX_OPERATOR_ADDITION_HPP
+#define STAN_MATH_PRIM_CPLX_OPERATOR_ADDITION_HPP
+
+#include <stan/math/prim/scal/fun/copysign.hpp>
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
+#include <stan/math/prim/scal/fun/isfinite.hpp>
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the sum of the specified arguments
+ *
+ * @tparam T type of complex first argument
+ * @tparam U type of second argument
+ * @param t complex first argument
+ * @param u second argument
+ * @return sum
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && (is_complex<U>::value
+                               || is_arith_like<U>::value)>* = nullptr>
+inline auto operator+(std::complex<T> const& t, U const& u) {
+  return complex_promote<U>(t) += u;
+}
+
+/**
+ * Return the sum of the specified arguments
+ *
+ * @tparam T type of first argument
+ * @tparam U type of complex second argument
+ * @param t first argument
+ * @param u complex second argument
+ * @return sum
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && is_arith_like<T>::value>* = nullptr>
+inline auto operator+(T const& t, std::complex<U> const& u) {
+  return complex_promote<U>(t) += u;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/cplx/operator_division.hpp b/stan/math/prim/cplx/operator_division.hpp
new file mode 100644
index 00000000000..2433a2baa1d
--- /dev/null
+++ b/stan/math/prim/cplx/operator_division.hpp
@@ -0,0 +1,89 @@
+#ifndef STAN_MATH_PRIM_CPLX_OPERATOR_DIVISION_HPP
+#define STAN_MATH_PRIM_CPLX_OPERATOR_DIVISION_HPP
+
+#include <stan/math/prim/scal/fun/copysign.hpp>
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
+#include <stan/math/prim/scal/fun/isfinite.hpp>
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <complex>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the result of dividing the first argument by the second.
+ *
+ * @tparam T type of complex first argument
+ * @tparam U type of second argument
+ * @param t complex first argument
+ * @param u second argument
+ * @return first argument divided by second argument
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && (is_complex<U>::value
+                               || is_arith_like<U>::value)>* = nullptr>
+inline auto operator/(std::complex<T> const& t, U const& u) {
+  return complex_promote<U>(t) /= u;
+}
+
+/**
+ * Return the result of dividing the first argument by the second.
+ *
+ * @tparam T type of first argument
+ * @tparam U type of complex second argument
+ * @param t first argument
+ * @param u complex second argument
+ * @return first argument divided by second argument
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && is_arith_like<T>::value>* = nullptr>
+inline auto operator/(T const& t, std::complex<U> const& u) {
+  return complex_promote<U>(t) /= u;
+}
+
+/**
+ * Return the result of dividing the first argument by the second.
+ * This function exists because libc++ uses logb and scalbn for
+ * complex division, which aren't defined for fvars. Four tightly
+ * scoped high priority operator/ overloads for fvar and var
+ * forward to this function for implementation.
+ *
+ * @tparam T type of complex value
+ * @param t first argument
+ * @param u second argument
+ * @return first argument divided by second argument
+ */
+template <class T>
+inline auto operator_division(std::complex<T> const& t,
+                              std::complex<T> const& u) {
+  T const n(norm(u));
+  T const r((t.real() * u.real() + t.imag() * u.imag()) / n);
+  T const i((t.imag() * u.real() - t.real() * u.imag()) / n);
+  return std::complex<rm_zeroing_t<T>>(r, i);
+}
+
+/**
+ * Return the result of dividing the first argument by the second.
+ *
+ * This overload exists because libc++ uses logb and scalbn for
+ * complex division, which aren't defined for fvars.
+ *
+ * @tparam T type of complex fvar value and tangent
+ * @param t first argument
+ * @param u second argument
+ * @return first argument divided by second argument
+ */
+template <class T>
+inline std::complex<rm_zeroing_t<T>> operator/(
+    std::complex<zeroing<T>> const& t, std::complex<zeroing<T>> const& u) {
+  return stan::math::operator_division(t, u);  // no recursion
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/cplx/operator_multiplication.hpp b/stan/math/prim/cplx/operator_multiplication.hpp
new file mode 100644
index 00000000000..c34fff99fdb
--- /dev/null
+++ b/stan/math/prim/cplx/operator_multiplication.hpp
@@ -0,0 +1,50 @@
+#ifndef STAN_MATH_PRIM_CPLX_OPERATOR_MULTIPLICATION_HPP
+#define STAN_MATH_PRIM_CPLX_OPERATOR_MULTIPLICATION_HPP
+
+#include <stan/math/prim/scal/fun/copysign.hpp>
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
+#include <stan/math/prim/scal/fun/isfinite.hpp>
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the product of the specified arguments
+ *
+ * @tparam T type of complex first argument
+ * @tparam U type of second argument
+ * @param t complex first argument
+ * @param u second argument
+ * @return product
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && (is_complex<U>::value
+                               || is_arith_like<U>::value)>* = nullptr>
+inline auto operator*(std::complex<T> const& t, U const& u) {
+  return complex_promote<U>(t) *= u;
+}
+
+/**
+ * Return the product of the specified arguments
+ *
+ * @tparam T type of first argument
+ * @tparam U type of complex second argument
+ * @param t first argument
+ * @param u complex second argument
+ * @return product
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && is_arith_like<T>::value>* = nullptr>
+inline auto operator*(T const& t, std::complex<U> const& u) {
+  return complex_promote<U>(t) *= u;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/cplx/operator_subtraction.hpp b/stan/math/prim/cplx/operator_subtraction.hpp
new file mode 100644
index 00000000000..49e255f9963
--- /dev/null
+++ b/stan/math/prim/cplx/operator_subtraction.hpp
@@ -0,0 +1,50 @@
+#ifndef STAN_MATH_PRIM_CPLX_OPERATOR_SUBTRACTION_HPP
+#define STAN_MATH_PRIM_CPLX_OPERATOR_SUBTRACTION_HPP
+
+#include <stan/math/prim/scal/fun/copysign.hpp>
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
+#include <stan/math/prim/scal/fun/isfinite.hpp>
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the difference of the specified arguments
+ *
+ * @tparam T type of complex first argument
+ * @tparam U type of second argument
+ * @param t complex first argument
+ * @param u second argument
+ * @return difference
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && (is_complex<U>::value
+                               || is_arith_like<U>::value)>* = nullptr>
+inline auto operator-(std::complex<T> const& t, U const& u) {
+  return complex_promote<U>(t) -= u;
+}
+
+/**
+ * Return the difference of the specified arguments
+ *
+ * @tparam T type of first argument
+ * @tparam U type of complex second argument
+ * @param t first argument
+ * @param u complex second argument
+ * @return difference
+ */
+template <class T, class U,
+          std::enable_if_t<is_fr_var<T, U>::value
+                           && is_arith_like<T>::value>* = nullptr>
+inline auto operator-(T const& t, std::complex<U> const& u) {
+  return complex_promote<U>(t) -= u;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/cplx/zeroing.hpp b/stan/math/prim/cplx/zeroing.hpp
new file mode 100644
index 00000000000..45859454508
--- /dev/null
+++ b/stan/math/prim/cplx/zeroing.hpp
@@ -0,0 +1,28 @@
+#ifndef STAN_MATH_PRIM_CPLX_ZEROING_HPP
+#define STAN_MATH_PRIM_CPLX_ZEROING_HPP
+
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ *  zeroing<T> is a zeroing T that std::complex<T> implicitly
+ *  uses  when inheriting from stan::math::complex. It has the
+ *  property that zeroing() results in T(0). See stan's
+ *  complex class for details. Without this, in libstdc++'s
+ *  complex class, when T() is called, vars (and fvar<vars>)
+ *  will have a null pointer, which is bad because libstdc++
+ *  immediately uses such expressions as if they were
+ *  initialized to zero..
+ */
+template <class T = double>
+struct zeroing : T {
+  template <class U = T,
+            std::enable_if_t<stan::is_arith_like<U>::value>* = nullptr>
+  zeroing(U const& t = 0) : T(t) {}  // NOLINT(runtime/explicit)
+};
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp
index 34bfdb22f1e..c06015697a3 100644
--- a/stan/math/prim/mat.hpp
+++ b/stan/math/prim/mat.hpp
@@ -12,6 +12,7 @@
 #include <stan/math/prim/arr/meta/length.hpp>
 
 #include <stan/math/prim/mat/meta/broadcast_array.hpp>
+#include <stan/math/prim/mat/meta/Eigen_ScalarBinaryOpTraits.hpp>
 #include <stan/math/prim/mat/meta/get.hpp>
 #include <stan/math/prim/mat/meta/index_type.hpp>
 #include <stan/math/prim/mat/meta/is_constant_struct.hpp>
diff --git a/stan/math/prim/mat/fun/gp_dot_prod_cov.hpp b/stan/math/prim/mat/fun/gp_dot_prod_cov.hpp
index 73b03aea68e..830f6215c67 100644
--- a/stan/math/prim/mat/fun/gp_dot_prod_cov.hpp
+++ b/stan/math/prim/mat/fun/gp_dot_prod_cov.hpp
@@ -8,7 +8,7 @@
 #include <stan/math/prim/scal/err/check_not_nan.hpp>
 #include <stan/math/prim/scal/fun/square.hpp>
 #include <stan/math/prim/scal/meta/return_type.hpp>
-#include <stan/math/rev/mat/fun/dot_product.hpp>
+#include <stan/math/prim/mat/fun/dot_product.hpp>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/prim/mat/meta/Eigen_ScalarBinaryOpTraits.hpp b/stan/math/prim/mat/meta/Eigen_ScalarBinaryOpTraits.hpp
new file mode 100644
index 00000000000..a1583acbd23
--- /dev/null
+++ b/stan/math/prim/mat/meta/Eigen_ScalarBinaryOpTraits.hpp
@@ -0,0 +1,38 @@
+#ifndef STAN_MATH_PRIM_MAT_META_EIGEN_SCALARBINARYOPTRAITS_HPP
+#define STAN_MATH_PRIM_MAT_META_EIGEN_SCALARBINARYOPTRAITS_HPP
+
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/prim/scal/meta/rm_complex.hpp>
+#include <stan/math/prim/scal/meta/rm_zeroing.hpp>
+#include <Eigen/Dense>
+#include <complex>
+#include <type_traits>
+
+namespace Eigen {
+/// @cond DO_NOT_DOCUMENT
+// Eigen scalar op traits specialization for complex variables
+template <class T1, class T2, template <class, class> class OP>
+struct ScalarBinaryOpTraits<
+    T1,
+    std::enable_if_t<
+        // !VectorBlock !is_eigen
+        (stan::is_complex<T1>::value || stan::is_arith_like<T1>::value)
+            && (stan::is_complex<T2>::value || stan::is_arith_like<T2>::value)
+            && !std::is_same<T1, T2>::value &&  // avoid Eigen's template
+            ((stan::is_complex<T1>::value &&    // next boolean avoids Eigen
+              !std::is_same<stan::rm_complex_t<T1>,
+                            stan::rm_zeroing_t<T2>>::value)
+             || (stan::is_complex<T2>::value &&  // next bool avoids Eigen
+                 !std::is_same<stan::rm_zeroing_t<T1>,
+                               stan::rm_complex_t<T2>>::value)),
+        T2>,
+    OP<T1, T2>> {
+  typedef std::complex<typename stan::return_type<stan::rm_complex_t<T1>,
+                                                  stan::rm_complex_t<T2>>::type>
+      ReturnType;
+};
+/// @endcond
+}  // namespace Eigen
+#endif
diff --git a/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp b/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp
index 020a7595602..a4b84e55724 100644
--- a/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp
+++ b/stan/math/prim/mat/vectorize/apply_scalar_unary.hpp
@@ -1,7 +1,9 @@
 #ifndef STAN_MATH_PRIM_MAT_VECTORIZE_APPLY_SCALAR_UNARY_HPP
 #define STAN_MATH_PRIM_MAT_VECTORIZE_APPLY_SCALAR_UNARY_HPP
 
+#include <stan/math/prim/scal/meta/rm_zeroing.hpp>
 #include <Eigen/Dense>
+#include <complex>
 #include <vector>
 
 namespace stan {
@@ -33,6 +35,17 @@ namespace math {
  */
 template <typename F, typename T>
 struct apply_scalar_unary {
+  static_assert(
+      std::is_base_of<
+          Eigen::EigenBase<typename Eigen::internal::remove_all<T>::type>,
+          typename Eigen::internal::remove_all<T>::type>::value,
+      "If this fires, it is likely that apply_scalar_unary has no "
+      "specialization "
+      "for T *OR* that the function F isn't defined for T and has thus tried "
+      "the "
+      "matrix overload. In any case, this unconstrained template struct "
+      "only works for Eigen types.");
+
   /**
    * Type of underlying scalar for the matrix type T.
    */
@@ -114,6 +127,84 @@ struct apply_scalar_unary<F, int> {
   static inline return_t apply(int x) { return F::fun(static_cast<double>(x)); }
 };
 
+/**
+ * Template specialization for vectorized functions applying to
+ * complex arguments.
+ *
+ * @tparam F Type of function defining static apply function.
+ * @tparam T Underlying type of the complex type
+ */
+template <typename F, class T>
+struct apply_scalar_unary<F, std::complex<T>> {
+  /**
+   * The return type, complex.
+   */
+  typedef std::complex<rm_zeroing_t<T>> return_t;
+
+  /**
+   * Apply the function specified by F to the specified argument.
+   * This is defined through a direct application of
+   * <code>F::fun()</code>, which must be defined for complex
+   * arguments.
+   *
+   * @param x Argument scalar.
+   * @return Result of applying F to the scalar.
+   */
+  static inline return_t apply(return_t x) { return F::fun(x); }
+};
+
+template <class>
+struct complex;  // forward declare stan's complex
+
+/**
+ * Template specialization for vectorized functions applying to
+ * stan's complex arguments.
+ *
+ * @tparam F Type of function defining static apply function.
+ * @tparam T Underlying type of the complex type
+ */
+template <typename F, class T>
+struct apply_scalar_unary<F, stan::math::complex<T>> {
+  /**
+   * The return type, complex.
+   */
+  typedef std::complex<rm_zeroing_t<T>> return_t;
+
+  /**
+   * Apply the function specified by F to the specified argument.
+   * This is defined through a direct application of
+   * <code>F::fun()</code>, which must be defined for
+   * std::complex arguments.
+   *
+   * @param x Argument scalar. (implictly upcast through ctor)
+   * @return Result of applying F to the scalar.
+   */
+  static inline return_t apply(return_t x) { return F::fun(x); }
+};
+
+/**
+ * Template specialization for vectorized functions applying
+ * to zeroing variable arguments.
+ *
+ * @tparam F Class defining a static apply function
+ * @tparam T wrapped type of the zeroing variable
+ */
+template <typename F, typename T>
+struct apply_scalar_unary<F, zeroing<T>> {
+  /**
+   * Function return type, <code>T</code>.
+   */
+  typedef T return_t;
+
+  /**
+   * Apply the function specified by F to the specified argument.
+   *
+   * @param x Argument variable
+   * @return Function applied to the variable
+   */
+  static inline return_t apply(const T& x) { return F::fun(x); }
+};
+
 /**
  * Template specialization for vectorized functions applying to
  * standard vector containers.  The lowest-level scalar type of
@@ -124,7 +215,7 @@ struct apply_scalar_unary<F, int> {
  * @tparam T Type of element contained in standard vector.
  */
 template <typename F, typename T>
-struct apply_scalar_unary<F, std::vector<T> > {
+struct apply_scalar_unary<F, std::vector<T>> {
   /**
    * Return type, which is calculated recursively as a standard
    * vector of the return type of the contained type T.
diff --git a/stan/math/prim/scal.hpp b/stan/math/prim/scal.hpp
index ed09053123a..a80deae3627 100644
--- a/stan/math/prim/scal.hpp
+++ b/stan/math/prim/scal.hpp
@@ -3,6 +3,8 @@
 
 #include <stan/math/version.hpp>
 
+#include <stan/math/prim/cplx.hpp>
+
 #include <stan/math/prim/scal/meta/ad_promotable.hpp>
 #include <stan/math/prim/scal/meta/child_type.hpp>
 #include <stan/math/prim/scal/meta/contains_fvar.hpp>
@@ -12,8 +14,11 @@
 #include <stan/math/prim/scal/meta/get.hpp>
 #include <stan/math/prim/scal/meta/include_summand.hpp>
 #include <stan/math/prim/scal/meta/index_type.hpp>
+#include <stan/math/prim/scal/meta/is_arith_like.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
 #include <stan/math/prim/scal/meta/is_constant.hpp>
 #include <stan/math/prim/scal/meta/is_constant_struct.hpp>
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
 #include <stan/math/prim/scal/meta/is_fvar.hpp>
 #include <stan/math/prim/scal/meta/is_var.hpp>
 #include <stan/math/prim/scal/meta/is_var_or_arithmetic.hpp>
@@ -29,6 +34,8 @@
 #include <stan/math/prim/scal/meta/partials_return_type.hpp>
 #include <stan/math/prim/scal/meta/partials_type.hpp>
 #include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/prim/scal/meta/rm_complex.hpp>
+#include <stan/math/prim/scal/meta/rm_zeroing.hpp>
 #include <stan/math/prim/scal/meta/scalar_seq_view.hpp>
 #include <stan/math/prim/scal/meta/scalar_type.hpp>
 #include <stan/math/prim/scal/meta/scalar_type_pre.hpp>
@@ -72,7 +79,9 @@
 #include <stan/math/prim/scal/fun/cbrt.hpp>
 #include <stan/math/prim/scal/fun/choose.hpp>
 #include <stan/math/prim/scal/fun/common_type.hpp>
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
 #include <stan/math/prim/scal/fun/constants.hpp>
+#include <stan/math/prim/scal/fun/copysign.hpp>
 #include <stan/math/prim/scal/fun/corr_constrain.hpp>
 #include <stan/math/prim/scal/fun/corr_free.hpp>
 #include <stan/math/prim/scal/fun/digamma.hpp>
@@ -114,6 +123,7 @@
 #include <stan/math/prim/scal/fun/is_inf.hpp>
 #include <stan/math/prim/scal/fun/is_nan.hpp>
 #include <stan/math/prim/scal/fun/is_uninitialized.hpp>
+#include <stan/math/prim/scal/fun/isfinite.hpp>
 #include <stan/math/prim/scal/fun/lb_constrain.hpp>
 #include <stan/math/prim/scal/fun/lb_free.hpp>
 #include <stan/math/prim/scal/fun/lbeta.hpp>
@@ -154,17 +164,21 @@
 #include <stan/math/prim/scal/fun/owens_t.hpp>
 #include <stan/math/prim/scal/fun/Phi.hpp>
 #include <stan/math/prim/scal/fun/Phi_approx.hpp>
+#include <stan/math/prim/scal/fun/polar.hpp>
 #include <stan/math/prim/scal/fun/positive_constrain.hpp>
 #include <stan/math/prim/scal/fun/positive_free.hpp>
+#include <stan/math/prim/scal/fun/pow.hpp>
 #include <stan/math/prim/scal/fun/primitive_value.hpp>
 #include <stan/math/prim/scal/fun/prob_constrain.hpp>
 #include <stan/math/prim/scal/fun/prob_free.hpp>
+#include <stan/math/prim/scal/fun/proj.hpp>
 #include <stan/math/prim/scal/fun/promote_elements.hpp>
 #include <stan/math/prim/scal/fun/promote_scalar.hpp>
 #include <stan/math/prim/scal/fun/promote_scalar_type.hpp>
 #include <stan/math/prim/scal/fun/rising_factorial.hpp>
 #include <stan/math/prim/scal/fun/round.hpp>
 #include <stan/math/prim/scal/fun/sign.hpp>
+#include <stan/math/prim/scal/fun/signbit.hpp>
 #include <stan/math/prim/scal/fun/size_zero.hpp>
 #include <stan/math/prim/scal/fun/square.hpp>
 #include <stan/math/prim/scal/fun/squared_distance.hpp>
@@ -174,6 +188,7 @@
 #include <stan/math/prim/scal/fun/trunc.hpp>
 #include <stan/math/prim/scal/fun/ub_constrain.hpp>
 #include <stan/math/prim/scal/fun/ub_free.hpp>
+#include <stan/math/prim/scal/fun/val.hpp>
 #include <stan/math/prim/scal/fun/value_of.hpp>
 #include <stan/math/prim/scal/fun/value_of_rec.hpp>
 
diff --git a/stan/math/prim/scal/fun/complex_promote.hpp b/stan/math/prim/scal/fun/complex_promote.hpp
new file mode 100644
index 00000000000..2a42f4e3797
--- /dev/null
+++ b/stan/math/prim/scal/fun/complex_promote.hpp
@@ -0,0 +1,30 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_COMPLEX_PROMOTE_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_COMPLEX_PROMOTE_HPP
+
+#include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/prim/scal/meta/rm_complex.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+/**
+ * Complex auto-diff (AD) promotion for arguments
+ * to complex functions. Handles cases where types
+ * may or may not be wrapped in complex.
+ *
+ * @tparam N non-deduced input type
+ * @tparam D deduced input type
+ * @param d argument
+ * @return value of promoted type
+ */
+template <class N, class D,
+          class P
+          = typename return_type<rm_complex_t<N>, rm_complex_t<D>>::type>
+inline auto complex_promote(D const& d) {
+  return std::complex<P>(d);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/fun/copysign.hpp b/stan/math/prim/scal/fun/copysign.hpp
new file mode 100644
index 00000000000..6f6c916038a
--- /dev/null
+++ b/stan/math/prim/scal/fun/copysign.hpp
@@ -0,0 +1,34 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_COPYSIGN_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_COPYSIGN_HPP
+
+#include <stan/math/prim/scal/fun/signbit.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the magnitude of t with the sign of u
+ *
+ * Either T or U must be a stan type for ADL
+ *
+ * Needed for libc++'s implementation of
+ * complex multiplication on stan types
+ *
+ * @tparam T type of object
+ * @tparam U type of object
+ * @param t magnitude reference
+ * @param u sign reference
+ * @return magnitude of t with the sign of u
+ */
+template <class T, class U = T>
+inline auto copysign(T const& t, U const& u) {
+  using std::fabs;
+  using std::signbit;
+  auto sign(signbit(u) ? -1 : 1);
+  return fabs(t) * sign;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/fun/grad_reg_lower_inc_gamma.hpp b/stan/math/prim/scal/fun/grad_reg_lower_inc_gamma.hpp
index c989761f905..c2096470f35 100644
--- a/stan/math/prim/scal/fun/grad_reg_lower_inc_gamma.hpp
+++ b/stan/math/prim/scal/fun/grad_reg_lower_inc_gamma.hpp
@@ -108,6 +108,7 @@ typename return_type<T1, T2>::type grad_reg_lower_inc_gamma(
   using std::exp;
   using std::log;
   using std::pow;
+  using std::sqrt;
   typedef typename return_type<T1, T2>::type TP;
 
   if (is_nan(a) || is_nan(z))
diff --git a/stan/math/prim/scal/fun/isfinite.hpp b/stan/math/prim/scal/fun/isfinite.hpp
new file mode 100644
index 00000000000..504cb80da4b
--- /dev/null
+++ b/stan/math/prim/scal/fun/isfinite.hpp
@@ -0,0 +1,36 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_ISFINITE_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_ISFINITE_HPP
+
+#include <stan/math/prim/scal/fun/val.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Checks if the given number has finite value.
+ *
+ * Return <code>true</code> if the specified variable's
+ * value is finite.
+ *
+ * Termination happens when a non-stan type is
+ * returned from val, preventing ADL from finding
+ * this declaration of isfinite again.
+ *
+ * Needed for Eigen's fullPivLu() method on
+ * non-Hermitian  matrix types. Called on stan
+ * types from Eigen::complex::isfinite_impl via ADL.
+ *
+ * @tparam T type of stan object
+ * @param t Variable to test.
+ * @return <code>true</code> if variable is finite.
+ */
+template <class T>
+inline auto isfinite(T const& t) {
+  using std::isfinite;
+  return isfinite(val(t));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/fun/polar.hpp b/stan/math/prim/scal/fun/polar.hpp
new file mode 100644
index 00000000000..295243bc99e
--- /dev/null
+++ b/stan/math/prim/scal/fun/polar.hpp
@@ -0,0 +1,34 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_POLAR_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_POLAR_HPP
+
+#include <stan/math/prim/scal/meta/return_type.hpp>
+#include <stan/math/prim/scal/meta/rm_zeroing.hpp>
+#include <complex>
+#include <cmath>
+#include <type_traits>
+
+namespace std {
+/**
+ * Returns a complex number with magnitude r
+ * and phase theta
+ *
+ * This is in the std namespace because the
+ * relevant call we are catching is fully
+ * qualified with std, making ADL useless.
+ *
+ * @tparam R magnitude type
+ * @tparam TH phase type
+ * @param r magnitude
+ * @param theta angle, radians
+ * @return complex with magnitude r & phase theta
+ */
+template <class R, class TH,
+          class RT = typename stan::return_type<stan::rm_zeroing_t<R>,
+                                                stan::rm_zeroing_t<TH>>::type,
+          std::enable_if_t<!std::is_same<R, TH>::value>* = nullptr>
+inline std::complex<RT> polar(R const& r, TH const& theta) {
+  return polar(RT(r), RT(theta));
+}
+
+}  // namespace std
+#endif
diff --git a/stan/math/prim/scal/fun/pow.hpp b/stan/math/prim/scal/fun/pow.hpp
new file mode 100644
index 00000000000..187e260ad84
--- /dev/null
+++ b/stan/math/prim/scal/fun/pow.hpp
@@ -0,0 +1,38 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_POW_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_POW_HPP
+
+#include <stan/math/prim/scal/fun/complex_promote.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <stan/math/prim/scal/meta/rm_complex.hpp>
+#include <cmath>
+#include <complex>
+#include <type_traits>
+
+namespace stan {
+namespace math {
+
+template <class T, class U
+          // in libc++ complex doesn't fire for this case
+          // so avoding T = U causes failure
+          // libstd++ may require avoiding that case
+          // ,
+          //        std::enable_if_t<!std::is_same<rm_complex_t<T>,
+          //                                       rm_complex_t<U>>::value>* =
+          //                                       nullptr
+          >
+inline auto pow(std::complex<T> const& t, U const& u) {
+  return pow(complex_promote<U>(t), complex_promote<T>(u));
+}
+
+template <class T, class U,
+          std::enable_if_t<
+              // avoid declaration above in case both t & u are complex
+              // !std::is_same<rm_complex_t<T>, rm_complex_t<U>>::value &&
+              !is_complex<T>::value>* = nullptr>
+inline auto pow(T const& t, std::complex<U> const& u) {
+  return pow(complex_promote<U>(t), complex_promote<T>(u));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/fun/proj.hpp b/stan/math/prim/scal/fun/proj.hpp
new file mode 100644
index 00000000000..aff76757189
--- /dev/null
+++ b/stan/math/prim/scal/fun/proj.hpp
@@ -0,0 +1,33 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_PROJ_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_PROJ_HPP
+
+#include <stan/math/prim/scal/fun/copysign.hpp>
+#include <stan/math/prim/scal/fun/is_inf.hpp>
+#include <complex>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return a projection of a complex number onto
+ * the Riemann sphere.
+ *
+ * This overloads an erroneous definition in
+ * libstdc++ for general types.
+ *
+ * @tparam T complex type
+ * @param[in] t complex variable input
+ * @return projection onto Riemann sphere
+ */
+
+template <class T>
+inline std::complex<T> proj(std::complex<T> const& t) {
+  if (is_inf(t.real()) || is_inf(t.imag()))
+    return std::complex<T>(INFINITY, copysign(T(0), t.imag()));
+  return t;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/fun/signbit.hpp b/stan/math/prim/scal/fun/signbit.hpp
new file mode 100644
index 00000000000..c8f3faf21a4
--- /dev/null
+++ b/stan/math/prim/scal/fun/signbit.hpp
@@ -0,0 +1,28 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_SIGNBIT_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_SIGNBIT_HPP
+
+#include <stan/math/prim/scal/fun/val.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Determines if the argument is negative
+ *
+ * Needed for libc++'s implementation of
+ * trig on complex values
+ *
+ * @tparam T type of stan AD object
+ * @param t stan AD object
+ * @return true if t is negative
+ */
+template <class T>
+inline auto signbit(T const& t) {
+  using std::signbit;
+  return signbit(val(t));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/fun/val.hpp b/stan/math/prim/scal/fun/val.hpp
new file mode 100644
index 00000000000..5d024a74fe7
--- /dev/null
+++ b/stan/math/prim/scal/fun/val.hpp
@@ -0,0 +1,83 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUN_VAL_HPP
+#define STAN_MATH_PRIM_SCAL_FUN_VAL_HPP
+
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <complex>
+#include <type_traits>
+
+namespace stan {
+namespace math {
+
+// helper descent function arithmetic case
+template <class T, std::enable_if_t<std::is_arithmetic<T>::value>* = nullptr>
+inline auto val_helper(T const& t) {
+  return t;
+}
+
+// helper descent function for stan object
+template <class T, std::enable_if_t<!std::is_arithmetic<T>::value>* = nullptr>
+inline auto val_helper(T const& t) {
+  return t.val();
+}
+
+/**
+ * Returns the underlying value of the
+ * stan type, zero or multiple levels deep
+ * if required by template parameter S.
+ *
+ * This is the termination condition
+ * of the function template.
+ *
+ * @tparam T type argument
+ * @tparam S return type sought
+ * @param t argument
+ * @param s optional return type arg
+ * @return underlying value
+ */
+template <class T, class S = decltype(val_helper(T())),
+          std::enable_if_t<std::is_same<T, S>::value>* = nullptr>
+inline S val(T const& t, S const& = S()) {
+  return t;  // terminate descent regardless of type T
+}
+
+/**
+ * Returns the underlying value of the
+ * (f)var type, zero or multiple levels deep
+ * if required by template parameter S.
+ *
+ * @tparam T type argument
+ * @tparam S return type sought
+ * @param t argument.
+ * @param s optional return type arg
+ * @return value of variable.
+ */
+template <class T, class S = decltype(val_helper(T())),
+          std::enable_if_t<!std::is_same<T, S>::value>* = nullptr>
+inline S val(T const& t, S const& s = S()) {
+  // fundamental types require
+  // val to be fully qualified
+  return stan::math::val(val_helper(t), s);
+}
+
+/**
+ * Returns the underlying value of the
+ * stan type, multiple levels deep if
+ * required by template parameter S.
+ *
+ * @tparam T complex type argument
+ * @tparam S complex return type sought
+ * @param t argument
+ * @param s optional return type arg
+ * @return underlying value
+ */
+template <class T, class S = decltype(val_helper(T()))>
+inline std::complex<S> val(std::complex<T> const& t, S const& s = S()) {
+  return std::complex<S>(
+      // fundamental types require
+      // val to be fully qualified
+      stan::math::val(t.real(), s), stan::math::val(t.imag(), s));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/meta/is_arith_like.hpp b/stan/math/prim/scal/meta/is_arith_like.hpp
new file mode 100644
index 00000000000..af4db9da2d7
--- /dev/null
+++ b/stan/math/prim/scal/meta/is_arith_like.hpp
@@ -0,0 +1,37 @@
+#ifndef STAN_MATH_PRIM_SCAL_META_IS_ARITH_LIKE_HPP
+#define STAN_MATH_PRIM_SCAL_META_IS_ARITH_LIKE_HPP
+
+#include <stan/math/prim/scal/meta/is_fr_var.hpp>
+#include <stan/math/prim/scal/meta/is_complex.hpp>
+#include <type_traits>
+
+namespace stan {
+
+/**
+ * std::true_type if T is fvar, var, or arithmetic
+ * std::false_type if T is complex or other
+ * restriction on complex matches std::is_arithmetic
+ */
+template <class T>
+struct is_arith_like_helper
+    : std::integral_constant<bool, !is_complex<T>::value
+                                       && (is_fr_var<T>::value
+                                           || std::is_arithmetic<T>::value)> {};
+
+/**
+ * std::true_type if any parameter is fvar, var,
+ * or arithmetic. Restriction on complex matches
+ * std::is_arithmetic.
+ */
+template <class T, class U = void, class V = void, class W = void,
+          class X = void, class Y = void>
+struct is_arith_like
+    : std::integral_constant<bool, is_arith_like_helper<T>::value
+                                       || is_arith_like_helper<U>::value
+                                       || is_arith_like_helper<V>::value
+                                       || is_arith_like_helper<W>::value
+                                       || is_arith_like_helper<X>::value
+                                       || is_arith_like_helper<Y>::value> {};
+
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/meta/is_complex.hpp b/stan/math/prim/scal/meta/is_complex.hpp
new file mode 100644
index 00000000000..fe7f6831294
--- /dev/null
+++ b/stan/math/prim/scal/meta/is_complex.hpp
@@ -0,0 +1,25 @@
+#ifndef STAN_MATH_PRIM_SCAL_META_IS_COMPLEX_HPP
+#define STAN_MATH_PRIM_SCAL_META_IS_COMPLEX_HPP
+
+#include <complex>
+
+namespace stan {
+namespace math {
+
+template <class T>
+struct complex;  // forward declare stan's complex
+
+}  // namespace math
+
+/// trait to see if the template parameter is complex
+template <class>
+struct is_complex : std::false_type {};
+
+template <class T>
+struct is_complex<std::complex<T>> : std::true_type {};
+
+template <class T>
+struct is_complex<stan::math::complex<T>> : std::true_type {};
+
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/meta/is_fr_var.hpp b/stan/math/prim/scal/meta/is_fr_var.hpp
new file mode 100644
index 00000000000..be1db43ce46
--- /dev/null
+++ b/stan/math/prim/scal/meta/is_fr_var.hpp
@@ -0,0 +1,35 @@
+#ifndef STAN_MATH_PRIM_SCAL_META_IS_FR_VAR_HPP
+#define STAN_MATH_PRIM_SCAL_META_IS_FR_VAR_HPP
+
+#include <stan/math/prim/scal/meta/is_fvar.hpp>
+#include <stan/math/prim/scal/meta/is_var.hpp>
+#include <stan/math/prim/scal/meta/rm_complex.hpp>
+#include <type_traits>
+
+namespace stan {
+
+/**
+ * Is std::true_type if T is fvar or var.
+ * Strips outer complex types.
+ * Note rm_complex also does rm_zeroing.
+ */
+template <class T>
+struct is_fr_var_helper
+    : std::integral_constant<bool, is_fvar<rm_complex_t<T>>::value
+                                       || is_var<rm_complex_t<T>>::value> {};
+
+/**
+ * Is std::true_type if any parameter is fvar or var.
+ * Strips outer complex types..
+ */
+template <class T, class U = double, class V = double, class W = double,
+          class X = double, class Y = double>
+struct is_fr_var
+    : std::integral_constant<
+          bool, is_fr_var_helper<T>::value || is_fr_var_helper<U>::value
+                    || is_fr_var_helper<V>::value || is_fr_var_helper<W>::value
+                    || is_fr_var_helper<X>::value
+                    || is_fr_var_helper<Y>::value> {};
+
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/meta/rm_complex.hpp b/stan/math/prim/scal/meta/rm_complex.hpp
new file mode 100644
index 00000000000..fb27c562415
--- /dev/null
+++ b/stan/math/prim/scal/meta/rm_complex.hpp
@@ -0,0 +1,32 @@
+#ifndef STAN_MATH_PRIM_SCAL_META_RM_COMPLEX_HPP
+#define STAN_MATH_PRIM_SCAL_META_RM_COMPLEX_HPP
+
+#include <stan/math/prim/scal/meta/rm_zeroing.hpp>
+#include <complex>
+
+namespace stan {
+namespace math {
+
+template <class T>
+struct complex;  // forward declare stan's complex
+
+}  // namespace math
+
+/// trait to remove the complex wrapper around a type
+template <class T>
+struct rm_complex {
+  typedef rm_zeroing_t<T> type;
+};
+template <class T>
+struct rm_complex<std::complex<T>> {
+  typedef rm_zeroing_t<T> type;
+};
+template <class T>
+struct rm_complex<stan::math::complex<T>> {
+  typedef rm_zeroing_t<T> type;
+};
+template <class T>
+using rm_complex_t = typename rm_complex<T>::type;
+
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/meta/rm_zeroing.hpp b/stan/math/prim/scal/meta/rm_zeroing.hpp
new file mode 100644
index 00000000000..738e31253b7
--- /dev/null
+++ b/stan/math/prim/scal/meta/rm_zeroing.hpp
@@ -0,0 +1,34 @@
+#ifndef STAN_MATH_PRIM_SCAL_META_RM_ZEROING_HPP
+#define STAN_MATH_PRIM_SCAL_META_RM_ZEROING_HPP
+
+namespace stan {
+namespace math {
+
+template <class T>
+struct zeroing;  // forward declare
+
+}  // namespace math
+
+/**
+ * rm_zeroing replaces <code>zeroing<T></code> with
+ * T. This is the general version of the trait.
+ */
+template <class T>
+struct rm_zeroing {
+  typedef T type;
+};
+
+/**
+ * rm_zeroing replaces <code>zeroing<t></code> with
+ * T. This is the removal version of the trait.
+ */
+template <class T>
+struct rm_zeroing<math::zeroing<T>> {
+  typedef T type;
+};
+
+template <class T>
+using rm_zeroing_t = typename rm_zeroing<T>::type;
+
+}  // namespace stan
+#endif
diff --git a/stan/math/prim/scal/prob/beta_binomial_cdf.hpp b/stan/math/prim/scal/prob/beta_binomial_cdf.hpp
index 5a923a675a9..636de40a908 100644
--- a/stan/math/prim/scal/prob/beta_binomial_cdf.hpp
+++ b/stan/math/prim/scal/prob/beta_binomial_cdf.hpp
@@ -68,7 +68,6 @@ typename return_type<T_size1, T_size2>::type beta_binomial_cdf(
   scalar_seq_view<T_size2> beta_vec(beta);
   size_t size = max_size(n, N, alpha, beta);
 
-  using std::exp;
   using std::exp;
 
   operands_and_partials<T_size1, T_size2> ops_partials(alpha, beta);
diff --git a/stan/math/prim/scal/prob/beta_binomial_lccdf.hpp b/stan/math/prim/scal/prob/beta_binomial_lccdf.hpp
index 750d9c80f5c..3f156132ce3 100644
--- a/stan/math/prim/scal/prob/beta_binomial_lccdf.hpp
+++ b/stan/math/prim/scal/prob/beta_binomial_lccdf.hpp
@@ -68,7 +68,6 @@ typename return_type<T_size1, T_size2>::type beta_binomial_lccdf(
   scalar_seq_view<T_size2> beta_vec(beta);
   size_t size = max_size(n, N, alpha, beta);
 
-  using std::exp;
   using std::exp;
   using std::log;
 
diff --git a/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp b/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp
index 39af5b28cc6..9d88e6e499d 100644
--- a/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp
+++ b/stan/math/prim/scal/prob/beta_binomial_lcdf.hpp
@@ -68,7 +68,6 @@ typename return_type<T_size1, T_size2>::type beta_binomial_lcdf(
   scalar_seq_view<T_size2> beta_vec(beta);
   size_t size = max_size(n, N, alpha, beta);
 
-  using std::exp;
   using std::exp;
   using std::log;
 
diff --git a/stan/math/prim/scal/prob/beta_lccdf.hpp b/stan/math/prim/scal/prob/beta_lccdf.hpp
index 660edcad5b9..41482cda9d1 100644
--- a/stan/math/prim/scal/prob/beta_lccdf.hpp
+++ b/stan/math/prim/scal/prob/beta_lccdf.hpp
@@ -80,7 +80,6 @@ typename return_type<T_y, T_scale_succ, T_scale_fail>::type beta_lccdf(
   operands_and_partials<T_y, T_scale_succ, T_scale_fail> ops_partials(y, alpha,
                                                                       beta);
 
-  using std::exp;
   using std::exp;
   using std::log;
   using std::pow;
diff --git a/stan/math/prim/scal/prob/beta_lcdf.hpp b/stan/math/prim/scal/prob/beta_lcdf.hpp
index 0af7208eccc..b1d64dc92c0 100644
--- a/stan/math/prim/scal/prob/beta_lcdf.hpp
+++ b/stan/math/prim/scal/prob/beta_lcdf.hpp
@@ -80,7 +80,6 @@ typename return_type<T_y, T_scale_succ, T_scale_fail>::type beta_lcdf(
   operands_and_partials<T_y, T_scale_succ, T_scale_fail> ops_partials(y, alpha,
                                                                       beta);
 
-  using std::exp;
   using std::exp;
   using std::log;
   using std::pow;
diff --git a/stan/math/prim/scal/prob/binomial_cdf.hpp b/stan/math/prim/scal/prob/binomial_cdf.hpp
index 745238b7654..3cbfd06d896 100644
--- a/stan/math/prim/scal/prob/binomial_cdf.hpp
+++ b/stan/math/prim/scal/prob/binomial_cdf.hpp
@@ -69,7 +69,6 @@ typename return_type<T_prob>::type binomial_cdf(const T_n& n, const T_N& N,
   scalar_seq_view<T_prob> theta_vec(theta);
   size_t size = max_size(n, N, theta);
 
-  using std::exp;
   using std::exp;
   using std::pow;
 
diff --git a/stan/math/prim/scal/prob/binomial_lccdf.hpp b/stan/math/prim/scal/prob/binomial_lccdf.hpp
index e7c0c850687..c920e0fa709 100644
--- a/stan/math/prim/scal/prob/binomial_lccdf.hpp
+++ b/stan/math/prim/scal/prob/binomial_lccdf.hpp
@@ -69,7 +69,6 @@ typename return_type<T_prob>::type binomial_lccdf(const T_n& n, const T_N& N,
   scalar_seq_view<T_prob> theta_vec(theta);
   size_t size = max_size(n, N, theta);
 
-  using std::exp;
   using std::exp;
   using std::log;
   using std::pow;
diff --git a/stan/math/prim/scal/prob/binomial_lcdf.hpp b/stan/math/prim/scal/prob/binomial_lcdf.hpp
index 4d4552857fb..928e8e902c2 100644
--- a/stan/math/prim/scal/prob/binomial_lcdf.hpp
+++ b/stan/math/prim/scal/prob/binomial_lcdf.hpp
@@ -69,7 +69,6 @@ typename return_type<T_prob>::type binomial_lcdf(const T_n& n, const T_N& N,
   scalar_seq_view<T_prob> theta_vec(theta);
   size_t size = max_size(n, N, theta);
 
-  using std::exp;
   using std::exp;
   using std::log;
   using std::pow;
diff --git a/stan/math/prim/scal/prob/chi_square_cdf.hpp b/stan/math/prim/scal/prob/chi_square_cdf.hpp
index 2fac95668e3..2ae911c9b4a 100644
--- a/stan/math/prim/scal/prob/chi_square_cdf.hpp
+++ b/stan/math/prim/scal/prob/chi_square_cdf.hpp
@@ -72,7 +72,6 @@ typename return_type<T_y, T_dof>::type chi_square_cdf(const T_y& y,
 
   using boost::math::tgamma;
   using std::exp;
-  using std::exp;
   using std::pow;
 
   VectorBuilder<!is_constant_struct<T_dof>::value, T_partials_return, T_dof>
diff --git a/stan/math/prim/scal/prob/chi_square_lccdf.hpp b/stan/math/prim/scal/prob/chi_square_lccdf.hpp
index c4e17602e6b..31f195e97c0 100644
--- a/stan/math/prim/scal/prob/chi_square_lccdf.hpp
+++ b/stan/math/prim/scal/prob/chi_square_lccdf.hpp
@@ -72,7 +72,6 @@ typename return_type<T_y, T_dof>::type chi_square_lccdf(const T_y& y,
 
   using boost::math::tgamma;
   using std::exp;
-  using std::exp;
   using std::log;
   using std::pow;
 
diff --git a/stan/math/prim/scal/prob/chi_square_lcdf.hpp b/stan/math/prim/scal/prob/chi_square_lcdf.hpp
index e6dfd3d6cf2..6c958001dd7 100644
--- a/stan/math/prim/scal/prob/chi_square_lcdf.hpp
+++ b/stan/math/prim/scal/prob/chi_square_lcdf.hpp
@@ -72,7 +72,6 @@ typename return_type<T_y, T_dof>::type chi_square_lcdf(const T_y& y,
 
   using boost::math::tgamma;
   using std::exp;
-  using std::exp;
   using std::log;
   using std::pow;
 
diff --git a/stan/math/prim/scal/prob/double_exponential_lccdf.hpp b/stan/math/prim/scal/prob/double_exponential_lccdf.hpp
index 5aa62e7060a..c53f5aa63a4 100644
--- a/stan/math/prim/scal/prob/double_exponential_lccdf.hpp
+++ b/stan/math/prim/scal/prob/double_exponential_lccdf.hpp
@@ -55,7 +55,6 @@ typename return_type<T_y, T_loc, T_scale>::type double_exponential_lccdf(
   check_consistent_sizes(function, "Random variable", y, "Location parameter",
                          mu, "Scale Parameter", sigma);
 
-  using std::exp;
   using std::exp;
   using std::log;
 
diff --git a/stan/math/prim/scal/prob/double_exponential_lcdf.hpp b/stan/math/prim/scal/prob/double_exponential_lcdf.hpp
index bf87e42806c..57cd09d0360 100644
--- a/stan/math/prim/scal/prob/double_exponential_lcdf.hpp
+++ b/stan/math/prim/scal/prob/double_exponential_lcdf.hpp
@@ -54,7 +54,6 @@ typename return_type<T_y, T_loc, T_scale>::type double_exponential_lcdf(
   check_consistent_sizes(function, "Random variable", y, "Location parameter",
                          mu, "Scale Parameter", sigma);
 
-  using std::exp;
   using std::exp;
   using std::log;
 
diff --git a/stan/math/prim/scal/prob/double_exponential_lpdf.hpp b/stan/math/prim/scal/prob/double_exponential_lpdf.hpp
index 8e417d37809..58c681de7da 100644
--- a/stan/math/prim/scal/prob/double_exponential_lpdf.hpp
+++ b/stan/math/prim/scal/prob/double_exponential_lpdf.hpp
@@ -48,7 +48,6 @@ typename return_type<T_y, T_loc, T_scale>::type double_exponential_lpdf(
   using stan::is_constant_struct;
   using std::fabs;
   using std::log;
-  using std::log;
 
   if (size_zero(y, mu, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/exp_mod_normal_lccdf.hpp b/stan/math/prim/scal/prob/exp_mod_normal_lccdf.hpp
index f482a40ac7e..7d00149c429 100644
--- a/stan/math/prim/scal/prob/exp_mod_normal_lccdf.hpp
+++ b/stan/math/prim/scal/prob/exp_mod_normal_lccdf.hpp
@@ -49,7 +49,6 @@ exp_mod_normal_lccdf(const T_y& y, const T_loc& mu, const T_scale& sigma,
 
   using std::exp;
   using std::log;
-  using std::log;
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_loc> mu_vec(mu);
diff --git a/stan/math/prim/scal/prob/exp_mod_normal_lcdf.hpp b/stan/math/prim/scal/prob/exp_mod_normal_lcdf.hpp
index eaceef2b08b..b6fb3cc865e 100644
--- a/stan/math/prim/scal/prob/exp_mod_normal_lcdf.hpp
+++ b/stan/math/prim/scal/prob/exp_mod_normal_lcdf.hpp
@@ -49,7 +49,6 @@ exp_mod_normal_lcdf(const T_y& y, const T_loc& mu, const T_scale& sigma,
 
   using std::exp;
   using std::log;
-  using std::log;
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_loc> mu_vec(mu);
diff --git a/stan/math/prim/scal/prob/frechet_cdf.hpp b/stan/math/prim/scal/prob/frechet_cdf.hpp
index 89ac31a9e60..d0a30084512 100644
--- a/stan/math/prim/scal/prob/frechet_cdf.hpp
+++ b/stan/math/prim/scal/prob/frechet_cdf.hpp
@@ -37,6 +37,7 @@ typename return_type<T_y, T_shape, T_scale>::type frechet_cdf(
   using boost::math::tools::promote_args;
   using std::exp;
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 1.0;
diff --git a/stan/math/prim/scal/prob/frechet_lccdf.hpp b/stan/math/prim/scal/prob/frechet_lccdf.hpp
index d32d8f2b77b..a96a7a1b72d 100644
--- a/stan/math/prim/scal/prob/frechet_lccdf.hpp
+++ b/stan/math/prim/scal/prob/frechet_lccdf.hpp
@@ -48,6 +48,7 @@ typename return_type<T_y, T_shape, T_scale>::type frechet_lccdf(
 
   using std::exp;
   using std::log;
+  using std::pow;
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_scale> sigma_vec(sigma);
   scalar_seq_view<T_shape> alpha_vec(alpha);
diff --git a/stan/math/prim/scal/prob/frechet_lcdf.hpp b/stan/math/prim/scal/prob/frechet_lcdf.hpp
index a5fe330caf1..ce439f9bb06 100644
--- a/stan/math/prim/scal/prob/frechet_lcdf.hpp
+++ b/stan/math/prim/scal/prob/frechet_lcdf.hpp
@@ -36,6 +36,7 @@ typename return_type<T_y, T_shape, T_scale>::type frechet_lcdf(
 
   using boost::math::tools::promote_args;
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/frechet_lpdf.hpp b/stan/math/prim/scal/prob/frechet_lpdf.hpp
index 0b1ab5e96ab..bd29a80cce5 100644
--- a/stan/math/prim/scal/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/scal/prob/frechet_lpdf.hpp
@@ -36,6 +36,7 @@ typename return_type<T_y, T_shape, T_scale>::type frechet_lpdf(
       T_partials_return;
 
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/lognormal_lpdf.hpp b/stan/math/prim/scal/prob/lognormal_lpdf.hpp
index 8159e961576..c5f08e767de 100644
--- a/stan/math/prim/scal/prob/lognormal_lpdf.hpp
+++ b/stan/math/prim/scal/prob/lognormal_lpdf.hpp
@@ -59,7 +59,6 @@ typename return_type<T_y, T_loc, T_scale>::type lognormal_lpdf(
 
   operands_and_partials<T_y, T_loc, T_scale> ops_partials(y, mu, sigma);
 
-  using std::log;
   using std::log;
 
   VectorBuilder<include_summand<propto, T_scale>::value, T_partials_return,
diff --git a/stan/math/prim/scal/prob/neg_binomial_2_lpmf.hpp b/stan/math/prim/scal/prob/neg_binomial_2_lpmf.hpp
index e45cf567809..6b140aace48 100644
--- a/stan/math/prim/scal/prob/neg_binomial_2_lpmf.hpp
+++ b/stan/math/prim/scal/prob/neg_binomial_2_lpmf.hpp
@@ -54,7 +54,6 @@ typename return_type<T_location, T_precision>::type neg_binomial_2_lpmf(
     return 0.0;
 
   using std::log;
-  using std::log;
 
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_location> mu_vec(mu);
diff --git a/stan/math/prim/scal/prob/neg_binomial_lccdf.hpp b/stan/math/prim/scal/prob/neg_binomial_lccdf.hpp
index 0c177bd4ded..232c868f541 100644
--- a/stan/math/prim/scal/prob/neg_binomial_lccdf.hpp
+++ b/stan/math/prim/scal/prob/neg_binomial_lccdf.hpp
@@ -51,7 +51,6 @@ typename return_type<T_shape, T_inv_scale>::type neg_binomial_lccdf(
   scalar_seq_view<T_inv_scale> beta_vec(beta);
   size_t size = max_size(n, alpha, beta);
 
-  using std::exp;
   using std::exp;
   using std::log;
   using std::pow;
diff --git a/stan/math/prim/scal/prob/neg_binomial_lcdf.hpp b/stan/math/prim/scal/prob/neg_binomial_lcdf.hpp
index 4892d6bbf8b..37de5f111f9 100644
--- a/stan/math/prim/scal/prob/neg_binomial_lcdf.hpp
+++ b/stan/math/prim/scal/prob/neg_binomial_lcdf.hpp
@@ -51,7 +51,6 @@ typename return_type<T_shape, T_inv_scale>::type neg_binomial_lcdf(
   scalar_seq_view<T_inv_scale> beta_vec(beta);
   size_t size = max_size(n, alpha, beta);
 
-  using std::exp;
   using std::exp;
   using std::log;
   using std::pow;
diff --git a/stan/math/prim/scal/prob/neg_binomial_lpmf.hpp b/stan/math/prim/scal/prob/neg_binomial_lpmf.hpp
index 0f486094e4a..1983cbd3874 100644
--- a/stan/math/prim/scal/prob/neg_binomial_lpmf.hpp
+++ b/stan/math/prim/scal/prob/neg_binomial_lpmf.hpp
@@ -54,7 +54,6 @@ typename return_type<T_shape, T_inv_scale>::type neg_binomial_lpmf(
     return 0.0;
 
   using std::log;
-  using std::log;
 
   scalar_seq_view<T_n> n_vec(n);
   scalar_seq_view<T_shape> alpha_vec(alpha);
diff --git a/stan/math/prim/scal/prob/normal_lpdf.hpp b/stan/math/prim/scal/prob/normal_lpdf.hpp
index 9bacc746736..0028a14902c 100644
--- a/stan/math/prim/scal/prob/normal_lpdf.hpp
+++ b/stan/math/prim/scal/prob/normal_lpdf.hpp
@@ -50,7 +50,6 @@ typename return_type<T_y, T_loc, T_scale>::type normal_lpdf(
 
   using stan::is_constant_struct;
   using std::log;
-  using std::log;
 
   if (size_zero(y, mu, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/normal_sufficient_lpdf.hpp b/stan/math/prim/scal/prob/normal_sufficient_lpdf.hpp
index 7d4961ada41..cf447234cf0 100644
--- a/stan/math/prim/scal/prob/normal_sufficient_lpdf.hpp
+++ b/stan/math/prim/scal/prob/normal_sufficient_lpdf.hpp
@@ -64,6 +64,7 @@ typename return_type<T_y, T_s, T_loc, T_scale>::type normal_sufficient_lpdf(
   using stan::math::include_summand;
   using stan::math::value_of;
   using std::log;
+  using std::pow;
 
   // check if any vectors are zero length
   if (size_zero(y_bar, s_squared, n_obs, mu, sigma))
diff --git a/stan/math/prim/scal/prob/pareto_type_2_cdf.hpp b/stan/math/prim/scal/prob/pareto_type_2_cdf.hpp
index dfd1bd6daa3..00202a79129 100644
--- a/stan/math/prim/scal/prob/pareto_type_2_cdf.hpp
+++ b/stan/math/prim/scal/prob/pareto_type_2_cdf.hpp
@@ -37,6 +37,7 @@ typename return_type<T_y, T_loc, T_scale, T_shape>::type pareto_type_2_cdf(
   static const char* function = "pareto_type_2_cdf";
 
   using std::log;
+  using std::pow;
 
   T_partials_return P(1.0);
 
diff --git a/stan/math/prim/scal/prob/pareto_type_2_lcdf.hpp b/stan/math/prim/scal/prob/pareto_type_2_lcdf.hpp
index 6b2408cf9dc..1eafff04890 100644
--- a/stan/math/prim/scal/prob/pareto_type_2_lcdf.hpp
+++ b/stan/math/prim/scal/prob/pareto_type_2_lcdf.hpp
@@ -36,6 +36,7 @@ typename return_type<T_y, T_loc, T_scale, T_shape>::type pareto_type_2_lcdf(
   static const char* function = "pareto_type_2_lcdf";
 
   using std::log;
+  using std::pow;
 
   T_partials_return P(0.0);
 
diff --git a/stan/math/prim/scal/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/scal/prob/pareto_type_2_lpdf.hpp
index 95aa89f4c61..a132130fe8d 100644
--- a/stan/math/prim/scal/prob/pareto_type_2_lpdf.hpp
+++ b/stan/math/prim/scal/prob/pareto_type_2_lpdf.hpp
@@ -35,7 +35,6 @@ typename return_type<T_y, T_loc, T_scale, T_shape>::type pareto_type_2_lpdf(
           T_partials_return;
 
   using std::log;
-  using std::log;
 
   if (size_zero(y, mu, lambda, alpha))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/poisson_cdf.hpp b/stan/math/prim/scal/prob/poisson_cdf.hpp
index 7a3daa796d8..e7b7d946fb1 100644
--- a/stan/math/prim/scal/prob/poisson_cdf.hpp
+++ b/stan/math/prim/scal/prob/poisson_cdf.hpp
@@ -45,7 +45,6 @@ typename return_type<T_rate>::type poisson_cdf(const T_n& n,
   scalar_seq_view<T_rate> lambda_vec(lambda);
   size_t size = max_size(n, lambda);
 
-  using std::exp;
   using std::exp;
   using std::pow;
 
diff --git a/stan/math/prim/scal/prob/rayleigh_lpdf.hpp b/stan/math/prim/scal/prob/rayleigh_lpdf.hpp
index 31ef5a29e32..e699f83cdd3 100644
--- a/stan/math/prim/scal/prob/rayleigh_lpdf.hpp
+++ b/stan/math/prim/scal/prob/rayleigh_lpdf.hpp
@@ -33,7 +33,6 @@ typename return_type<T_y, T_scale>::type rayleigh_lpdf(const T_y& y,
 
   using stan::is_constant_struct;
   using std::log;
-  using std::log;
 
   if (size_zero(y, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/student_t_lpdf.hpp b/stan/math/prim/scal/prob/student_t_lpdf.hpp
index 277477c4144..833e8f0d448 100644
--- a/stan/math/prim/scal/prob/student_t_lpdf.hpp
+++ b/stan/math/prim/scal/prob/student_t_lpdf.hpp
@@ -86,7 +86,6 @@ typename return_type<T_y, T_dof, T_loc, T_scale>::type student_t_lpdf(
   scalar_seq_view<T_scale> sigma_vec(sigma);
   size_t N = max_size(y, nu, mu, sigma);
 
-  using std::log;
   using std::log;
 
   VectorBuilder<include_summand<propto, T_y, T_dof, T_loc, T_scale>::value,
diff --git a/stan/math/prim/scal/prob/weibull_cdf.hpp b/stan/math/prim/scal/prob/weibull_cdf.hpp
index 170fff97736..0651ac76b4c 100644
--- a/stan/math/prim/scal/prob/weibull_cdf.hpp
+++ b/stan/math/prim/scal/prob/weibull_cdf.hpp
@@ -49,6 +49,7 @@ typename return_type<T_y, T_shape, T_scale>::type weibull_cdf(
   using boost::math::tools::promote_args;
   using std::exp;
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 1.0;
diff --git a/stan/math/prim/scal/prob/weibull_lccdf.hpp b/stan/math/prim/scal/prob/weibull_lccdf.hpp
index f402e65b400..0a061aa603a 100644
--- a/stan/math/prim/scal/prob/weibull_lccdf.hpp
+++ b/stan/math/prim/scal/prob/weibull_lccdf.hpp
@@ -48,6 +48,7 @@ typename return_type<T_y, T_shape, T_scale>::type weibull_lccdf(
 
   using boost::math::tools::promote_args;
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/weibull_lcdf.hpp b/stan/math/prim/scal/prob/weibull_lcdf.hpp
index a942af044b3..2f9c40fbaef 100644
--- a/stan/math/prim/scal/prob/weibull_lcdf.hpp
+++ b/stan/math/prim/scal/prob/weibull_lcdf.hpp
@@ -49,6 +49,7 @@ typename return_type<T_y, T_shape, T_scale>::type weibull_lcdf(
   using boost::math::tools::promote_args;
   using std::exp;
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 0.0;
diff --git a/stan/math/prim/scal/prob/weibull_lpdf.hpp b/stan/math/prim/scal/prob/weibull_lpdf.hpp
index 3d5d794126b..234e76c293b 100644
--- a/stan/math/prim/scal/prob/weibull_lpdf.hpp
+++ b/stan/math/prim/scal/prob/weibull_lpdf.hpp
@@ -46,6 +46,7 @@ typename return_type<T_y, T_shape, T_scale>::type weibull_lpdf(
       T_partials_return;
 
   using std::log;
+  using std::pow;
 
   if (size_zero(y, alpha, sigma))
     return 0.0;
diff --git a/stan/math/rev/core/std_isinf.hpp b/stan/math/rev/core/std_isinf.hpp
index 006105f4ddb..dfdf2740607 100644
--- a/stan/math/rev/core/std_isinf.hpp
+++ b/stan/math/rev/core/std_isinf.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_REV_CORE_STD_ISINF_HPP
 
 #include <stan/math/prim/scal/fun/is_inf.hpp>
+#include <stan/math/rev/core/std_numeric_limits.hpp>
 #include <stan/math/rev/core/var.hpp>
 
 namespace std {
@@ -18,4 +19,13 @@ inline int isinf(const stan::math::var& a) {
 }
 
 }  // namespace std
+
+namespace stan {
+namespace math {
+
+// forwarding for ADL
+inline auto isinf(const var& a) { return std::isinf(a); }
+
+}  // namespace math
+}  // namespace stan
 #endif
diff --git a/stan/math/rev/core/std_isnan.hpp b/stan/math/rev/core/std_isnan.hpp
index 4288b316821..8c658d7f4cb 100644
--- a/stan/math/rev/core/std_isnan.hpp
+++ b/stan/math/rev/core/std_isnan.hpp
@@ -18,4 +18,13 @@ namespace std {
 inline int isnan(const stan::math::var& a) { return isnan(a.val()); }
 
 }  // namespace std
+
+namespace stan {
+namespace math {
+
+// forwarding for ADL
+inline auto isnan(const var& a) { return std::isnan(a); }
+
+}  // namespace math
+}  // namespace stan
 #endif
diff --git a/stan/math/rev/core/var.hpp b/stan/math/rev/core/var.hpp
index 124f4562410..eda85431539 100644
--- a/stan/math/rev/core/var.hpp
+++ b/stan/math/rev/core/var.hpp
@@ -7,7 +7,6 @@
 #include <boost/math/tools/config.hpp>
 #include <ostream>
 #include <vector>
-#include <complex>
 #include <string>
 #include <exception>
 
@@ -182,69 +181,6 @@ class var {
   // NOLINTNEXTLINE
   var(unsigned long x) : vi_(new vari(static_cast<double>(x))) {}  // 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));
-    } 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)));
-    } 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)));
-    } 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
diff --git a/stan/math/rev/cplx.hpp b/stan/math/rev/cplx.hpp
new file mode 100644
index 00000000000..d7bfd029a2c
--- /dev/null
+++ b/stan/math/rev/cplx.hpp
@@ -0,0 +1,7 @@
+#ifndef STAN_MATH_REV_CPLX_HPP
+#define STAN_MATH_REV_CPLX_HPP
+
+#include <stan/math/rev/cplx/var.hpp>
+#include <stan/math/rev/cplx/operator_division.hpp>
+
+#endif
diff --git a/stan/math/rev/cplx/operator_division.hpp b/stan/math/rev/cplx/operator_division.hpp
new file mode 100644
index 00000000000..dc6b46bc02c
--- /dev/null
+++ b/stan/math/rev/cplx/operator_division.hpp
@@ -0,0 +1,28 @@
+#ifndef STAN_MATH_REV_CPLX_OPERATOR_DIVISION_HPP
+#define STAN_MATH_REV_CPLX_OPERATOR_DIVISION_HPP
+
+#include <stan/math/rev/cplx/var.hpp>
+#include <stan/math/rev/scal/fun/pow.hpp>
+#include <stan/math/prim/cplx/operator_division.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the result of dividing the first argument by the second.
+ * This overload exists because libc++ uses logb and scalbn for
+ * complex division, which aren't defined for fvars.
+ *
+ * @tparam T type of complex var
+ * @param t first argument
+ * @param u second argument
+ * @return first argument divided by second argument
+ */
+inline std::complex<var> operator/(std::complex<var> const& t,
+                                   std::complex<var> const& u) {
+  return stan::math::operator_division(t, u);  // no recursion
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/rev/cplx/var.hpp b/stan/math/rev/cplx/var.hpp
new file mode 100644
index 00000000000..f678a669251
--- /dev/null
+++ b/stan/math/rev/cplx/var.hpp
@@ -0,0 +1,22 @@
+#ifndef STAN_MATH_REV_CPLX_VAR_HPP
+#define STAN_MATH_REV_CPLX_VAR_HPP
+
+#include <stan/math/rev/scal/meta/ad_promotable.hpp>
+#include <stan/math/rev/scal/meta/is_var.hpp>
+#include <stan/math/rev/core/var.hpp>
+#include <stan/math/prim/cplx/complex.hpp>
+
+namespace std {
+
+/**
+ *  std::complex<var> inherits from stan::math::complex. Details
+ *  are provided there.
+ */
+template <>
+struct complex<stan::math::var> : stan::math::complex<stan::math::var> {
+  using stan::math::complex<stan::math::var>::complex;
+};
+
+}  // namespace std
+
+#endif
diff --git a/stan/math/rev/scal.hpp b/stan/math/rev/scal.hpp
index 1b001ed20e5..287d558320b 100644
--- a/stan/math/rev/scal.hpp
+++ b/stan/math/rev/scal.hpp
@@ -1,6 +1,8 @@
 #ifndef STAN_MATH_REV_SCAL_HPP
 #define STAN_MATH_REV_SCAL_HPP
 
+#include <stan/math/rev/cplx.hpp>
+
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/scal/meta/ad_promotable.hpp>
 #include <stan/math/rev/scal/meta/is_var.hpp>
@@ -87,6 +89,7 @@
 #include <stan/math/rev/scal/fun/Phi_approx.hpp>
 #include <stan/math/rev/scal/fun/pow.hpp>
 #include <stan/math/rev/scal/fun/primitive_value.hpp>
+#include <stan/math/rev/scal/fun/proj.hpp>
 #include <stan/math/rev/scal/fun/rising_factorial.hpp>
 #include <stan/math/rev/scal/fun/round.hpp>
 #include <stan/math/rev/scal/fun/sin.hpp>
diff --git a/stan/math/rev/scal/fun/abs.hpp b/stan/math/rev/scal/fun/abs.hpp
index 2aa7271d92b..3419c4ec75c 100644
--- a/stan/math/rev/scal/fun/abs.hpp
+++ b/stan/math/rev/scal/fun/abs.hpp
@@ -1,7 +1,9 @@
 #ifndef STAN_MATH_REV_SCAL_FUN_ABS_HPP
 #define STAN_MATH_REV_SCAL_FUN_ABS_HPP
 
+#include <stan/math/rev/cplx.hpp>
 #include <stan/math/rev/scal/fun/fabs.hpp>
+#include <stan/math/rev/scal/fun/hypot.hpp>
 
 namespace stan {
 namespace math {
@@ -36,4 +38,11 @@ inline var abs(const var& a) { return fabs(a); }
 
 }  // namespace math
 }  // namespace stan
+
+namespace std {
+// constrained complex overload to forward zeroing var to std::abs
+inline stan::math::var abs(std::complex<stan::math::var> const& t) {
+  return abs(std::complex<stan::math::zeroing<stan::math::var>>(t));
+}
+}  // namespace std
 #endif
diff --git a/stan/math/rev/scal/fun/exp.hpp b/stan/math/rev/scal/fun/exp.hpp
index 1ca0706f4db..eb28b4ef42c 100644
--- a/stan/math/rev/scal/fun/exp.hpp
+++ b/stan/math/rev/scal/fun/exp.hpp
@@ -2,6 +2,8 @@
 #define STAN_MATH_REV_SCAL_FUN_EXP_HPP
 
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/scal/fun/cos.hpp>
+#include <stan/math/rev/scal/fun/sin.hpp>
 #include <cmath>
 
 namespace stan {
@@ -39,6 +41,14 @@ class exp_vari : public op_v_vari {
  */
 inline var exp(const var& a) { return var(new exp_vari(a.vi_)); }
 
+#if defined(__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ == 9
+inline std::complex<stan::math::var> exp(
+    const std::complex<stan::math::var>& a) {
+  return exp(real(a))
+         * std::complex<stan::math::var>(cos(imag(a)), sin(imag(a)));
+}
+#endif
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/scal/fun/pow.hpp b/stan/math/rev/scal/fun/pow.hpp
index 792c0702e62..8edda2cf461 100644
--- a/stan/math/rev/scal/fun/pow.hpp
+++ b/stan/math/rev/scal/fun/pow.hpp
@@ -2,6 +2,12 @@
 #define STAN_MATH_REV_SCAL_FUN_POW_HPP
 
 #include <stan/math/rev/core.hpp>
+#include <stan/math/rev/scal/fun/abs.hpp>
+#include <stan/math/rev/scal/fun/atan2.hpp>
+#include <stan/math/rev/scal/fun/cos.hpp>
+#include <stan/math/rev/scal/fun/exp.hpp>
+#include <stan/math/rev/scal/fun/log.hpp>
+#include <stan/math/rev/scal/fun/sin.hpp>
 #include <stan/math/rev/scal/fun/inv.hpp>
 #include <stan/math/rev/scal/fun/inv_sqrt.hpp>
 #include <stan/math/rev/scal/fun/inv_square.hpp>
@@ -149,6 +155,17 @@ inline var pow(double base, const var& exponent) {
   return var(new pow_dv_vari(base, exponent.vi_));
 }
 
+#if defined(__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ == 9
+inline std::complex<stan::math::var> pow(
+    const std::complex<stan::math::var>& base,
+    const std::complex<stan::math::var>& exponent) {
+  auto z = exponent * log(base);
+  auto a = real(z);
+  auto b = imag(z);
+  return exp(a) * std::complex<stan::math::var>(cos(b), sin(b));
+}
+#endif
+
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/rev/scal/fun/proj.hpp b/stan/math/rev/scal/fun/proj.hpp
new file mode 100644
index 00000000000..5955b707cb6
--- /dev/null
+++ b/stan/math/rev/scal/fun/proj.hpp
@@ -0,0 +1,27 @@
+#ifndef STAN_MATH_REV_SCAL_FUN_PROJ_HPP
+#define STAN_MATH_REV_SCAL_FUN_PROJ_HPP
+
+#include <stan/math/rev/cplx.hpp>
+#include <stan/math/rev/scal/fun/is_inf.hpp>
+#include <stan/math/prim/scal/fun/proj.hpp>
+
+namespace std {
+
+/**
+ * Return a projection of a complex number onto
+ * the Riemann sphere.
+ *
+ * This overloads an erroneous definition in
+ * libstdc++ for general types.
+ *
+ * @param[in] t Variable input.
+ * @return projection onto Riemann sphere
+ */
+
+inline std::complex<stan::math::var> proj(
+    std::complex<stan::math::var> const& t) {
+  return stan::math::proj(t);
+}
+
+}  // namespace std
+#endif
diff --git a/test/prob/chi_square/chi_square_cdf_log_test.hpp b/test/prob/chi_square/chi_square_cdf_log_test.hpp
index bb96f7c35c7..6f1b9b68825 100644
--- a/test/prob/chi_square/chi_square_cdf_log_test.hpp
+++ b/test/prob/chi_square/chi_square_cdf_log_test.hpp
@@ -60,7 +60,6 @@ class AgradCdfLogChiSquare : public AgradCdfLogTest {
       const T_y& y, const T_dof& nu, const T2&, const T3&, const T4&,
       const T5&) {
     using stan::math::gamma_p;
-    using stan::math::gamma_p;
     using std::log;
 
     return log(gamma_p(nu * 0.5, y * 0.5));
diff --git a/test/prob/chi_square/chi_square_cdf_test.hpp b/test/prob/chi_square/chi_square_cdf_test.hpp
index 05f06a979f7..57f77744949 100644
--- a/test/prob/chi_square/chi_square_cdf_test.hpp
+++ b/test/prob/chi_square/chi_square_cdf_test.hpp
@@ -58,7 +58,6 @@ class AgradCdfChiSquare : public AgradCdfTest {
       const T_y& y, const T_dof& nu, const T2&, const T3&, const T4&,
       const T5&) {
     using stan::math::gamma_p;
-    using stan::math::gamma_p;
 
     return gamma_p(nu * 0.5, y * 0.5);
   }
diff --git a/test/prob/test_fixture_cdf.hpp b/test/prob/test_fixture_cdf.hpp
index 74500587bd4..8fece70add7 100644
--- a/test/prob/test_fixture_cdf.hpp
+++ b/test/prob/test_fixture_cdf.hpp
@@ -689,7 +689,6 @@ class AgradCdfTestFixture : public ::testing::Test {
   }
 
   void test_lower_bound() {
-    using stan::math::value_of;
     using stan::math::value_of;
     const size_t N_REPEAT = 3;
     vector<double> expected_cdf;
@@ -727,7 +726,6 @@ class AgradCdfTestFixture : public ::testing::Test {
   }
 
   void test_upper_bound() {
-    using stan::math::value_of;
     using stan::math::value_of;
     const size_t N_REPEAT = 3;
     vector<double> expected_cdf;
diff --git a/test/unit/math/fwd/mat/fun/elt_multiply_test.cpp b/test/unit/math/fwd/mat/fun/elt_multiply_test.cpp
index 698f3568320..7f6730a4458 100644
--- a/test/unit/math/fwd/mat/fun/elt_multiply_test.cpp
+++ b/test/unit/math/fwd/mat/fun/elt_multiply_test.cpp
@@ -195,7 +195,6 @@ TEST(AgradFwdMatrixEltMultiply, ffd_vec_vv) {
   using stan::math::elt_multiply;
   using stan::math::fvar;
   using stan::math::vector_ffd;
-  using stan::math::vector_ffd;
 
   fvar<fvar<double> > a, b, c, d;
   a.val_.val_ = 2.0;
diff --git a/test/unit/math/fwd/mat/fun/log_sum_exp_test.cpp b/test/unit/math/fwd/mat/fun/log_sum_exp_test.cpp
index 6fecc266695..41aa218a067 100644
--- a/test/unit/math/fwd/mat/fun/log_sum_exp_test.cpp
+++ b/test/unit/math/fwd/mat/fun/log_sum_exp_test.cpp
@@ -3,7 +3,6 @@
 
 using stan::math::fvar;
 using stan::math::log_sum_exp;
-using stan::math::log_sum_exp;
 
 TEST(AgradFwdMatrixLogSumExp, vector_fd) {
   using stan::math::vector_fd;
diff --git a/test/unit/math/fwd/mat/fun/mdivide_right_spd_test.cpp b/test/unit/math/fwd/mat/fun/mdivide_right_spd_test.cpp
index 33a559e9cfe..1cf1be1989f 100644
--- a/test/unit/math/fwd/mat/fun/mdivide_right_spd_test.cpp
+++ b/test/unit/math/fwd/mat/fun/mdivide_right_spd_test.cpp
@@ -156,8 +156,6 @@ TEST(AgradFwdMatrixMdivideRightSPD, fd_exceptions) {
   using stan::math::matrix_fd;
   using stan::math::mdivide_right_spd;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_fd;
   using stan::math::row_vector_fd;
 
   matrix_fd fv1(3, 3), fv2(4, 4);
@@ -357,8 +355,6 @@ TEST(AgradFwdMatrixMdivideRightSPD, ffd_exceptions) {
   using stan::math::matrix_ffd;
   using stan::math::mdivide_right_spd;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_ffd;
   using stan::math::row_vector_ffd;
 
   matrix_ffd fv1(3, 3), fv2(4, 4);
diff --git a/test/unit/math/fwd/mat/fun/mdivide_right_tri_test.cpp b/test/unit/math/fwd/mat/fun/mdivide_right_tri_test.cpp
index e67a140bbbc..1ea1d6c6944 100644
--- a/test/unit/math/fwd/mat/fun/mdivide_right_tri_test.cpp
+++ b/test/unit/math/fwd/mat/fun/mdivide_right_tri_test.cpp
@@ -156,8 +156,6 @@ TEST(AgradFwdMatrixMdivideRightTri, fd_exceptions_lower) {
   using stan::math::matrix_fd;
   using stan::math::mdivide_right_tri;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_fd;
   using stan::math::row_vector_fd;
 
   matrix_fd fv1(3, 3), fv2(4, 4);
@@ -373,8 +371,6 @@ TEST(AgradFwdMatrixMdivideRightTri, ffd_exceptions_lower) {
   using stan::math::matrix_ffd;
   using stan::math::mdivide_right_tri;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_ffd;
   using stan::math::row_vector_ffd;
 
   matrix_ffd fv1(3, 3), fv2(4, 4);
diff --git a/test/unit/math/fwd/mat/fun/sort_test.cpp b/test/unit/math/fwd/mat/fun/sort_test.cpp
index abd8e5cf36e..9938d625bf8 100644
--- a/test/unit/math/fwd/mat/fun/sort_test.cpp
+++ b/test/unit/math/fwd/mat/fun/sort_test.cpp
@@ -10,7 +10,6 @@ using stan::math::fvar;
 
 void test_sort_asc(VEC val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   AVEC x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -31,7 +30,6 @@ void test_sort_asc(VEC val) {
 }
 void test_sort_asc3(std::vector<double> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   std::vector<fvar<fvar<double> > > x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -52,7 +50,6 @@ void test_sort_asc3(std::vector<double> val) {
 }
 void test_sort_desc(VEC val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   AVEC x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -73,7 +70,6 @@ void test_sort_desc(VEC val) {
 }
 void test_sort_desc3(VEC val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   std::vector<fvar<fvar<double> > > x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -95,7 +91,6 @@ void test_sort_desc3(VEC val) {
 template <typename T, int R, int C>
 void test_sort_asc(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   typedef Eigen::Matrix<AVAR, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -122,7 +117,6 @@ void test_sort_asc(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_asc3(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   typedef Eigen::Matrix<fvar<fvar<double> >, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -149,7 +143,6 @@ void test_sort_asc3(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_desc(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   typedef Eigen::Matrix<AVAR, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -176,7 +169,6 @@ void test_sort_desc(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_desc3(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   typedef Eigen::Matrix<fvar<fvar<double> >, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
diff --git a/test/unit/math/fwd/mat/prob/categorical_logit_test.cpp b/test/unit/math/fwd/mat/prob/categorical_logit_test.cpp
index e433004c39d..d735ed6581d 100644
--- a/test/unit/math/fwd/mat/prob/categorical_logit_test.cpp
+++ b/test/unit/math/fwd/mat/prob/categorical_logit_test.cpp
@@ -6,7 +6,6 @@
 using Eigen::Dynamic;
 using Eigen::Matrix;
 using stan::math::log_softmax;
-using stan::math::log_softmax;
 
 TEST(ProbDistributionsCategoricalLogit, fvar_double) {
   using stan::math::fvar;
diff --git a/test/unit/math/fwd/scal/fun/abs_test.cpp b/test/unit/math/fwd/scal/fun/abs_test.cpp
index ef919c5b6e1..48cbf27c79b 100644
--- a/test/unit/math/fwd/scal/fun/abs_test.cpp
+++ b/test/unit/math/fwd/scal/fun/abs_test.cpp
@@ -49,15 +49,45 @@ TEST(AgradFwdAbs, Fvar) {
   EXPECT_FLOAT_EQ(0.0, h.d_);
 }
 
+TEST(AgradFwdAbs, complexNotNullIssue123) {
+  using stan::math::fvar;
+  fvar<double> x(3, 1);
+  std::complex<fvar<double>> z = x;
+  EXPECT_TRUE(z.real().val());
+  EXPECT_FALSE(z.imag().val());
+  auto y(z.real() + 1);
+  // gradient propagates partly through imaginary <- y <- real <- x
+  z.imag(y);
+  EXPECT_TRUE(z.imag().val());
+  auto zabs(abs(z));
+  EXPECT_FLOAT_EQ(zabs.val(), 5);
+  // z^2 = x^2 + (x+1)^2 = 2x^2+2x+1 ==> 2z dz/dx = 4x + 2
+  // dz/dx = (4x + 2)/(2 z) = (4 * 3 + 2)/(2 * 5) = 1.4
+  EXPECT_FLOAT_EQ(zabs.tangent(), 1.4);
+
+  // the following are complex var left and right op parameter coercions
+  // (i.e. checks for correct instantiations)
+
+  // 8/((1+3+4i)*2) = 1/(1+1i) * (1-i)/(1-i) = (1-i)/2
+  auto q(8 / ((1 + z) * 2));
+  q += std::complex<double>(.5, 1.5);             // 0.5-0.5i + .5+1.5i = 1+1i
+  std::complex<fvar<double>> r(2 * (z + 1) / 8);  // 2*(3+4i+1)/8  = 1+1i
+  EXPECT_FLOAT_EQ(abs(q - std::complex<double>(1, 1)).val(), 0.0);
+  EXPECT_FLOAT_EQ(abs(r - std::complex<double>(1, 1)).val(), 0.0);
+
+  EXPECT_NEAR((abs(q) * abs(r)).val(), 2, 1E-8);
+  EXPECT_NEAR((abs(r) * abs(q)).val(), 2, 1E-8);
+}
+
 TEST(AgradFwdAbs, FvarFvarDouble) {
   using stan::math::fvar;
   using std::abs;
 
-  fvar<fvar<double> > x;
+  fvar<fvar<double>> x;
   x.val_.val_ = 4.0;
   x.val_.d_ = 1.0;
 
-  fvar<fvar<double> > a = abs(x);
+  fvar<fvar<double>> a = abs(x);
 
   EXPECT_FLOAT_EQ(4.0, a.val_.val_);
   EXPECT_FLOAT_EQ(1.0, a.val_.d_);
diff --git a/test/unit/math/fwd/scal/fun/binomial_coefficient_log_test.cpp b/test/unit/math/fwd/scal/fun/binomial_coefficient_log_test.cpp
index 67d9195f01e..888c318ee81 100644
--- a/test/unit/math/fwd/scal/fun/binomial_coefficient_log_test.cpp
+++ b/test/unit/math/fwd/scal/fun/binomial_coefficient_log_test.cpp
@@ -17,7 +17,6 @@ TEST(AgradFwdBinomialCoefficientLog, Fvar) {
 }
 
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarDouble) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
 
diff --git a/test/unit/math/fwd/scal/fun/value_of_test.cpp b/test/unit/math/fwd/scal/fun/value_of_test.cpp
index 90c346203c3..e7a7e8712c7 100644
--- a/test/unit/math/fwd/scal/fun/value_of_test.cpp
+++ b/test/unit/math/fwd/scal/fun/value_of_test.cpp
@@ -6,7 +6,6 @@
 TEST(AgradFwd, value_of) {
   using stan::math::fvar;
   using stan::math::value_of;
-  using stan::math::value_of;
 
   fvar<double> a = 5.0;
   EXPECT_FLOAT_EQ(5.0, value_of(a));
@@ -18,7 +17,6 @@ TEST(AgradFwd, value_of) {
 TEST(AgradFwd, value_of_nan) {
   using stan::math::fvar;
   using stan::math::value_of;
-  using stan::math::value_of;
   double nan = std::numeric_limits<double>::quiet_NaN();
 
   fvar<double> a = nan;
diff --git a/test/unit/math/mix/mat/fun/elt_multiply_test.cpp b/test/unit/math/mix/mat/fun/elt_multiply_test.cpp
index 9015dccdfe3..f87e467ccd3 100644
--- a/test/unit/math/mix/mat/fun/elt_multiply_test.cpp
+++ b/test/unit/math/mix/mat/fun/elt_multiply_test.cpp
@@ -8,7 +8,6 @@ TEST(AgradMixMatrixEltMultiply, fv_vec_vv_1stDeriv) {
   using stan::math::fvar;
   using stan::math::var;
   using stan::math::vector_fv;
-  using stan::math::vector_fv;
 
   fvar<var> a(2.0, 1.0);
   fvar<var> b(5.0, 1.0);
@@ -39,7 +38,6 @@ TEST(AgradMixMatrixEltMultiply, fv_vec_vv_2ndDeriv) {
   using stan::math::fvar;
   using stan::math::var;
   using stan::math::vector_fv;
-  using stan::math::vector_fv;
 
   fvar<var> a(2.0, 1.0);
   fvar<var> b(5.0, 1.0);
@@ -599,7 +597,6 @@ TEST(AgradMixMatrixEltMultiply, ffv_vec_vv_1stDeriv) {
   using stan::math::fvar;
   using stan::math::var;
   using stan::math::vector_ffv;
-  using stan::math::vector_ffv;
 
   fvar<fvar<var> > a(2.0, 1.0);
   fvar<fvar<var> > b(5.0, 1.0);
@@ -631,7 +628,6 @@ TEST(AgradMixMatrixEltMultiply, ffv_vec_vv_2ndDeriv_1) {
   using stan::math::fvar;
   using stan::math::var;
   using stan::math::vector_ffv;
-  using stan::math::vector_ffv;
 
   fvar<fvar<var> > a(2.0, 1.0);
   fvar<fvar<var> > b(5.0, 1.0);
@@ -659,7 +655,6 @@ TEST(AgradMixMatrixEltMultiply, ffv_vec_vv_2ndDeriv_2) {
   using stan::math::fvar;
   using stan::math::var;
   using stan::math::vector_ffv;
-  using stan::math::vector_ffv;
 
   fvar<fvar<var> > a(2.0, 1.0);
   fvar<fvar<var> > b(5.0, 1.0);
@@ -687,7 +682,6 @@ TEST(AgradMixMatrixEltMultiply, ffv_vec_vv_3rdDeriv) {
   using stan::math::fvar;
   using stan::math::var;
   using stan::math::vector_ffv;
-  using stan::math::vector_ffv;
 
   fvar<fvar<var> > a(2.0, 1.0);
   fvar<fvar<var> > b(5.0, 1.0);
diff --git a/test/unit/math/mix/mat/fun/log_sum_exp_test.cpp b/test/unit/math/mix/mat/fun/log_sum_exp_test.cpp
index 4957e4f2cfd..9504a3698f4 100644
--- a/test/unit/math/mix/mat/fun/log_sum_exp_test.cpp
+++ b/test/unit/math/mix/mat/fun/log_sum_exp_test.cpp
@@ -5,7 +5,6 @@
 
 using stan::math::fvar;
 using stan::math::log_sum_exp;
-using stan::math::log_sum_exp;
 using stan::math::var;
 
 TEST(AgradMixMatrixLogSumExp, vector_fv_1st_deriv) {
diff --git a/test/unit/math/mix/mat/fun/mdivide_right_spd_test.cpp b/test/unit/math/mix/mat/fun/mdivide_right_spd_test.cpp
index cf4ca49a0af..a4579edcd8b 100644
--- a/test/unit/math/mix/mat/fun/mdivide_right_spd_test.cpp
+++ b/test/unit/math/mix/mat/fun/mdivide_right_spd_test.cpp
@@ -448,8 +448,6 @@ TEST(AgradMixMatrixMdivideRightSPD, fv_exceptions) {
   using stan::math::matrix_fv;
   using stan::math::mdivide_right_spd;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_fv;
   using stan::math::row_vector_fv;
 
   matrix_fv fv1(3, 3), fv2(4, 4);
@@ -1389,8 +1387,6 @@ TEST(AgradMixMatrixMdivideRightSPD, ffv_exceptions) {
   using stan::math::matrix_ffv;
   using stan::math::mdivide_right_spd;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_ffv;
   using stan::math::row_vector_ffv;
 
   matrix_ffv fv1(3, 3), fv2(4, 4);
diff --git a/test/unit/math/mix/mat/fun/mdivide_right_tri_test.cpp b/test/unit/math/mix/mat/fun/mdivide_right_tri_test.cpp
index 585562d65f9..b7f10a0878a 100644
--- a/test/unit/math/mix/mat/fun/mdivide_right_tri_test.cpp
+++ b/test/unit/math/mix/mat/fun/mdivide_right_tri_test.cpp
@@ -460,8 +460,6 @@ TEST(AgradMixMatrixMdivideRightTri, fv_exceptions_lower) {
   using stan::math::matrix_fv;
   using stan::math::mdivide_right_tri;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_fv;
   using stan::math::row_vector_fv;
 
   matrix_fv fv1(3, 3), fv2(4, 4);
@@ -1437,8 +1435,6 @@ TEST(AgradMixMatrixMdivideRightTri, ffv_exceptions_lower) {
   using stan::math::matrix_ffv;
   using stan::math::mdivide_right_tri;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_ffv;
   using stan::math::row_vector_ffv;
 
   matrix_ffv fv1(3, 3), fv2(4, 4);
@@ -1963,8 +1959,6 @@ TEST(AgradMixMatrixMdivideRightTri, fv_exceptions_upper) {
   using stan::math::matrix_fv;
   using stan::math::mdivide_right_tri;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_fv;
   using stan::math::row_vector_fv;
 
   matrix_fv fv1(3, 3), fv2(4, 4);
@@ -2944,8 +2938,6 @@ TEST(AgradMixMatrixMdivideRightTri, ffv_exceptions_upper) {
   using stan::math::matrix_ffv;
   using stan::math::mdivide_right_tri;
   using stan::math::row_vector_d;
-  using stan::math::row_vector_d;
-  using stan::math::row_vector_ffv;
   using stan::math::row_vector_ffv;
 
   matrix_ffv fv1(3, 3), fv2(4, 4);
diff --git a/test/unit/math/mix/mat/fun/sort_test.cpp b/test/unit/math/mix/mat/fun/sort_test.cpp
index d02f6be9535..8de9e75402c 100644
--- a/test/unit/math/mix/mat/fun/sort_test.cpp
+++ b/test/unit/math/mix/mat/fun/sort_test.cpp
@@ -9,7 +9,6 @@ using stan::math::var;
 
 void test_sort_asc2(std::vector<double> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   std::vector<fvar<var> > x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -30,7 +29,6 @@ void test_sort_asc2(std::vector<double> val) {
 }
 void test_sort_asc4(std::vector<double> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   std::vector<fvar<fvar<var> > > x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -51,7 +49,6 @@ void test_sort_asc4(std::vector<double> val) {
 }
 void test_sort_desc2(VEC val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   std::vector<fvar<var> > x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -72,7 +69,6 @@ void test_sort_desc2(VEC val) {
 }
 void test_sort_desc4(VEC val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   std::vector<fvar<fvar<var> > > x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -94,7 +90,6 @@ void test_sort_desc4(VEC val) {
 template <typename T, int R, int C>
 void test_sort_asc2(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   typedef Eigen::Matrix<fvar<var>, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -121,7 +116,6 @@ void test_sort_asc2(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_asc4(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   typedef Eigen::Matrix<fvar<fvar<var> >, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -148,7 +142,6 @@ void test_sort_asc4(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_desc2(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   typedef Eigen::Matrix<fvar<var>, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -175,7 +168,6 @@ void test_sort_desc2(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_desc4(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   typedef Eigen::Matrix<fvar<fvar<var> >, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
diff --git a/test/unit/math/mix/mat/functor/autodiff_test.cpp b/test/unit/math/mix/mat/functor/autodiff_test.cpp
index 80b864cdeb9..044a04d78ea 100644
--- a/test/unit/math/mix/mat/functor/autodiff_test.cpp
+++ b/test/unit/math/mix/mat/functor/autodiff_test.cpp
@@ -3,6 +3,7 @@
 #include <test/unit/math/rev/mat/fun/util.hpp>
 #include <iostream>
 #include <stdexcept>
+#include <type_traits>
 #include <vector>
 
 using Eigen::Dynamic;
@@ -238,11 +239,11 @@ TEST(AgradAutoDiff, GradientHessian) {
   Matrix<double, Dynamic, Dynamic> norm_hess_analytic;
   Matrix<double, Dynamic, Dynamic> poly_hess_analytic;
 
-  std::vector<Matrix<double, Dynamic, Dynamic> > norm_grad_hess_agrad;
-  std::vector<Matrix<double, Dynamic, Dynamic> > poly_grad_hess_agrad;
+  std::vector<Matrix<double, Dynamic, Dynamic>> norm_grad_hess_agrad;
+  std::vector<Matrix<double, Dynamic, Dynamic>> poly_grad_hess_agrad;
 
-  std::vector<Matrix<double, Dynamic, Dynamic> > norm_grad_hess_analytic;
-  std::vector<Matrix<double, Dynamic, Dynamic> > poly_grad_hess_analytic;
+  std::vector<Matrix<double, Dynamic, Dynamic>> norm_grad_hess_analytic;
+  std::vector<Matrix<double, Dynamic, Dynamic>> poly_grad_hess_analytic;
 
   normal_eval_analytic = log_normal_density(normal_eval_vec);
   poly_eval_analytic = mixed_third_poly(poly_eval_vec);
@@ -278,3 +279,68 @@ TEST(AgradAutoDiff, GradientHessian) {
                         poly_grad_hess_agrad[i](j, k));
       }
 }
+
+// using declarations for static_assert checks
+// do not use anonymous namespace (unique per TU & injects names)
+using cScal = std::complex<double>;
+using cT = Eigen::VectorBlock<
+    Eigen::Block<Eigen::Matrix<std::complex<double>, 2, 2, 0, 2, 2>, 1, 2,
+                 false>,
+    -1>;
+static_assert(!(stan::is_complex<cT>::value || stan::is_arith_like<cT>::value),
+              "");
+// ensure Eigen still thinks complex VectorBlocks aren't scalar operands
+static_assert(
+    !Eigen::internal::has_ReturnType<Eigen::ScalarBinaryOpTraits<
+        cScal, cT, Eigen::internal::scalar_product_op<cScal, cT>>>::value,
+    "");
+
+TEST(AgradAutoDiff, ComplexEigenvalueOfRotationGradientHessian) {
+  auto tol([](auto d) {
+    return pow(2., -53. / 2.) * (fabs(d) + 1.0);
+  });  // tolerance
+  auto equal(
+      [tol](auto l, auto r) { return fabs(l - r) < tol(l); });  // equality
+  auto dbl([&](auto e) { return stan::math::val(e, double()); });
+
+  // return an eigenvalue of rotation by angle alpha
+  auto rotation_eigenvalue([equal, dbl](auto alpha) {
+    auto c(cos(alpha));  // cosine of rotation angle alpha
+    auto s(sin(alpha));  // sine of rotation angle alpha
+    // rotation matrix
+    auto r((Eigen::Matrix<decltype(c), 2, 2>() << c, -s, s, c).finished());
+    Eigen::EigenSolver<decltype(r)> es(r);           // complex eigensolution
+    auto vtr(es.eigenvectors().transpose().eval());  // eigenvectors v^T
+    // rv=vD; (rv)^T=(vD)^T; v^T r^T = D^T v^T = D v^T; r^T = v^T^-1 D v^T
+    auto r_check((vtr.fullPivLu().solve(es.eigenvalues().asDiagonal() * vtr))
+                     .transpose()
+                     .eval());  // ScalarType=std::complex<...>
+    auto r_err_mat((r.unaryExpr(dbl) - r_check.unaryExpr(dbl)).eval());
+    auto r_err(fabs(r_err_mat.norm()));  // check re-composition error
+    EXPECT_TRUE(equal(r_err, 0.0));      // tight tolerance residual
+    return es.eigenvalues()[0];          // return whatever eigenvalue is first
+  });                                    // end of rotation_eigenvalue
+
+  auto f([&](auto x) {  // objective critical min=0 when real = imag
+    auto rv(rotation_eigenvalue(x[0]));    // cplx eigenvalue of rotation
+    return pow(rv.real() - rv.imag(), 2);  // delta^2 of real and imag parts
+  });                                      // end of f
+
+  auto x((Eigen::VectorXd(1) << 3. * M_PI / 8.).finished());  // note: cplx eig
+  double fx(f(x));  // objective function value
+  EXPECT_FALSE(equal(fx, 0.));
+  auto rv_ini(rotation_eigenvalue(x[0]));
+  EXPECT_FALSE(equal(rv_ini.real(), rv_ini.imag()));
+  auto g((Eigen::VectorXd(1) << 1.).finished());     // gradient
+  auto h((Eigen::MatrixXd(1, 1) << 0.).finished());  // hessian
+  auto fx_old(fx + 2. * tol(fx));  // old f(x) value offset to enter loop
+  for (auto i(0u); i < 9u && !(equal(fx, fx_old) && equal(g.norm(), 0.)); i++) {
+    fx_old = fx;
+    stan::math::hessian(f, x, fx, g, h);
+    x -= h.fullPivLu().solve(g);  // newton's method for minimization
+  }
+  EXPECT_TRUE(equal(x[0], M_PI / 4.));
+  EXPECT_TRUE(equal(fx, 0.));
+  auto rv_fin(rotation_eigenvalue(x[0]));
+  EXPECT_TRUE(equal(rv_fin.real(), rv_fin.imag()));
+}
diff --git a/test/unit/math/mix/mat/prob/categorical_logit_test.cpp b/test/unit/math/mix/mat/prob/categorical_logit_test.cpp
index 03578a81e83..939c801d1de 100644
--- a/test/unit/math/mix/mat/prob/categorical_logit_test.cpp
+++ b/test/unit/math/mix/mat/prob/categorical_logit_test.cpp
@@ -6,7 +6,6 @@
 using Eigen::Dynamic;
 using Eigen::Matrix;
 using stan::math::log_softmax;
-using stan::math::log_softmax;
 
 TEST(ProbDistributionsCategoricalLogit, fvar_var) {
   using stan::math::fvar;
diff --git a/test/unit/math/mix/scal/fun/binomial_coefficient_log_test.cpp b/test/unit/math/mix/scal/fun/binomial_coefficient_log_test.cpp
index 08aa49da595..b0a93d6c318 100644
--- a/test/unit/math/mix/scal/fun/binomial_coefficient_log_test.cpp
+++ b/test/unit/math/mix/scal/fun/binomial_coefficient_log_test.cpp
@@ -117,7 +117,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarVar_FvarVar_2ndDeriv) {
 }
 
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_1stDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -144,7 +143,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_1stDeriv) {
   EXPECT_NEAR(0, g[1], 1e-8);
 }
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_Double_1stDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -168,7 +166,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_Double_1stDeriv) {
   EXPECT_FLOAT_EQ(0.69289774, g[0]);
 }
 TEST(AgradFwdBinomialCoefficientLog, Double_FvarFvarVar_1stDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -193,7 +190,6 @@ TEST(AgradFwdBinomialCoefficientLog, Double_FvarFvarVar_1stDeriv) {
 }
 
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_2ndDeriv_x) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -220,7 +216,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_2ndDeriv_x) {
   EXPECT_FLOAT_EQ(0.00099750615170258105, g[1]);
 }
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_2ndDeriv_y) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -247,7 +242,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_2ndDeriv_y) {
   EXPECT_FLOAT_EQ(-0.0019950123034051291, g[1]);
 }
 TEST(AgradFwdBinomialCoefficientLog, Double_FvarFvarVar_2ndDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -271,7 +265,6 @@ TEST(AgradFwdBinomialCoefficientLog, Double_FvarFvarVar_2ndDeriv) {
   EXPECT_NEAR(-0.00199501230340513, g[0], 1e-8);
 }
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_Double_2ndDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -295,7 +288,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_Double_2ndDeriv) {
   EXPECT_FLOAT_EQ(-0.00049862863648177515, g[0]);
 }
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_3rdDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -322,7 +314,6 @@ TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_FvarFvarVar_3rdDeriv) {
   EXPECT_FLOAT_EQ(9.9501847e-07, g[1]);
 }
 TEST(AgradFwdBinomialCoefficientLog, Double_FvarFvarVar_3rdDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
@@ -347,7 +338,6 @@ TEST(AgradFwdBinomialCoefficientLog, Double_FvarFvarVar_3rdDeriv) {
   EXPECT_NEAR(0, g[0], 1e-8);
 }
 TEST(AgradFwdBinomialCoefficientLog, FvarFvarVar_Double_3rdDeriv) {
-  using stan::math::binomial_coefficient_log;
   using stan::math::binomial_coefficient_log;
   using stan::math::fvar;
   using stan::math::var;
diff --git a/test/unit/math/prim/mat/fun/add_diag_test.cpp b/test/unit/math/prim/mat/fun/add_diag_test.cpp
index 76a2cc54fbe..f65310cf5db 100644
--- a/test/unit/math/prim/mat/fun/add_diag_test.cpp
+++ b/test/unit/math/prim/mat/fun/add_diag_test.cpp
@@ -59,29 +59,29 @@ TEST(MathPrimMat, double_mat_double_rvec_add_diag) {
 }
 
 TEST(MathPrimMat, var_mat_double_vec_add_diag) {
-  Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> mat(2, 3);
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> mat(2, 3);
   mat << 1, 1, 1, 1, 1, 1;
 
   Eigen::Matrix<double, 1, -1> to_add(2);
   to_add << 0, 1;
 
-  Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> out_mat;
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> out_mat;
   EXPECT_NO_THROW(out_mat = stan::math::add_diag(mat, to_add));
   for (int i = 0; i < 2; ++i)
-    EXPECT_FLOAT_EQ(1 + to_add[i], value_of(out_mat(i, i)))
+    EXPECT_FLOAT_EQ(1 + to_add[i], out_mat(i, i))
         << "index: ( " << i << ", " << i << ")";
 }
 
 TEST(MathPrimMat, var_mat_double_rvec_add_diag) {
-  Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> mat(2, 3);
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> mat(2, 3);
   mat << 1, 1, 1, 1, 1, 1;
 
   Eigen::Matrix<double, -1, 1> to_add(2);
   to_add << 0, 1;
 
-  Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> out_mat;
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> out_mat;
   EXPECT_NO_THROW(out_mat = stan::math::add_diag(mat, to_add));
   for (int i = 0; i < 2; ++i)
-    EXPECT_FLOAT_EQ(1 + to_add[i], value_of(out_mat(i, i)))
+    EXPECT_FLOAT_EQ(1 + to_add[i], out_mat(i, i))
         << "index: ( " << i << ", " << i << ")";
 }
diff --git a/test/unit/math/rev/arr/fun/log_sum_exp_test.cpp b/test/unit/math/rev/arr/fun/log_sum_exp_test.cpp
index 0baa401eadc..d69fbc5eecd 100644
--- a/test/unit/math/rev/arr/fun/log_sum_exp_test.cpp
+++ b/test/unit/math/rev/arr/fun/log_sum_exp_test.cpp
@@ -76,7 +76,6 @@ TEST(AgradRev, log_sum_exp_vector) {
 }
 
 TEST(AgradRev, log_sum_exp_vec_1) {
-  using stan::math::log_sum_exp;
   using stan::math::log_sum_exp;
   AVAR a(5.0);
   AVEC as = createAVEC(a);
@@ -90,7 +89,6 @@ TEST(AgradRev, log_sum_exp_vec_1) {
 }
 
 TEST(AgradRev, log_sum_exp_vec_2) {
-  using stan::math::log_sum_exp;
   using stan::math::log_sum_exp;
   AVAR a(5.0);
   AVAR b(-7.0);
@@ -119,7 +117,6 @@ TEST(AgradRev, log_sum_exp_vec_2) {
 }
 
 TEST(AgradRev, log_sum_exp_vec_3) {
-  using stan::math::log_sum_exp;
   using stan::math::log_sum_exp;
   AVAR a(5.0);
   AVAR b(-7.0);
diff --git a/test/unit/math/rev/core/var_test.cpp b/test/unit/math/rev/core/var_test.cpp
index dd174f64e6e..6c8ef1ba322 100644
--- a/test/unit/math/rev/core/var_test.cpp
+++ b/test/unit/math/rev/core/var_test.cpp
@@ -69,16 +69,35 @@ 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, complexNotNullIssue123) {
+  stan::math::start_nested();
+  AVAR x = 3;
+  std::complex<AVAR> z = x;
+  ASSERT_TRUE(z.imag().operator->() != nullptr);
+  EXPECT_TRUE(z.real().val());
+  EXPECT_FALSE(z.imag().val());
+  auto y(z.real() + 1);
+  // gradient back propagates partly through imaginary -> y -> real -> x
+  z.imag(y);
+  EXPECT_TRUE(z.imag().val());
+  auto zabs(abs(z));
+  EXPECT_FLOAT_EQ(zabs.val(), 5);
+  zabs.grad();  // z^2 = x^2 + (x+1)^2 = 2x^2+2x+1 ==> 2z dz/dx = 4x + 2
+  // dz/dx = (4x + 2)/(2 z) = (4 * 3 + 2)/(2 * 5) = 1.4
+  EXPECT_FLOAT_EQ(x.adj(), 1.4);
+
+  // the following are complex var left and right op parameter coercions
+  // (i.e. checks for correct instantiations)
+
+  // 8/((1+3+4i)*2) = 1/(1+1i) * (1-i)/(1-i) = (1-i)/2
+  auto q(8 / ((1 + z) * 2));
+  q += std::complex<double>(.5, 1.5);     // 0.5-0.5i + .5+1.5i = 1+1i
+  std::complex<AVAR> r(2 * (z + 1) / 8);  // 2*(3+4i+1)/8  = 1+1i
+  EXPECT_FLOAT_EQ(abs(q - std::complex<double>(1, 1)).val(), 0.0);
+  EXPECT_FLOAT_EQ(abs(r - std::complex<double>(1, 1)).val(), 0.0);
+  stan::math::recover_memory_nested();
 }
 
 TEST_F(AgradRev, a_eq_x) {
diff --git a/test/unit/math/rev/mat/fun/complex_eigenvalues_test.cpp b/test/unit/math/rev/mat/fun/complex_eigenvalues_test.cpp
new file mode 100644
index 00000000000..80c54b272e4
--- /dev/null
+++ b/test/unit/math/rev/mat/fun/complex_eigenvalues_test.cpp
@@ -0,0 +1,47 @@
+#include <stan/math/rev/mat.hpp>
+#include <gtest/gtest.h>
+
+struct make_I {
+  template <typename T>
+  Eigen::Matrix<T, Eigen::Dynamic, 1> operator()(
+      const Eigen::Matrix<T, Eigen::Dynamic, 1>& a) const {
+    typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
+    const int K = static_cast<int>(sqrt(a.rows()));
+    matrix_t A(K, K);
+    int pos = 0;
+    for (int i = 0; i < K; i++)
+      for (int j = 0; j < K; j++)
+        A(i, j) = a[pos++];
+    Eigen::EigenSolver<matrix_t> es(A);
+    matrix_t I = (es.eigenvectors().inverse() * A * es.eigenvectors()
+                  * es.eigenvalues().asDiagonal().inverse())
+                     .real();
+    I.resize(K * K, 1);
+    return I;
+  }
+};
+
+TEST(MathMatrix, eigen_decomposition) {
+  using stan::math::matrix_d;
+  int K = 9;
+  matrix_d A(K, K);
+  A.setRandom();
+  matrix_d J;
+  Eigen::VectorXd f_x;
+  A.resize(K * K, 1);
+  make_I f;
+  stan::math::jacobian(f, A, f_x, J);
+  const double TOL = 1e-12;
+  int pos = 0;
+  for (int j = 0; j < K; j++) {
+    for (int i = 0; i < j; i++)
+      EXPECT_NEAR(f_x[pos++], 0, TOL);
+    EXPECT_NEAR(f_x[pos++], 1, TOL);
+    for (int i = j + 1; i < K; i++)
+      EXPECT_NEAR(f_x[pos++], 0, TOL);
+  }
+  for (int i = 0; i < J.rows(); i++)
+    for (int j = 0; j < J.cols(); j++)
+      EXPECT_NEAR(J(i, j), 0, TOL);
+  EXPECT_TRUE(J.array().any());
+}
diff --git a/test/unit/math/rev/mat/fun/complex_schur_test.cpp b/test/unit/math/rev/mat/fun/complex_schur_test.cpp
new file mode 100644
index 00000000000..280ece0a0e0
--- /dev/null
+++ b/test/unit/math/rev/mat/fun/complex_schur_test.cpp
@@ -0,0 +1,46 @@
+#include <stan/math/rev/mat.hpp>
+#include <gtest/gtest.h>
+
+struct make_I {
+  template <typename T>
+  Eigen::Matrix<T, Eigen::Dynamic, 1> operator()(
+      const Eigen::Matrix<T, Eigen::Dynamic, 1>& a) const {
+    typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
+    const int K = static_cast<int>(sqrt(a.rows()));
+    matrix_t A(K, K);
+    int pos = 0;
+    for (int i = 0; i < K; i++)
+      for (int j = 0; j < K; j++)
+        A(i, j) = a[pos++];
+    Eigen::ComplexSchur<matrix_t> cs(A);
+    auto M = (cs.matrixU().adjoint() * cs.matrixU()).eval();
+    matrix_t I = M.real() + M.imag();
+    I.resize(K * K, 1);
+    return I;
+  }
+};
+
+TEST(MathMatrix, complex_schur_decomposition) {
+  using stan::math::matrix_d;
+  int K = 9;
+  matrix_d A(K, K);
+  A.setRandom();
+  matrix_d J;
+  Eigen::VectorXd f_x;
+  A.resize(K * K, 1);
+  make_I f;
+  stan::math::jacobian(f, A, f_x, J);
+  const double TOL = 7e-13;
+  int pos = 0;
+  for (int j = 0; j < K; j++) {
+    for (int i = 0; i < j; i++)
+      EXPECT_NEAR(f_x[pos++], 0, TOL);
+    EXPECT_NEAR(f_x[pos++], 1, TOL);
+    for (int i = j + 1; i < K; i++)
+      EXPECT_NEAR(f_x[pos++], 0, TOL);
+  }
+  for (int i = 0; i < J.rows(); i++)
+    for (int j = 0; j < J.cols(); j++)
+      EXPECT_NEAR(J(i, j), 0, TOL);
+  EXPECT_TRUE(J.array().any());
+}
diff --git a/test/unit/math/rev/mat/fun/multiply_test.cpp b/test/unit/math/rev/mat/fun/multiply_test.cpp
index 83c64498e88..741e39c1314 100644
--- a/test/unit/math/rev/mat/fun/multiply_test.cpp
+++ b/test/unit/math/rev/mat/fun/multiply_test.cpp
@@ -667,7 +667,6 @@ TEST(AgradRevMatrix, multiply_scalar_matrix_vc) {
 TEST(AgradRevMatrix, multiply_vector_int) {
   // test namespace resolution
   using stan::math::multiply;
-  using stan::math::multiply;
   using stan::math::vector_d;
   using stan::math::vector_v;
 
diff --git a/test/unit/math/rev/mat/fun/pseudo_eigendecomposition_test.cpp b/test/unit/math/rev/mat/fun/pseudo_eigendecomposition_test.cpp
new file mode 100644
index 00000000000..c2c3bbcae94
--- /dev/null
+++ b/test/unit/math/rev/mat/fun/pseudo_eigendecomposition_test.cpp
@@ -0,0 +1,47 @@
+#include <stan/math/rev/mat.hpp>
+#include <gtest/gtest.h>
+
+struct make_I {
+  template <typename T>
+  Eigen::Matrix<T, Eigen::Dynamic, 1> operator()(
+      const Eigen::Matrix<T, Eigen::Dynamic, 1>& a) const {
+    typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrix_t;
+    const int K = static_cast<int>(sqrt(a.rows()));
+    matrix_t A(K, K);
+    int pos = 0;
+    for (int i = 0; i < K; i++)
+      for (int j = 0; j < K; j++)
+        A(i, j) = a[pos++];
+    Eigen::EigenSolver<matrix_t> es(A);
+    matrix_t D = es.pseudoEigenvalueMatrix();
+    matrix_t V = es.pseudoEigenvectors();
+    matrix_t I = V.inverse() * A * V * D.inverse();
+    I.resize(K * K, 1);
+    return I;
+  }
+};
+
+TEST(MathMatrix, pseudo_eigen_decomposition) {
+  using stan::math::matrix_d;
+  int K = 9;
+  matrix_d A(K, K);
+  A.setRandom();
+  matrix_d J;
+  Eigen::VectorXd f_x;
+  A.resize(K * K, 1);
+  make_I f;
+  stan::math::jacobian(f, A, f_x, J);
+  const double TOL = 1e-11;
+  int pos = 0;
+  for (int j = 0; j < K; j++) {
+    for (int i = 0; i < j; i++)
+      EXPECT_NEAR(f_x[pos++], 0, TOL);
+    EXPECT_NEAR(f_x[pos++], 1, TOL);
+    for (int i = j + 1; i < K; i++)
+      EXPECT_NEAR(f_x[pos++], 0, TOL);
+  }
+  for (int i = 0; i < J.rows(); i++)
+    for (int j = 0; j < J.cols(); j++)
+      EXPECT_NEAR(J(i, j), 0, TOL);
+  EXPECT_TRUE(J.array().any());
+}
diff --git a/test/unit/math/rev/mat/fun/sort_test.cpp b/test/unit/math/rev/mat/fun/sort_test.cpp
index ffb99b68c1b..a09f1b7f88d 100644
--- a/test/unit/math/rev/mat/fun/sort_test.cpp
+++ b/test/unit/math/rev/mat/fun/sort_test.cpp
@@ -7,7 +7,6 @@
 
 void test_sort_asc(VEC val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   AVEC x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -29,7 +28,6 @@ void test_sort_asc(VEC val) {
 
 void test_sort_desc(VEC val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   AVEC x;
   for (size_t i = 0U; i < val.size(); i++)
@@ -52,7 +50,6 @@ void test_sort_desc(VEC val) {
 template <typename T, int R, int C>
 void test_sort_asc(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_asc;
-  using stan::math::sort_asc;
 
   typedef Eigen::Matrix<AVAR, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
@@ -80,7 +77,6 @@ void test_sort_asc(Eigen::Matrix<T, R, C> val) {
 template <typename T, int R, int C>
 void test_sort_desc(Eigen::Matrix<T, R, C> val) {
   using stan::math::sort_desc;
-  using stan::math::sort_desc;
 
   typedef Eigen::Matrix<AVAR, R, C> AVEC;
   typedef Eigen::Matrix<double, R, C> VEC;
diff --git a/test/unit/math/rev/scal/fun/abs_test.cpp b/test/unit/math/rev/scal/fun/abs_test.cpp
index 1789568a751..4870529381c 100644
--- a/test/unit/math/rev/scal/fun/abs_test.cpp
+++ b/test/unit/math/rev/scal/fun/abs_test.cpp
@@ -78,3 +78,13 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::abs(a));
 }
+
+TEST(AgradRev, abs_complex) {
+  std::complex<stan::math::var> z = std::complex<stan::math::var>(3, 4);
+  auto f = abs(z);
+  EXPECT_EQ(5, f.val());
+  AVEC x = createAVEC(real(z));
+  VEC g;
+  f.grad(x, g);
+  EXPECT_FLOAT_EQ(0.6, g[0]);
+}
diff --git a/test/unit/math/rev/scal/fun/acos_test.cpp b/test/unit/math/rev/scal/fun/acos_test.cpp
index 830e92d00b5..b9918e75a3a 100644
--- a/test/unit/math/rev/scal/fun/acos_test.cpp
+++ b/test/unit/math/rev/scal/fun/acos_test.cpp
@@ -74,3 +74,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::acos(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = acos(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/acosh_test.cpp b/test/unit/math/rev/scal/fun/acosh_test.cpp
index 0d70add0e8a..21f2cc2d98a 100644
--- a/test/unit/math/rev/scal/fun/acosh_test.cpp
+++ b/test/unit/math/rev/scal/fun/acosh_test.cpp
@@ -61,3 +61,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 1.3;
   test::check_varis_on_stack(stan::math::acosh(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 3.0 / 2;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = acosh(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/arg_test.cpp b/test/unit/math/rev/scal/fun/arg_test.cpp
new file mode 100644
index 00000000000..d0dc22244b5
--- /dev/null
+++ b/test/unit/math/rev/scal/fun/arg_test.cpp
@@ -0,0 +1,10 @@
+#include <stan/math/rev/scal.hpp>
+#include <gtest/gtest.h>
+
+TEST(AgradRev, arg_test) {
+  std::complex<stan::math::var> z
+      = std::complex<stan::math::var>(1, stan::math::pi());
+  auto f = arg(z);
+  using stan::math::atan2;
+  EXPECT_EQ(f.val(), atan2(std::imag(z), std::real(z)));
+}
diff --git a/test/unit/math/rev/scal/fun/asin_test.cpp b/test/unit/math/rev/scal/fun/asin_test.cpp
index 80f72444ed3..7bf3fc2b251 100644
--- a/test/unit/math/rev/scal/fun/asin_test.cpp
+++ b/test/unit/math/rev/scal/fun/asin_test.cpp
@@ -74,3 +74,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::asin(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = asin(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/asinh_test.cpp b/test/unit/math/rev/scal/fun/asinh_test.cpp
index e215c6ab5f7..af248bf0ce2 100644
--- a/test/unit/math/rev/scal/fun/asinh_test.cpp
+++ b/test/unit/math/rev/scal/fun/asinh_test.cpp
@@ -61,3 +61,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.2;
   test::check_varis_on_stack(stan::math::asinh(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = asinh(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/atan_test.cpp b/test/unit/math/rev/scal/fun/atan_test.cpp
index f1c62c404ca..f3729fe6125 100644
--- a/test/unit/math/rev/scal/fun/atan_test.cpp
+++ b/test/unit/math/rev/scal/fun/atan_test.cpp
@@ -62,3 +62,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 1;
   test::check_varis_on_stack(stan::math::atan(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = atan(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/atanh_test.cpp b/test/unit/math/rev/scal/fun/atanh_test.cpp
index 35847145ba8..1265f1b5b5a 100644
--- a/test/unit/math/rev/scal/fun/atanh_test.cpp
+++ b/test/unit/math/rev/scal/fun/atanh_test.cpp
@@ -61,3 +61,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.3;
   test::check_varis_on_stack(stan::math::atanh(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = atanh(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/conj_test.cpp b/test/unit/math/rev/scal/fun/conj_test.cpp
new file mode 100644
index 00000000000..e3497664ed3
--- /dev/null
+++ b/test/unit/math/rev/scal/fun/conj_test.cpp
@@ -0,0 +1,9 @@
+#include <stan/math/rev/scal.hpp>
+#include <gtest/gtest.h>
+
+TEST(AgradRev, conj_test) {
+  std::complex<stan::math::var> z(1, 2);
+  auto f = conj(z);
+  EXPECT_EQ(real(f).val(), 1);
+  EXPECT_EQ(imag(f).val(), -2);
+}
diff --git a/test/unit/math/rev/scal/fun/cos_test.cpp b/test/unit/math/rev/scal/fun/cos_test.cpp
index 172d18a95e8..ea127349b93 100644
--- a/test/unit/math/rev/scal/fun/cos_test.cpp
+++ b/test/unit/math/rev/scal/fun/cos_test.cpp
@@ -51,3 +51,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.49;
   test::check_varis_on_stack(stan::math::cos(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = cos(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/cosh_test.cpp b/test/unit/math/rev/scal/fun/cosh_test.cpp
index e2a83ec0990..82703960d18 100644
--- a/test/unit/math/rev/scal/fun/cosh_test.cpp
+++ b/test/unit/math/rev/scal/fun/cosh_test.cpp
@@ -67,3 +67,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::cosh(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = cosh(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/exp_test.cpp b/test/unit/math/rev/scal/fun/exp_test.cpp
index 98a9dfd3e99..d24ec8df9e8 100644
--- a/test/unit/math/rev/scal/fun/exp_test.cpp
+++ b/test/unit/math/rev/scal/fun/exp_test.cpp
@@ -30,3 +30,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a(6.0);
   test::check_varis_on_stack(stan::math::exp(a));
 }
+
+// this test used to not build on g++-4.9 and lower
+TEST(AgradRev, euler) {
+  std::complex<stan::math::var> z
+      = std::complex<stan::math::var>(0.0, stan::math::pi());
+  std::complex<stan::math::var> f = exp(z);
+  EXPECT_FLOAT_EQ(-1, real(f).val());
+  EXPECT_FLOAT_EQ(1, 1 + imag(f).val());
+  AVEC x = createAVEC(imag(z));
+  VEC g;
+  real(f).grad(x, g);
+  EXPECT_FLOAT_EQ(1 + g[0], 1);
+}
diff --git a/test/unit/math/rev/scal/fun/ibeta_test.cpp b/test/unit/math/rev/scal/fun/ibeta_test.cpp
index 13c0dcab921..74fad115a58 100644
--- a/test/unit/math/rev/scal/fun/ibeta_test.cpp
+++ b/test/unit/math/rev/scal/fun/ibeta_test.cpp
@@ -6,7 +6,6 @@
 #include <limits>
 
 TEST(AgradRev, ibeta_vvv) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
@@ -37,7 +36,6 @@ TEST(AgradRev, ibeta_vvv) {
   EXPECT_FLOAT_EQ(ibeta_derivative(a.val(), b.val(), c.val()), grad_f[2]);
 }
 TEST(AgradRev, ibeta_vvd) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
@@ -66,7 +64,6 @@ TEST(AgradRev, ibeta_vvd) {
   EXPECT_FLOAT_EQ(0.02507405, grad_f[1]);
 }
 TEST(AgradRev, ibeta_vdv) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
@@ -95,7 +92,6 @@ TEST(AgradRev, ibeta_vdv) {
   EXPECT_FLOAT_EQ(ibeta_derivative(a.val(), b, c.val()), grad_f[1]);
 }
 TEST(AgradRev, ibeta_vdd) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
@@ -122,7 +118,6 @@ TEST(AgradRev, ibeta_vdd) {
   EXPECT_FLOAT_EQ(-0.03737671, grad_f[0]);
 }
 TEST(AgradRev, ibeta_dvv) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
@@ -151,7 +146,6 @@ TEST(AgradRev, ibeta_dvv) {
   EXPECT_FLOAT_EQ(ibeta_derivative(a, b.val(), c.val()), grad_f[1]);
 }
 TEST(AgradRev, ibeta_dvd) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
@@ -178,7 +172,6 @@ TEST(AgradRev, ibeta_dvd) {
   EXPECT_FLOAT_EQ(0.02507405, grad_f[0]);
 }
 TEST(AgradRev, ibeta_ddv) {
-  using stan::math::ibeta;
   using stan::math::ibeta;
   using stan::math::var;
 
diff --git a/test/unit/math/rev/scal/fun/if_else_test.cpp b/test/unit/math/rev/scal/fun/if_else_test.cpp
index 76c88b78077..c1f678d73dc 100644
--- a/test/unit/math/rev/scal/fun/if_else_test.cpp
+++ b/test/unit/math/rev/scal/fun/if_else_test.cpp
@@ -5,7 +5,6 @@
 #include <limits>
 
 TEST(AgradRev, if_else) {
-  using stan::math::if_else;
   using stan::math::if_else;
   using stan::math::var;
 
diff --git a/test/unit/math/rev/scal/fun/log10_test.cpp b/test/unit/math/rev/scal/fun/log10_test.cpp
index 9ab86536e3e..18b463c2c43 100644
--- a/test/unit/math/rev/scal/fun/log10_test.cpp
+++ b/test/unit/math/rev/scal/fun/log10_test.cpp
@@ -29,3 +29,16 @@ TEST(AgradRev, check_varis_on_stack) {
   stan::math::var x = 4.0;
   test::check_varis_on_stack(stan::math::log10(x));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = stan::math::pi();
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = log10(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/log1m_exp_test.cpp b/test/unit/math/rev/scal/fun/log1m_exp_test.cpp
index 147878122cc..6850cf607f8 100644
--- a/test/unit/math/rev/scal/fun/log1m_exp_test.cpp
+++ b/test/unit/math/rev/scal/fun/log1m_exp_test.cpp
@@ -39,7 +39,6 @@ TEST(AgradRev, log1m_exp) {
 }
 
 TEST(AgradRev, log1m_exp_exception) {
-  using stan::math::log1m_exp;
   using stan::math::log1m_exp;
   EXPECT_NO_THROW(log1m_exp(AVAR(-3)));
   EXPECT_NO_THROW(log1m_exp(AVAR(3)));
diff --git a/test/unit/math/rev/scal/fun/log1p_exp_test.cpp b/test/unit/math/rev/scal/fun/log1p_exp_test.cpp
index 8c447f8615c..965bce8282b 100644
--- a/test/unit/math/rev/scal/fun/log1p_exp_test.cpp
+++ b/test/unit/math/rev/scal/fun/log1p_exp_test.cpp
@@ -6,7 +6,6 @@
 void test_log1p_exp(double val) {
   using stan::math::exp;
   using stan::math::log1p_exp;
-  using stan::math::log1p_exp;
   using std::exp;
 
   AVAR a(val);
diff --git a/test/unit/math/rev/scal/fun/log_test.cpp b/test/unit/math/rev/scal/fun/log_test.cpp
index 8fc555ea7f5..b9d12141323 100644
--- a/test/unit/math/rev/scal/fun/log_test.cpp
+++ b/test/unit/math/rev/scal/fun/log_test.cpp
@@ -46,3 +46,18 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a(5.0);
   test::check_varis_on_stack(stan::math::log(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = stan::math::pi();
+  std::complex<stan::math::var> z(x, 2.0 / 3);
+  EXPECT_TRUE(log(conj(z)) == conj(log(z)));
+
+  double h = 1e-8;
+  z = std::complex<stan::math::var>(x, h);
+  auto f = log(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/norm_test.cpp b/test/unit/math/rev/scal/fun/norm_test.cpp
new file mode 100644
index 00000000000..937d526ecfb
--- /dev/null
+++ b/test/unit/math/rev/scal/fun/norm_test.cpp
@@ -0,0 +1,15 @@
+#include <stan/math/rev/scal.hpp>
+#include <test/unit/math/rev/scal/fun/nan_util.hpp>
+#include <test/unit/math/rev/scal/util.hpp>
+#include <gtest/gtest.h>
+
+TEST(AgradRev, norm_test) {
+  std::complex<stan::math::var> z = std::complex<stan::math::var>(3, 4);
+  auto f = norm(z);
+  EXPECT_EQ(f.val(), 25);
+
+  AVEC x = createAVEC(imag(z));
+  VEC g;
+  f.grad(x, g);
+  EXPECT_FLOAT_EQ(g[0], 8);
+}
diff --git a/test/unit/math/rev/scal/fun/polar_test.cpp b/test/unit/math/rev/scal/fun/polar_test.cpp
new file mode 100644
index 00000000000..6560e537f71
--- /dev/null
+++ b/test/unit/math/rev/scal/fun/polar_test.cpp
@@ -0,0 +1,8 @@
+#include <stan/math/rev/scal.hpp>
+#include <gtest/gtest.h>
+
+TEST(AgradRev, polar_test) {
+  std::complex<stan::math::var> z = std::polar(1, 0);
+  auto f = arg(z);
+  EXPECT_EQ(f.val(), 0.0);
+}
diff --git a/test/unit/math/rev/scal/fun/pow_test.cpp b/test/unit/math/rev/scal/fun/pow_test.cpp
index 7406f66eda2..b391f60aa27 100644
--- a/test/unit/math/rev/scal/fun/pow_test.cpp
+++ b/test/unit/math/rev/scal/fun/pow_test.cpp
@@ -102,3 +102,10 @@ TEST(AgradRev, check_varis_on_stack) {
   test::check_varis_on_stack(stan::math::pow(3.0, b));
   test::check_varis_on_stack(stan::math::pow(a, 4.0));
 }
+
+// this test used to not build with g++-4.9 and lower
+TEST(AgradRev, complex) {
+  std::complex<stan::math::var> i(0, 1);
+  auto f = pow(i, i);
+  EXPECT_EQ(real(f).val(), exp(-stan::math::pi() / 2));
+}
diff --git a/test/unit/math/rev/scal/fun/proj_test.cpp b/test/unit/math/rev/scal/fun/proj_test.cpp
new file mode 100644
index 00000000000..c6bcfe9d741
--- /dev/null
+++ b/test/unit/math/rev/scal/fun/proj_test.cpp
@@ -0,0 +1,22 @@
+#include <stan/math/rev/scal.hpp>
+#include <gtest/gtest.h>
+
+// This does not match the example at
+// https://en.cppreference.com/w/cpp/numeric/complex/proj
+TEST(AgradRev, proj_test) {
+  std::complex<stan::math::var> z1(1, 2);
+  auto f = proj(z1);
+  EXPECT_TRUE(f == z1) << f << std::endl;
+
+  using stan::math::positive_infinity;
+  std::complex<stan::math::var> z2(positive_infinity(), -1);
+  f = proj(z2);
+  EXPECT_TRUE(f == std::complex<stan::math::var>(positive_infinity(), -0))
+      << f << std::endl;
+
+  using stan::math::negative_infinity;
+  std::complex<stan::math::var> z3(0, negative_infinity());
+  f = proj(z3);
+  EXPECT_TRUE(f == std::complex<stan::math::var>(positive_infinity(), -0))
+      << f << std::endl;
+}
diff --git a/test/unit/math/rev/scal/fun/sin_test.cpp b/test/unit/math/rev/scal/fun/sin_test.cpp
index 879d1a29167..97164b9a0f8 100644
--- a/test/unit/math/rev/scal/fun/sin_test.cpp
+++ b/test/unit/math/rev/scal/fun/sin_test.cpp
@@ -51,3 +51,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.49;
   test::check_varis_on_stack(stan::math::sin(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = sin(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/sinh_test.cpp b/test/unit/math/rev/scal/fun/sinh_test.cpp
index fb2f8ff1508..0d3e4531a88 100644
--- a/test/unit/math/rev/scal/fun/sinh_test.cpp
+++ b/test/unit/math/rev/scal/fun/sinh_test.cpp
@@ -66,3 +66,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::sinh(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = sinh(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/sqrt_test.cpp b/test/unit/math/rev/scal/fun/sqrt_test.cpp
index 305fb300dca..d84af43011e 100644
--- a/test/unit/math/rev/scal/fun/sqrt_test.cpp
+++ b/test/unit/math/rev/scal/fun/sqrt_test.cpp
@@ -62,3 +62,8 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a(5.0);
   test::check_varis_on_stack(stan::math::sqrt(a));
 }
+
+TEST(AgradRev, complex) {
+  std::complex<stan::math::var> z(4, -5);
+  EXPECT_TRUE(std::sqrt(std::conj(z)) == std::conj(std::sqrt(z)));
+}
diff --git a/test/unit/math/rev/scal/fun/tan_test.cpp b/test/unit/math/rev/scal/fun/tan_test.cpp
index 7632639bdf7..0cd48bd702b 100644
--- a/test/unit/math/rev/scal/fun/tan_test.cpp
+++ b/test/unit/math/rev/scal/fun/tan_test.cpp
@@ -53,3 +53,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::tan(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = tan(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/tanh_test.cpp b/test/unit/math/rev/scal/fun/tanh_test.cpp
index 2620da9d11d..838481218cb 100644
--- a/test/unit/math/rev/scal/fun/tanh_test.cpp
+++ b/test/unit/math/rev/scal/fun/tanh_test.cpp
@@ -66,3 +66,16 @@ TEST(AgradRev, check_varis_on_stack) {
   AVAR a = 0.68;
   test::check_varis_on_stack(stan::math::tanh(a));
 }
+
+TEST(AgradRev, complex) {
+  stan::math::var x = 2.0 / 3;
+
+  double h = 1e-8;
+  std::complex<stan::math::var> z(x, h);
+  auto f = tanh(z);
+
+  AVEC v = createAVEC(real(z));
+  VEC g;
+  real(f).grad(v, g);
+  EXPECT_FLOAT_EQ(g[0], imag(f).val() / h);
+}
diff --git a/test/unit/math/rev/scal/fun/value_of_test.cpp b/test/unit/math/rev/scal/fun/value_of_test.cpp
index 11ed379917e..2e369affc86 100644
--- a/test/unit/math/rev/scal/fun/value_of_test.cpp
+++ b/test/unit/math/rev/scal/fun/value_of_test.cpp
@@ -2,7 +2,6 @@
 #include <gtest/gtest.h>
 
 TEST(AgradRev, value_of) {
-  using stan::math::value_of;
   using stan::math::value_of;
   using stan::math::var;