Fix implementation of complex expm1. Add tests that fail with previous implementation, but pass with the current one.

Eigen/src/Core/MathFunctions.h

diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 4e6053b..dadf2cb 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -529,7 +529,24 @@
   EIGEN_DEVICE_FUNC static inline std::complex<RealScalar> run(
       const std::complex<RealScalar>& x) {
     EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
-    return std_fallback::expm1(x);
+    RealScalar xr = x.real();
+    RealScalar xi = x.imag();
+    // expm1(z) = exp(z) - 1
+    //          = exp(x +  i * y) - 1
+    //          = exp(x) * (cos(y) + i * sin(y)) - 1
+    //          = exp(x) * cos(y) - 1 + i * exp(x) * sin(y)
+    // Imag(expm1(z)) = exp(x) * sin(y)
+    // Real(expm1(z)) = exp(x) * cos(y) - 1
+    //          = exp(x) * cos(y) - 1.
+    //          = expm1(x) + exp(x) * (cos(y) - 1)
+    //          = expm1(x) + exp(x) * (2 * sin(y / 2) ** 2)
+    RealScalar erm1 = std_fallback::expm1(xr);
+    RealScalar er = erm1 + RealScalar(1.);
+    RealScalar sin2 = sin(xi / RealScalar(2.));
+    sin2 = sin2 * sin2;
+    RealScalar s = sin(xi);
+    RealScalar real_part = erm1 - RealScalar(2.) * er * sin2;
+    return std::complex<RealScalar>(real_part, er * s);
   }
 };
 

test/array_cwise.cpp

diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp
index ae0f9f9..48ebcc8 100644
--- a/test/array_cwise.cpp
+++ b/test/array_cwise.cpp
@@ -429,6 +429,11 @@
   VERIFY_IS_APPROX(m1.exp(), exp(m1));
   VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp());
 
+  VERIFY_IS_APPROX(m1.expm1(), expm1(m1));
+  VERIFY_IS_APPROX(expm1(m1), exp(m1) - 1.);
+  // Check for larger magnitude complex numbers that expm1 matches exp - 1.
+  VERIFY_IS_APPROX(expm1(10. * m1), exp(10. * m1) - 1.);
+
   VERIFY_IS_APPROX(sinh(m1), 0.5*(exp(m1)-exp(-m1)));
   VERIFY_IS_APPROX(cosh(m1), 0.5*(exp(m1)+exp(-m1)));
   VERIFY_IS_APPROX(tanh(m1), (0.5*(exp(m1)-exp(-m1)))/(0.5*(exp(m1)+exp(-m1))));