Avoid inefficient 2x2 LU. Move atanh to internal for maintainability.

2 files changed, 3 insertions, 36 deletions

diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
index 18bcf3d0d..e1e5b770c 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h
@@ -51,7 +51,6 @@ private:
 
   void compute2x2(const MatrixType& A, MatrixType& result);
   void computeBig(const MatrixType& A, MatrixType& result);
-  static Scalar atanh2(Scalar y, Scalar x);
   int getPadeDegree(float normTminusI);
   int getPadeDegree(double normTminusI);
   int getPadeDegree(long double normTminusI);
@@ -93,20 +92,6 @@ MatrixType MatrixLogarithmAtomic<MatrixType>::compute(const MatrixType& A)
   return result;
 }
 
-/** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */
-template <typename MatrixType>
-typename MatrixType::Scalar MatrixLogarithmAtomic<MatrixType>::atanh2(Scalar y, Scalar x)
-{
-  using std::abs;
-  using std::sqrt;
-
-  Scalar z = y / x;
-  if (abs(z) > sqrt(NumTraits<Scalar>::epsilon()))
-    return Scalar(0.5) * log((x + y) / (x - y));
-  else
-    return z + z*z*z / Scalar(3);
-}
-
 /** \brief Compute logarithm of 2x2 triangular matrix. */
 template <typename MatrixType>
 void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixType& result)
@@ -131,7 +116,7 @@ void MatrixLogarithmAtomic<MatrixType>::compute2x2(const MatrixType& A, MatrixTy
     // 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) * atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y;
+    result(0,1) = A(0,1) * (Scalar(2) * internal::atanh2(y,x) + Scalar(0,2*M_PI*unwindingNumber)) / y;
   }
 }
 
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 2a46d2cc0..7238501ed 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -111,9 +111,6 @@ class MatrixPower
      */
     void getFractionalExponent();
 
-    /** \brief Compute atanh (inverse hyperbolic tangent) for \f$ y / x \f$. */
-    static ComplexScalar atanh2(const ComplexScalar& y, const ComplexScalar& x);
-
     /** \brief Compute power of 2x2 triangular matrix. */
     void compute2x2(RealScalar p);
 
@@ -223,7 +220,7 @@ void MatrixPower<MatrixType,PlainObject>::computeChainProduct(ResultType& result
   int cost = computeCost(p);
 
   if (m_pInt < RealScalar(0)) {
-    if (p * m_dimb <= cost * m_dimA) {
+    if (p * m_dimb <= cost * m_dimA && m_dimA > 2) {
       partialPivLuSolve(result, p);
       return;
     } else {
@@ -297,21 +294,6 @@ void MatrixPower<MatrixType,PlainObject>::getFractionalExponent()
 }
 
 template<typename MatrixType, typename PlainObject>
-std::complex<typename MatrixType::RealScalar>
-MatrixPower<MatrixType,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x)
-{
-  using std::abs;
-  using std::log;
-  using std::sqrt;
-  const ComplexScalar z = y / x;
-
-  if (abs(z) > sqrt(NumTraits<RealScalar>::epsilon()))
-    return RealScalar(0.5) * log((x + y) / (x - y));
-  else
-    return z + z*z*z / RealScalar(3);
-}
-
-template<typename MatrixType, typename PlainObject>
 void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
 {
   using std::abs;
@@ -337,7 +319,7 @@ void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
     } else {
       // computation in previous branch is inaccurate if abs(m_T(j,j)) \approx abs(m_T(i,i))
       unwindingNumber = ceil((imag(m_logTdiag[j] - m_logTdiag[i]) - M_PI) / (2 * M_PI));
-      w = atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber);
+      w = internal::atanh2(m_T(j,j) - m_T(i,i), m_T(j,j) + m_T(i,i)) + ComplexScalar(0, M_PI * unwindingNumber);
       m_fT(i,j) = m_T(i,j) * RealScalar(2) * exp(RealScalar(0.5) * p * (m_logTdiag[j] + m_logTdiag[i])) *
 	  sinh(p * w) / (m_T(j,j) - m_T(i,i));
     }