GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-03 17:24:11 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #876, matrix_log_compute_2x2: directly use logp1 instead of atanh2

Browse Source
This commit is contained in:
Gael Guennebaud 2014-12-08 16:28:06 +01:00
parent bea36925db
commit 77294047d6
1 changed files with 11 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
View File

@ -53,15 +53,20 @@ void matrix_log_compute_2x2(const MatrixType& A, MatrixType& result)
result(1,0) = Scalar(0);
result(1,1) = logA11;
if (A(0,0) == A(1,1)) {
Scalar y = A(1,1) - A(0,0);
if (y==Scalar(0))
{
result(0,1) = A(0,1) / A(0,0);
} else if ((abs(A(0,0)) < 0.5*abs(A(1,1))) || (abs(A(0,0)) > 2*abs(A(1,1)))) {
result(0,1) = A(0,1) * (logA11 - logA00) / (A(1,1) - A(0,0));
} else {
}
else if ((abs(A(0,0)) < 0.5*abs(A(1,1))) || (abs(A(0,0)) > 2*abs(A(1,1))))
{
result(0,1) = A(0,1) * (logA11 - logA00) / y;
}
else
{
// computation in previous branch is inaccurate if A(1,1) \approx A(0,0)
int unwindingNumber = static_cast<int>(ceil((imag(logA11 - logA00) - M_PI) / (2*M_PI)));
Scalar y = A(1,1) - A(0,0), x = A(1,1) + A(0,0);
result(0,1) = A(0,1) * (Scalar(2) * numext::atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y;
result(0,1) = A(0,1) * (numext::log1p(y/A(0,0)) + Scalar(0,2*M_PI*unwindingNumber)) / y;
}
}