Optimize getting exponent from IEEE floating points.

3 files changed, 24 insertions, 15 deletions

diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index fdb1dc6e2..34944e055 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h
@@ -353,7 +353,7 @@ struct check_transpose_aliasing_run_time_selector
 {
   static bool run(const Scalar* dest, const OtherDerived& src)
   {
-    return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src));
+    return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src));
   }
 };
 
@@ -362,8 +362,8 @@ struct check_transpose_aliasing_run_time_selector<Scalar,DestIsTransposed,CwiseB
 {
   static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src)
   {
-    return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src.lhs())))
-        || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src.rhs())));
+    return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.lhs())))
+        || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.rhs())));
   }
 };
 
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index acc5576fe..24c78b4b2 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -743,7 +743,16 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
 //   RealScalar e2 = abs2(subdiag[end-1]);
 //   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
   // This explain the following, somewhat more complicated, version:
-  RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e));
+  RealScalar mu = diag[end];
+  if(td==0)
+    mu -= abs(e);
+  else
+  {
+    RealScalar e2 = abs2(subdiag[end-1]);
+    RealScalar h = hypot(td,e);
+    if(e2==0)  mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
+    else       mu -= e2 / (td + (td>0 ? h : -h));
+  }
   
   RealScalar x = diag[start] - mu;
   RealScalar z = subdiag[start];
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
index 642916764..9947d9007 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h
@@ -13,12 +13,7 @@
 
 #include "StemFunction.h"
 
-namespace Eigen { 
-
-#if defined(_MSC_VER) || defined(__FreeBSD__)
-  template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
-#endif
-
+namespace Eigen {
 
 /** \ingroup MatrixFunctions_Module
   * \brief Class for computing the matrix exponential.
@@ -297,7 +292,8 @@ void MatrixExponential<MatrixType>::computeUV(float)
     pade5(m_M);
   } else {
     const float maxnorm = 3.925724783138660f;
-    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    frexp(m_l1norm / maxnorm, &m_squarings);
+    if (m_squarings < 0) m_squarings = 0;
     MatrixType A = m_M / pow(Scalar(2), m_squarings);
     pade7(A);
   }
@@ -319,7 +315,8 @@ void MatrixExponential<MatrixType>::computeUV(double)
     pade9(m_M);
   } else {
     const double maxnorm = 5.371920351148152;
-    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    frexp(m_l1norm / maxnorm, &m_squarings);
+    if (m_squarings < 0) m_squarings = 0;
     MatrixType A = m_M / pow(Scalar(2), m_squarings);
     pade13(A);
   }
@@ -344,7 +341,8 @@ void MatrixExponential<MatrixType>::computeUV(long double)
     pade9(m_M);
   } else {
     const long double maxnorm = 4.0246098906697353063L;
-    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    frexp(m_l1norm / maxnorm, &m_squarings);
+    if (m_squarings < 0) m_squarings = 0;
     MatrixType A = m_M / pow(Scalar(2), m_squarings);
     pade13(A);
   }
@@ -361,7 +359,8 @@ void MatrixExponential<MatrixType>::computeUV(long double)
     pade13(m_M);
   } else {
     const long double maxnorm = 3.2579440895405400856599663723517L;
-    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    frexp(m_l1norm / maxnorm, &m_squarings);
+    if (m_squarings < 0) m_squarings = 0;
     MatrixType A = m_M / pow(Scalar(2), m_squarings);
     pade17(A);
   }
@@ -378,7 +377,8 @@ void MatrixExponential<MatrixType>::computeUV(long double)
     pade13(m_M);
   } else {
     const long double maxnorm = 2.884233277829519311757165057717815L;
-    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    frexp(m_l1norm / maxnorm, &m_squarings);
+    if (m_squarings < 0) m_squarings = 0;
     MatrixType A = m_M / pow(Scalar(2), m_squarings);
     pade17(A);
   }