diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-12-08 16:28:06 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-12-08 16:28:06 +0100 |
commit | 77294047d68ea07c90ddabddbc4cee5dd97c5237 (patch) | |
tree | ebe4d709db8bc758dbbff91c498b46b28386690a /unsupported/Eigen/src/MatrixFunctions | |
parent | bea36925dbdd753c861d650623ae2692bb9de812 (diff) |
bug #876, matrix_log_compute_2x2: directly use logp1 instead of atanh2
Diffstat (limited to 'unsupported/Eigen/src/MatrixFunctions')
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h | 17 |
1 files changed, 11 insertions, 6 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 42b60b9b1..22bfdc4ac 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -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; } } |