GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Srinivas Vasudevan 2019-12-12 01:56:54 +00:00 committed by Rasmus Munk Larsen
parent 381f8f3139
commit 88062b7fed
2 changed files with 23 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Eigen/src/Core/MathFunctions.h
View File

@ -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);
}
};

5
test/array_cwise.cpp
View File

@ -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))));