diff options
| -rw-r--r-- | Eigen/src/Core/MathFunctions.h | 19 | ||||
| -rw-r--r-- | test/array_cwise.cpp | 5 |
2 files changed, 23 insertions, 1 deletions
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 4e6053b2e..dadf2cba2 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -529,7 +529,24 @@ struct expm1_impl<std::complex<RealScalar> > { 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); } }; diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp index ae0f9f97c..48ebcc88b 100644 --- a/test/array_cwise.cpp +++ b/test/array_cwise.cpp @@ -429,6 +429,11 @@ template<typename ArrayType> void array_complex(const ArrayType& m) 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)))); |