diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index bb553e778..3f28068ce 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -164,7 +164,7 @@ struct imag_ref_default_impl {
   typedef typename NumTraits<Scalar>::Real RealScalar;
   EIGEN_DEVICE_FUNC static inline RealScalar& run(Scalar& x) { return reinterpret_cast<RealScalar*>(&x)[1]; }
   EIGEN_DEVICE_FUNC static inline const RealScalar& run(const Scalar& x) {
-    return reinterpret_cast<RealScalar*>(&x)[1];
+    return reinterpret_cast<const RealScalar*>(&x)[1];
   }
 };
 
@@ -1541,6 +1541,25 @@ EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T exp(const T& x) {
   return exp(x);
 }
 
+// MSVC screws up some edge-cases for std::exp(complex).
+#ifdef EIGEN_COMP_MSVC
+template <typename RealScalar>
+EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex<RealScalar> exp(const std::complex<RealScalar>& x) {
+  EIGEN_USING_STD(exp);
+  // If z is (x,±∞) (for any finite x), the result is (NaN,NaN) and FE_INVALID is raised.
+  // If z is (x,NaN) (for any finite x), the result is (NaN,NaN) and FE_INVALID may be raised.
+  if ((isfinite)(real_ref(x)) && !(isfinite)(imag_ref(x))) {
+    return std::complex<RealScalar>(NumTraits<RealScalar>::quiet_NaN(), NumTraits<RealScalar>::quiet_NaN());
+  }
+  // If z is (+∞,±∞), the result is (±∞,NaN) and FE_INVALID is raised (the sign of the real part is unspecified)
+  // If z is (+∞,NaN), the result is (±∞,NaN) (the sign of the real part is unspecified)
+  if ((real_ref(x) == NumTraits<RealScalar>::infinity() && !(isfinite)(imag_ref(x)))) {
+    return std::complex<RealScalar>(NumTraits<RealScalar>::infinity(), NumTraits<RealScalar>::quiet_NaN());
+  }
+  return exp(x);
+}
+#endif
+
 #if defined(SYCL_DEVICE_ONLY)
 SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(exp, exp)
 #endif
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
index b44cb7893..78dbf207d 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -1068,40 +1068,39 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp_complex(const Pa
   typedef typename unpacket_traits<Packet>::type Scalar;
   typedef typename Scalar::value_type RealScalar;
   const RealPacket even_mask = peven_mask(a.v);
-  const Packet even_maskp = Packet(even_mask);
   const RealPacket odd_mask = pcplxflip(Packet(even_mask)).v;
 
-  Packet p0y = Packet(pand(odd_mask, a.v));
-  Packet py0 = pcplxflip(p0y);
-  Packet pyy = padd(p0y, py0);
+  // Let a = x + iy.
+  // exp(a) = exp(x) * cis(y), plus some special edge-case handling.
 
-  RealPacket sincos = psincos_float<false, RealPacket, true>(pyy.v);
-  RealPacket cossin = pcplxflip(Packet(sincos)).v;
+  // exp(x):
+  RealPacket x = pand(a.v, even_mask);
+  x = por(x, pcplxflip(Packet(x)).v);
+  RealPacket expx = pexp(x);  // exp(x);
+
+  // cis(y):
+  RealPacket y = pand(odd_mask, a.v);
+  y = por(y, pcplxflip(Packet(y)).v);
+  RealPacket cisy = psincos_float<false, RealPacket, true>(y);
+  cisy = pcplxflip(Packet(cisy)).v;  // cos(y) + i * sin(y)
 
   const RealPacket cst_pos_inf = pset1<RealPacket>(NumTraits<RealScalar>::infinity());
   const RealPacket cst_neg_inf = pset1<RealPacket>(-NumTraits<RealScalar>::infinity());
-  Packet x_is_inf = Packet(pcmp_eq(a.v, cst_pos_inf));
-  Packet x_is_minf = Packet(pcmp_eq(a.v, cst_neg_inf));
-  Packet x_is_zero = Packet(pcmp_eq(pzero(a).v, a.v));
-  Packet x_real_is_inf = pand(even_maskp, x_is_inf);
-  Packet x_real_is_minf = pand(even_maskp, x_is_minf);
-  Packet inf0 = pset1<Packet>(Scalar(NumTraits<RealScalar>::infinity(), RealScalar(0)));
-  Packet x_is_inf0 = pand(x_real_is_inf, pcplxflip(x_is_zero));
-  x_is_inf0 = por(x_is_inf0, pcplxflip(x_is_inf0));
-  Packet x_imag_goes_zero = pand(por(x_is_minf, x_is_inf), pcplxflip(x_real_is_minf));
-  Packet x_is_nan = Packet(pisnan(a.v));
-  Packet x_real_goes_zero = pand(x_is_nan, pcplxflip(x_real_is_minf));
 
-  RealPacket pexp_real = pexp(a.v);
-  Packet pexp_half = Packet(pand(even_mask, pexp_real));
-  RealPacket xexp_flip_rp = pcplxflip(pexp_half).v;
-  RealPacket xexp = padd(pexp_half.v, xexp_flip_rp);
-  Packet result(pmul(cossin, xexp));
+  // If x is -inf, we know that cossin(y) is bounded,
+  //   so the result is (0, +/-0), where the sign of the imaginary part comes
+  //   from the sign of cossin(y).
+  RealPacket cisy_sign = por(pandnot(cisy, pabs(cisy)), pset1<RealPacket>(RealScalar(1)));
+  cisy = pselect(pcmp_eq(x, cst_neg_inf), cisy_sign, cisy);
 
-  result = pselect(x_is_inf0, inf0, result);
-  result = pselect(x_real_is_minf, pzero(a), result);
-  result = pselect(x_imag_goes_zero, pzero(a), result);
-  result = pselect(x_real_goes_zero, pzero(a), result);
+  // If x is inf, and cos(y) has unknown sign (y is inf or NaN), the result
+  // is (+/-inf, NaN), where the signs are undetermined (take the sign of y).
+  RealPacket y_sign = por(pandnot(y, pabs(y)), pset1<RealPacket>(RealScalar(1)));
+  cisy = pselect(pand(pcmp_eq(x, cst_pos_inf), pisnan(cisy)), pand(y_sign, even_mask), cisy);
+  Packet result = Packet(pmul(expx, cisy));
+
+  // If y is +/- 0, the input is real, so take the real result for consistency.
+  result = pselect(Packet(pcmp_eq(y, pzero(y))), Packet(por(pand(expx, even_mask), pand(y, odd_mask))), result);
 
   return result;
 }
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index c05f87364..92e72cb6e 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -277,6 +277,7 @@ struct packetmath_pcast_ops_runner<Scalar, Packet, std::enable_if_t<NumTraits<Sc
 
 template <typename Scalar, typename Packet>
 void packetmath_boolean_mask_ops() {
+  using RealScalar = typename NumTraits<Scalar>::Real;
   const int PacketSize = internal::unpacket_traits<Packet>::size;
   const int size = 2 * PacketSize;
   EIGEN_ALIGN_MAX Scalar data1[size];
@@ -289,7 +290,7 @@ void packetmath_boolean_mask_ops() {
   CHECK_CWISE1(internal::ptrue, internal::ptrue);
   CHECK_CWISE2_IF(true, internal::pandnot, internal::pandnot);
   for (int i = 0; i < PacketSize; ++i) {
-    data1[i] = Scalar(i);
+    data1[i] = Scalar(RealScalar(i));
     data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
   }
 
@@ -1335,6 +1336,62 @@ void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval) {
 template <typename Scalar, typename Packet, bool HasExp = internal::packet_traits<Scalar>::HasExp>
 struct exp_complex_test_impl {
   typedef typename Scalar::value_type RealScalar;
+
+  static Scalar pexp1(const Scalar& x) {
+    Packet px = internal::pset1<Packet>(x);
+    Packet py = internal::pexp(px);
+    return internal::pfirst(py);
+  }
+
+  static Scalar cis(const RealScalar& x) { return Scalar(numext::cos(x), numext::sin(x)); }
+
+  // Verify equality with signed zero.
+  static bool is_exactly_equal(const RealScalar& a, const RealScalar& b) {
+    // NaNs are always unsigned, and always compare not equal directly.
+    if ((numext::isnan)(a)) {
+      return (numext::isnan)(b);
+    }
+    // Signed zero.
+    RealScalar zero(0);
+    if (a == zero) {
+      // Signs are either 0 or NaN, so verify that their comparisons to zero are equal.
+      return (a == b) && ((numext::signbit(a) == zero) == (numext::signbit(b) == zero));
+    }
+    // Allow _some_ tolerance.
+    return verifyIsApprox(a, b);
+  }
+
+  // Verify equality with signed zero.
+  static bool is_exactly_equal(const Scalar& a, const Scalar& b) {
+    bool result = is_exactly_equal(numext::real_ref(a), numext::real_ref(b)) &&
+                  is_exactly_equal(numext::imag_ref(a), numext::imag_ref(b));
+    if (!result) {
+      std::cout << a << " != " << b << std::endl;
+    }
+    return result;
+  }
+
+  static bool is_sign_exp_unspecified(const Scalar& z) {
+    const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
+    // If z is (-∞,±∞), the result is (±0,±0) (signs are unspecified)
+    if (numext::real_ref(z) == -inf && (numext::isinf)(numext::imag_ref(z))) {
+      return true;
+    }
+    // If z is (+∞,±∞), the result is (±∞,NaN) and FE_INVALID is raised (the sign of the real part is unspecified)
+    if (numext::real_ref(z) == +inf && (numext::isinf)(numext::imag_ref(z))) {
+      return true;
+    }
+    // If z is (-∞,NaN), the result is (±0,±0) (signs are unspecified)
+    if (numext::real_ref(z) == -inf && (numext::isnan)(numext::imag_ref(z))) {
+      return true;
+    }
+    // If z is (+∞,NaN), the result is (±∞,NaN) (the sign of the real part is unspecified)
+    if (numext::real_ref(z) == +inf && (numext::isnan)(numext::imag_ref(z))) {
+      return true;
+    }
+    return false;
+  }
+
   static void run(Scalar* data1, Scalar* data2, Scalar* ref, int size) {
     const int PacketSize = internal::unpacket_traits<Packet>::size;
 
@@ -1343,27 +1400,45 @@ struct exp_complex_test_impl {
     }
     CHECK_CWISE1_N(std::exp, internal::pexp, size);
 
-    // Test misc. corner cases.
-    const RealScalar zero = RealScalar(0);
-    const RealScalar one = RealScalar(1);
-    const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
-    const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN();
-    for (RealScalar x : {zero, one, inf}) {
-      for (RealScalar y : {zero, one, inf}) {
-        data1[0] = Scalar(x, y);
-        data1[1] = Scalar(-x, y);
-        data1[2] = Scalar(x, -y);
-        data1[3] = Scalar(-x, -y);
-        CHECK_CWISE1_N(std::exp, internal::pexp, 4);
+    // Test all corner cases (and more).
+    const RealScalar edges[] = {RealScalar(0),
+                                RealScalar(1),
+                                RealScalar(2),
+                                RealScalar(EIGEN_PI / 2),
+                                RealScalar(EIGEN_PI),
+                                RealScalar(3 * EIGEN_PI / 2),
+                                RealScalar(2 * EIGEN_PI),
+                                numext::log(NumTraits<RealScalar>::highest()) - 1,
+                                NumTraits<RealScalar>::highest(),
+                                std::numeric_limits<RealScalar>::infinity(),
+                                std::numeric_limits<RealScalar>::quiet_NaN(),
+                                -RealScalar(0),
+                                -RealScalar(1),
+                                -RealScalar(2),
+                                -RealScalar(EIGEN_PI / 2),
+                                -RealScalar(EIGEN_PI),
+                                -RealScalar(3 * EIGEN_PI / 2),
+                                -RealScalar(2 * EIGEN_PI),
+                                -numext::log(NumTraits<RealScalar>::highest()) + 1,
+                                -NumTraits<RealScalar>::highest(),
+                                -std::numeric_limits<RealScalar>::infinity(),
+                                -std::numeric_limits<RealScalar>::quiet_NaN()};
+
+    for (RealScalar x : edges) {
+      for (RealScalar y : edges) {
+        Scalar z = Scalar(x, y);
+        Scalar w = pexp1(z);
+        if (is_sign_exp_unspecified(z)) {
+          Scalar abs_w = Scalar(numext::abs(numext::real_ref(w)), numext::abs(numext::imag_ref(w)));
+          Scalar expected = numext::exp(z);
+          Scalar abs_expected =
+              Scalar(numext::abs(numext::real_ref(expected)), numext::abs(numext::imag_ref(expected)));
+          VERIFY(is_exactly_equal(abs_w, abs_expected));
+        } else {
+          VERIFY(is_exactly_equal(w, numext::exp(z)));
+        }
       }
     }
-    for (RealScalar x : {zero, one, inf}) {
-      data1[0] = Scalar(x, nan);
-      data1[1] = Scalar(-x, nan);
-      data1[2] = Scalar(nan, x);
-      data1[3] = Scalar(nan, -x);
-      CHECK_CWISE1_N(std::exp, internal::pexp, 4);
-    }
   }
 };