Tidy up and write dox.

5 files changed, 180 insertions, 293 deletions

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 576c5d25b..c00c1488c 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -454,8 +454,7 @@ template<typename Derived> class MatrixBase
     const MatrixFunctionReturnValue<Derived> sin() const;
     const MatrixSquareRootReturnValue<Derived> sqrt() const;
     const MatrixLogarithmReturnValue<Derived> log() const;
-    template <typename ExponentType>
-    const MatrixPowerReturnValue<Derived, ExponentType> pow(const ExponentType& p) const;
+    const MatrixPowerReturnValue<Derived> pow(RealScalar p) const;
 
 #ifdef EIGEN2_SUPPORT
     template<typename ProductDerived, typename Lhs, typename Rhs>
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index da79c9e0f..30d32e2dc 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -271,7 +271,7 @@ template<typename Derived> struct MatrixExponentialReturnValue;
 template<typename Derived> class MatrixFunctionReturnValue;
 template<typename Derived> class MatrixSquareRootReturnValue;
 template<typename Derived> class MatrixLogarithmReturnValue;
-template<typename Derived, typename ExponentType> class MatrixPowerReturnValue;
+template<typename Derived> class MatrixPowerReturnValue;
 
 namespace internal {
 template <typename Scalar>
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions
index 041d3b7ec..002c1f71c 100644
--- a/unsupported/Eigen/MatrixFunctions
+++ b/unsupported/Eigen/MatrixFunctions
@@ -223,8 +223,7 @@ Output: \verbinclude MatrixLogarithm.out
 Compute the matrix raised to arbitrary real power.
 
 \code
-template <typename ExponentType>
-const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const
+const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
 \endcode
 
 \param[in]  M  base of the matrix power, should be a square matrix.
@@ -247,6 +246,15 @@ diagonal and the first super-diagonal is directly computed.
 The actual work is done by the MatrixPower class, which can compute
 \f$ M^p v \f$, where \p v is another matrix with the same rows as
 \p M. The matrix \p v is set to be the identity matrix by default.
+Therefore, the expression <tt>M.pow(p) * v</tt> is specialized for
+this. No temporary storage is created for the result. The code below
+directly evaluates R-values into L-values without aliasing issue. Do
+\b NOT try to \a optimize with noalias(). It won't compile.
+\code
+v = m.pow(p) * v;
+m = m.pow(q);
+// v2.noalias() = m.pow(p) * v1; Won't compile!
+\endcode
 
 Details of the algorithm can be found in: Nicholas J. Higham and
 Lijing Lin, "A Schur-Pad&eacute; algorithm for fractional powers of a
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
index 4c9039cc5..2a46d2cc0 100644
--- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
+++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
@@ -23,10 +23,9 @@ namespace Eigen {
  *
  * \tparam MatrixType   type of the base, expected to be an instantiation
  * of the Matrix class template.
- * \tparam IsInteger    used internally to select correct specialization.
  * \tparam PlainObject  type of the multiplier.
  */
-template <typename MatrixType, int IsInteger, typename PlainObject = MatrixType>
+template<typename MatrixType, typename PlainObject = MatrixType>
 class MatrixPower
 {
   private:
@@ -65,11 +64,11 @@ class MatrixPower
      *
      * \param[out] result  \f$ A^p b \f$, as specified in the constructor.
      */
-    template <typename ResultType> void compute(ResultType& result);
+    template<typename ResultType> void compute(ResultType& result);
  
   private:
     /**
-     * \brief Compute the matrix power.
+     * \brief Compute the matrix to integral power.
      *
      * If \p b is \em fatter than \p A, it computes \f$ A^{p_{\textrm int}}
      * \f$ first, and then multiplies it with \p b. Otherwise,
@@ -77,7 +76,7 @@ class MatrixPower
      *
      * \sa computeChainProduct(ResultType&);
      */
-    template <typename ResultType>
+    template<typename ResultType>
     void computeIntPower(ResultType& result);
 
     /**
@@ -87,14 +86,14 @@ class MatrixPower
      * powering or matrix chain multiplication or solving the linear system
      * repetitively, according to which algorithm costs less.
      */
-    template <typename ResultType>
+    template<typename ResultType>
     void computeChainProduct(ResultType&);
 
     /** \brief Compute the cost of binary powering. */
     static int computeCost(RealScalar);
 
     /** \brief Solve the linear system repetitively. */
-    template <typename ResultType>
+    template<typename ResultType>
     void partialPivLuSolve(ResultType&, RealScalar);
 
     /** \brief Compute Schur decomposition of #m_A. */
@@ -170,76 +169,9 @@ class MatrixPower
     ComplexArray m_logTdiag; ///< \brief Logarithm of the main diagonal of #m_T.
 };
 
-/**
- * \internal \ingroup MatrixFunctions_Module
- * \brief Partial specialization for integral exponents.
- */
-template <typename MatrixType, typename PlainObject>
-class MatrixPower<MatrixType, 1, PlainObject>
-{
-  public:
-    /**
-     * \brief Constructor.
-     *
-     * \param[in] A  the base of the matrix power.
-     * \param[in] p  the exponent of the matrix power.
-     * \param[in] b  the multiplier.
-     */
-    MatrixPower(const MatrixType& A, int p, const PlainObject& b) :
-      m_A(A),
-      m_p(p),
-      m_b(b),
-      m_dimA(A.cols()),
-      m_dimb(b.cols())
-    { /* empty body */ }
-
-    /**
-     * \brief Compute the matrix power.
-     *
-     * If \p b is \em fatter than \p A, it computes \f$ A^p \f$ first, and
-     * then multiplies it with \p b. Otherwise, #computeChainProduct
-     * optimizes the expression.
-     *
-     * \param[out] result  \f$ A^p b \f$, as specified in the constructor.
-     *
-     * \sa computeChainProduct(ResultType&);
-     */
-    template <typename ResultType>
-    void compute(ResultType& result);
-
-  private:
-    typedef typename MatrixType::Index Index;
-
-    const MatrixType& m_A;  ///< \brief Reference to the matrix base.
-    const int m_p;          ///< \brief The integral exponent.
-    const PlainObject& m_b; ///< \brief Reference to the multiplier.
-    const Index m_dimA;     ///< \brief The dimension of #m_A, equivalent to %m_A.cols().
-    const Index m_dimb;     ///< \brief The dimension of #m_b, equivalent to %m_b.cols().
-    MatrixType m_tmp;       ///< \brief Used for temporary storage.
-
-    /**
-     * \brief Convert matrix power into chain product.
-     *
-     * This optimizes the matrix expression. It automatically chooses binary
-     * powering or matrix chain multiplication or solving the linear system
-     * repetitively, according to which algorithm costs less.
-     */
-    template <typename ResultType>
-    void computeChainProduct(ResultType& result);
-
-    /** \brief Compute the cost of binary powering. */
-    static int computeCost(int);
-
-    /** \brief Solve the linear system repetitively. */
-    template <typename ResultType>
-    void partialPivLuSolve(ResultType&, int);
-};
-
-/******* Specialized for real exponents *******/
-
-template <typename MatrixType, int IsInteger, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::compute(ResultType& result)
+template<typename MatrixType, typename PlainObject>
+template<typename ResultType>
+void MatrixPower<MatrixType,PlainObject>::compute(ResultType& result)
 {
   using std::floor;
   using std::pow;
@@ -255,19 +187,18 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::compute(ResultType& result)
     computeSchurDecomposition();
     getFractionalExponent();
     computeIntPower(result);
-    if (m_dimA == 2) {
+    if (m_dimA == 2)
       compute2x2(m_pFrac);
-    } else {
+    else
       computeBig();
-    }
     computeTmp(Scalar());
     result = m_tmp * result;
   }
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computeIntPower(ResultType& result)
+template<typename MatrixType, typename PlainObject>
+template<typename ResultType>
+void MatrixPower<MatrixType,PlainObject>::computeIntPower(ResultType& result)
 {
   MatrixType tmp;
   if (m_dimb > m_dimA) {
@@ -280,9 +211,9 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computeIntPower(ResultType&
   }
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computeChainProduct(ResultType& result)
+template<typename MatrixType, typename PlainObject>
+template<typename ResultType>
+void MatrixPower<MatrixType,PlainObject>::computeChainProduct(ResultType& result)
 {
   using std::abs;
   using std::fmod;
@@ -314,8 +245,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computeChainProduct(ResultTy
     result = m_tmp * result;
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-int MatrixPower<MatrixType,IsInteger,PlainObject>::computeCost(RealScalar p)
+template<typename MatrixType, typename PlainObject>
+int MatrixPower<MatrixType,PlainObject>::computeCost(RealScalar p)
 {
   using std::frexp;
   using std::ldexp;
@@ -329,25 +260,25 @@ int MatrixPower<MatrixType,IsInteger,PlainObject>::computeCost(RealScalar p)
   return cost;
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p)
+template<typename MatrixType, typename PlainObject>
+template<typename ResultType>
+void MatrixPower<MatrixType,PlainObject>::partialPivLuSolve(ResultType& result, RealScalar p)
 {
   const PartialPivLU<MatrixType> Asolver(m_A);
   for (; p >= RealScalar(1); p--)
     result = Asolver.solve(result);
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computeSchurDecomposition()
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::computeSchurDecomposition()
 {
   const ComplexSchur<MatrixType> schurOfA(m_A);
   m_T = schurOfA.matrixT();
   m_U = schurOfA.matrixU();
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::getFractionalExponent()
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::getFractionalExponent()
 {
   using std::pow;
   typedef Array<RealScalar, Rows, 1, ColMajor, MaxRows> RealArray;
@@ -365,9 +296,9 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::getFractionalExponent()
   }
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
+template<typename MatrixType, typename PlainObject>
 std::complex<typename MatrixType::RealScalar>
-MatrixPower<MatrixType,IsInteger,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x)
+MatrixPower<MatrixType,PlainObject>::atanh2(const ComplexScalar& y, const ComplexScalar& x)
 {
   using std::abs;
   using std::log;
@@ -380,8 +311,8 @@ MatrixPower<MatrixType,IsInteger,PlainObject>::atanh2(const ComplexScalar& y, co
     return z + z*z*z / RealScalar(3);
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::compute2x2(RealScalar p)
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::compute2x2(RealScalar p)
 {
   using std::abs;
   using std::ceil;
@@ -413,8 +344,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::compute2x2(RealScalar p)
   }
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computeBig()
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::computeBig()
 {
   using std::ldexp;
   const int digits = std::numeric_limits<RealScalar>::digits;
@@ -450,8 +381,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computeBig()
   compute2x2(m_pFrac);
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(float normIminusT)
+template<typename MatrixType, typename PlainObject>
+inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(float normIminusT)
 {
   const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f };
   int degree = 3;
@@ -461,8 +392,8 @@ inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(float no
   return degree;
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(double normIminusT)
+template<typename MatrixType, typename PlainObject>
+inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(double normIminusT)
 {
   const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2,
       1.239917516308172e-1, 1.999045567181744e-1, 2.789358995219730e-1 };
@@ -473,8 +404,8 @@ inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(double n
   return degree;
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(long double normIminusT)
+template<typename MatrixType, typename PlainObject>
+inline int MatrixPower<MatrixType,PlainObject>::getPadeDegree(long double normIminusT)
 {
 #if LDBL_MANT_DIG == 53
   const int maxPadeDegree = 7;
@@ -506,8 +437,8 @@ inline int MatrixPower<MatrixType,IsInteger,PlainObject>::getPadeDegree(long dou
       break;
   return degree;
 }
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT)
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::computePade(const int& degree, const ComplexMatrix& IminusT)
 {
   int i = degree << 1;
   m_fT = coeff(i) * IminusT;
@@ -518,8 +449,8 @@ void MatrixPower<MatrixType,IsInteger,PlainObject>::computePade(const int& degre
   m_fT += ComplexMatrix::Identity(m_dimA, m_dimA);
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-inline typename MatrixType::RealScalar MatrixPower<MatrixType,IsInteger,PlainObject>::coeff(const int& i)
+template<typename MatrixType, typename PlainObject>
+inline typename MatrixType::RealScalar MatrixPower<MatrixType,PlainObject>::coeff(const int& i)
 {
   if (i == 1)
     return -m_pFrac;
@@ -529,90 +460,79 @@ inline typename MatrixType::RealScalar MatrixPower<MatrixType,IsInteger,PlainObj
     return (m_pFrac - RealScalar(i >> 1)) / RealScalar((i - 1) << 1);
 }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(RealScalar)
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::computeTmp(RealScalar)
 { m_tmp = (m_U * m_fT * m_U.adjoint()).real(); }
 
-template <typename MatrixType, int IsInteger, typename PlainObject>
-void MatrixPower<MatrixType,IsInteger,PlainObject>::computeTmp(ComplexScalar)
+template<typename MatrixType, typename PlainObject>
+void MatrixPower<MatrixType,PlainObject>::computeTmp(ComplexScalar)
 { m_tmp = m_U * m_fT * m_U.adjoint(); }
 
-/******* Specialized for integral exponents *******/
-
-template <typename MatrixType, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,1,PlainObject>::compute(ResultType& result)
-{
-  MatrixType tmp;
-  if (m_dimb > m_dimA) {
-    tmp = MatrixType::Identity(m_dimA, m_dimA);
-    computeChainProduct(tmp);
-    result.noalias() = tmp * m_b;
-  } else {
-    result = m_b;
-    computeChainProduct(result);
-  }
-}
-
-template <typename MatrixType, typename PlainObject>
-int MatrixPower<MatrixType,1,PlainObject>::computeCost(int p)
-{
-  int cost = 0, tmp;
-  for (tmp = p; tmp; tmp >>= 1)
-    cost++;
-  for (tmp = 1; tmp <= p; tmp <<= 1)
-    if (tmp & p) cost++;
-  return cost;
-}
-
-template <typename MatrixType, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,1,PlainObject>::partialPivLuSolve(ResultType& result, int p)
+/**
+ * \ingroup MatrixFunctions_Module
+ *
+ * \brief Proxy for the matrix power multiplied by other matrix.
+ *
+ * \tparam Derived     type of the base, a matrix (expression).
+ * \tparam RhsDerived  type of the multiplier.
+ *
+ * This class holds the arguments to the matrix power until it is
+ * assigned or evaluated for some other reason (so the argument
+ * should not be changed in the meantime). It is the return type of
+ * MatrixPowerReturnValue::operator*() and related functions and most
+ * of the time this is the only way it is used.
+ */
+template<typename Derived, typename RhsDerived>
+class MatrixPowerProductReturnValue : public ReturnByValue<MatrixPowerProductReturnValue<Derived,RhsDerived> >
 {
-  const PartialPivLU<MatrixType> Asolver(m_A);
-  for(; p; p--)
-    result = Asolver.solve(result);
-}
+  private:
+    typedef typename Derived::PlainObject MatrixType;
+    typedef typename RhsDerived::PlainObject PlainObject;
+    typedef typename RhsDerived::RealScalar RealScalar;
+    typedef typename RhsDerived::Index Index;
 
-template <typename MatrixType, typename PlainObject>
-template <typename ResultType>
-void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& result)
-{
-  using std::abs;
-  int p = abs(m_p), cost = computeCost(p);
+  public:
+    /**
+     * \brief Constructor.
+     *
+     * \param[in] A  %Matrix (expression), the base of the matrix power.
+     * \param[in] p  scalar, the exponent of the matrix power.
+     * \prarm[in] b  %Matrix (expression), the multiplier.
+     */
+    MatrixPowerProductReturnValue(const Derived& A, RealScalar p, const RhsDerived& b)
+    : m_A(A), m_p(p), m_b(b) { }
 
-  if (m_p < 0) {
-    if (p * m_dimb <= cost * m_dimA) {
-      partialPivLuSolve(result, p);
-      return;
-    } else {
-      m_tmp = m_A.inverse();
-    }
-  } else {
-    m_tmp = m_A;
-  }
- 
-  while (p * m_dimb > cost * m_dimA) {
-    if (p & 1) {
-      result = m_tmp * result;
-      cost--;
+    /**
+     * \brief Compute the expression.
+     *
+     * \param[out] result  \f$ A^p b \f$ where \p A, \p p and \p bare as
+     * in the constructor.
+     */
+    template<typename ResultType>
+    inline void evalTo(ResultType& result) const
+    {
+      const MatrixType A = m_A;
+      const PlainObject b = m_b;
+      MatrixPower<MatrixType, PlainObject> mp(A, m_p, b);
+      mp.compute(result);
     }
-    m_tmp *= m_tmp;
-    cost--;
-    p >>= 1;
-  }
 
-  for (; p; p--)
-    result = m_tmp * result;
-}
+    Index rows() const { return m_b.rows(); }
+    Index cols() const { return m_b.cols(); }
+
+  private:
+    const Derived& m_A;
+    const RealScalar m_p;
+    const RhsDerived& m_b;
+    MatrixPowerProductReturnValue& operator=(const MatrixPowerProductReturnValue&);
+};
 
 /**
  * \ingroup MatrixFunctions_Module
  *
  * \brief Proxy for the matrix power of some matrix (expression).
  *
- * \tparam Derived       type of the base, a matrix (expression).
- * \tparam ExponentType  type of the exponent, a scalar.
+ * \tparam Derived  type of the base, a matrix (expression).
  *
  * This class holds the arguments to the matrix power until it is
  * assigned or evaluated for some other reason (so the argument
@@ -620,19 +540,21 @@ void MatrixPower<MatrixType,1,PlainObject>::computeChainProduct(ResultType& resu
  * MatrixBase::pow() and related functions and most of the
  * time this is the only way it is used.
  */
-template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
-: public ReturnByValue<MatrixPowerReturnValue<Derived, ExponentType> >
+template<typename Derived>
+class MatrixPowerReturnValue : public ReturnByValue<MatrixPowerReturnValue<Derived> >
 {
-  public:
+  private:
+    typedef typename Derived::RealScalar RealScalar;
     typedef typename Derived::Index Index;
 
+  public:
     /**
      * \brief Constructor.
      *
      * \param[in] A  %Matrix (expression), the base of the matrix power.
      * \param[in] p  scalar, the exponent of the matrix power.
      */
-    MatrixPowerReturnValue(const Derived& A, const ExponentType& p)
+    MatrixPowerReturnValue(const Derived& A, RealScalar p)
     : m_A(A), m_p(p) { }
 
     /**
@@ -641,22 +563,19 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
      * The %MatrixPower class can optimize \f$ A^p b \f$ computing, and
      * this method provides an elegant way to call it.
      *
-     * \param[in]  b       %Matrix (expression), the multiplier.
-     * \param[out] result  \f$ A^p b \f$ where \p A and \p p are as in
-     * the constructor.
+     * Unlike general matrix-matrix / matrix-vector product, this does
+     * \b NOT produce a temporary storage for the result. Therefore,
+     * the code below is \a already optimal:
+     * \code
+     * v = A.pow(p) * b;
+     * // v.noalias() = A.pow(p) * b; Won't compile!
+     * \endcode
+     *
+     * \param[in] b  %Matrix (expression), the multiplier.
      */
-    template <typename OtherDerived>
-    const typename OtherDerived::PlainObject operator*(const MatrixBase<OtherDerived>& b) const
-    {
-      typedef typename OtherDerived::PlainObject PlainObject;
-      const int IsInteger = NumTraits<ExponentType>::IsInteger;
-      const typename Derived::PlainObject Aevaluated = m_A.eval();
-      const PlainObject bevaluated = b.eval();
-      PlainObject result;
-      MatrixPower<Derived, IsInteger, PlainObject> mp(Aevaluated, m_p, bevaluated);
-      mp.compute(result);
-      return result;
-    }
+    template<typename RhsDerived>
+    const MatrixPowerProductReturnValue<Derived,RhsDerived> operator*(const MatrixBase<RhsDerived>& b) const
+    { return MatrixPowerProductReturnValue<Derived,RhsDerived>(m_A, m_p, b.derived()); }
 
     /**
      * \brief Compute the matrix power.
@@ -664,14 +583,13 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
      * \param[out] result  \f$ A^p \f$ where \p A and \p p are as in the
      * constructor.
      */
-    template <typename ResultType>
+    template<typename ResultType>
     inline void evalTo(ResultType& result) const
     {
       typedef typename Derived::PlainObject PlainObject;
-      const int IsInteger = NumTraits<ExponentType>::IsInteger;
-      const PlainObject Aevaluated = m_A.eval();
+      const PlainObject A = m_A;
       const PlainObject Identity = PlainObject::Identity(m_A.rows(), m_A.cols());
-      MatrixPower<PlainObject, IsInteger> mp(Aevaluated, m_p, Identity);
+      MatrixPower<PlainObject> mp(A, m_p, Identity);
       mp.compute(result);
     }
 
@@ -680,25 +598,29 @@ template<typename Derived, typename ExponentType> class MatrixPowerReturnValue
 
   private:
     const Derived& m_A;
-    const ExponentType& m_p;
-
+    const RealScalar m_p;
     MatrixPowerReturnValue& operator=(const MatrixPowerReturnValue&);
 };
 
 namespace internal {
-  template<typename Derived, typename ExponentType>
-  struct traits<MatrixPowerReturnValue<Derived, ExponentType> >
+  template<typename Derived>
+  struct traits<MatrixPowerReturnValue<Derived> >
   {
     typedef typename Derived::PlainObject ReturnType;
   };
+
+  template<typename Derived, typename RhsDerived>
+  struct traits<MatrixPowerProductReturnValue<Derived,RhsDerived> >
+  {
+    typedef typename RhsDerived::PlainObject ReturnType;
+  };
 }
 
-template <typename Derived>
-template <typename ExponentType>
-const MatrixPowerReturnValue<Derived, ExponentType> MatrixBase<Derived>::pow(const ExponentType& p) const
+template<typename Derived>
+const MatrixPowerReturnValue<Derived> MatrixBase<Derived>::pow(RealScalar p) const
 {
   eigen_assert(rows() == cols());
-  return MatrixPowerReturnValue<Derived, ExponentType>(derived(), p);
+  return MatrixPowerReturnValue<Derived>(derived(), p);
 }
 
 } // end namespace Eigen
diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp
index 3c0e4f356..d891641f4 100644
--- a/unsupported/test/matrix_power.cpp
+++ b/unsupported/test/matrix_power.cpp
@@ -9,7 +9,7 @@
 
 #include "matrix_functions.h"
 
-template <typename T>
+template<typename T>
 void test2dRotation(double tol)
 {
   Matrix<T,2,2> A, B, C;
@@ -28,7 +28,7 @@ void test2dRotation(double tol)
   }
 }
 
-template <typename T>
+template<typename T>
 void test2dHyperbolicRotation(double tol)
 {
   Matrix<std::complex<T>,2,2> A, B, C;
@@ -48,55 +48,7 @@ void test2dHyperbolicRotation(double tol)
   }
 }
 
-template <typename MatrixType>
-void testIntPowers(const MatrixType& m, double tol)
-{
-  typedef typename MatrixType::RealScalar RealScalar;
-  const MatrixType m1 = MatrixType::Random(m.rows(), m.cols());
-  const MatrixType identity = MatrixType::Identity(m.rows(), m.cols());
-  const PartialPivLU<MatrixType> solver(m1);
-  MatrixType m2, m3, m4;
-
-  m3 = m1.pow(0);
-  m4 = m1.pow(0.);
-  std::cout << "testIntPower: i = 0   error powerm = " << relerr(identity, m3) << "   " << relerr(identity, m4) << '\n';
-  VERIFY(identity == m3 && identity == m4);
-
-  m3 = m1.pow(1);
-  m4 = m1.pow(1.);
-  std::cout << "testIntPower: i = 1   error powerm = " << relerr(m1, m3) << "   " << relerr(m1, m4) << '\n';
-  VERIFY(m1 == m3 && m1 == m4);
-
-  m2.noalias() = m1 * m1;
-  m3 = m1.pow(2);
-  m4 = m1.pow(2.);
-  std::cout << "testIntPower: i = 2   error powerm = " << relerr(m2, m3) << "   " << relerr(m2, m4) << '\n';
-  VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
-
-  for (int i = 3; i <= 20; i++) {
-    m2 *= m1;
-    m3 = m1.pow(i);
-    m4 = m1.pow(RealScalar(i));
-    std::cout << "testIntPower: i = " << i << "   error powerm = " << relerr(m2, m3) << "   " << relerr (m2, m4) << '\n';
-    VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
-  }
-
-  m2 = solver.inverse();
-  m3 = m1.pow(-1);
-  m4 = m1.pow(-1.);
-  std::cout << "testIntPower: i = -1   error powerm = " << relerr(m2, m3) << "   " << relerr (m2, m4) << '\n';
-  VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
-
-  for (int i = -2; i >= -20; i--) {
-    m2 = solver.solve(m2);
-    m3 = m1.pow(i);
-    m4 = m1.pow(RealScalar(i));
-    std::cout << "testIntPower: i = " << i << "   error powerm = " << relerr(m2, m3) << "   " << relerr (m2, m4) << '\n';
-    VERIFY(m2.isApprox(m3, RealScalar(tol)) && m2.isApprox(m4, RealScalar(tol)));
-  }
-}
-
-template <typename MatrixType>
+template<typename MatrixType>
 void testExponentLaws(const MatrixType& m, double tol)
 {
   typedef typename MatrixType::RealScalar RealScalar;
@@ -129,31 +81,40 @@ void testExponentLaws(const MatrixType& m, double tol)
   }
 }
 
-template <typename MatrixType, typename VectorType>
+template<typename MatrixType, typename VectorType>
 void testMatrixVectorProduct(const MatrixType& m, const VectorType& v, double tol)
 {
   typedef typename MatrixType::RealScalar RealScalar;
   MatrixType m1;
   VectorType v1, v2, v3;
-  RealScalar pReal;
-  signed char pInt;
+  RealScalar p;
 
   for (int i = 0; i < g_repeat; i++) {
     generateTestMatrix<MatrixType>::run(m1, m.rows());
     v1 = VectorType::Random(v.rows(), v.cols());
-    pReal = internal::random<RealScalar>();
-    pInt = rand();
-    pInt >>= 2;
-
-    v2.noalias() = m1.pow(pReal).eval() * v1;
-    v3.noalias() = m1.pow(pReal) * v1;
-    std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3);
-    VERIFY(v2.isApprox(v3, RealScalar(tol)));
-
-    v2.noalias() = m1.pow(pInt).eval() * v1;
-    v3.noalias() = m1.pow(pInt) * v1;
-    std::cout << "   " << relerr(v2, v3) << '\n';
-    VERIFY(v2.isApprox(v3, RealScalar(tol)) || v2 == v3);
+    p = internal::random<RealScalar>();
+
+    v2.noalias() = m1.pow(p).eval() * v1;
+    v1 = m1.pow(p) * v1;
+    std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v1) << '\n';
+    VERIFY(v2.isApprox(v1, RealScalar(tol)));
+  }
+}
+
+template<typename MatrixType>
+void testAliasing(const MatrixType& m)
+{
+  typedef typename MatrixType::RealScalar RealScalar;
+  MatrixType m1, m2;
+  RealScalar p;
+
+  for (int i = 0; i < g_repeat; i++) {
+    generateTestMatrix<MatrixType>::run(m1, m.rows());
+    p = internal::random<RealScalar>();
+
+    m2 = m1.pow(p);
+    m1 = m1.pow(p);
+    VERIFY(m1 == m2);
   }
 }
 
@@ -166,15 +127,6 @@ void test_matrix_power()
   CALL_SUBTEST_1(test2dHyperbolicRotation<float>(1e-5));
   CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14));
 
-  CALL_SUBTEST_2(testIntPowers(Matrix2d(), 1e-13));
-  CALL_SUBTEST_7(testIntPowers(Matrix<double,3,3,RowMajor>(), 1e-13));
-  CALL_SUBTEST_3(testIntPowers(Matrix4cd(), 1e-13));
-  CALL_SUBTEST_4(testIntPowers(MatrixXd(8,8), 1e-13));
-  CALL_SUBTEST_1(testIntPowers(Matrix2f(), 1e-4));
-  CALL_SUBTEST_5(testIntPowers(Matrix3cf(), 1e-4));
-  CALL_SUBTEST_8(testIntPowers(Matrix4f(), 1e-4));
-  CALL_SUBTEST_6(testIntPowers(MatrixXf(8,8), 1e-4));
-
   CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13));
   CALL_SUBTEST_7(testExponentLaws(Matrix<double,3,3,RowMajor>(), 1e-13));
   CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13));
@@ -192,5 +144,11 @@ void test_matrix_power()
   CALL_SUBTEST_5(testMatrixVectorProduct(Matrix3cf(), Vector3cf(), 1e-4));
   CALL_SUBTEST_8(testMatrixVectorProduct(Matrix4f(), Vector4f(), 1e-4));
   CALL_SUBTEST_6(testMatrixVectorProduct(MatrixXf(8,8), VectorXf(8), 1e-4));
-  CALL_SUBTEST_10(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13));
+  CALL_SUBTEST_9(testMatrixVectorProduct(Matrix<long double,Dynamic,Dynamic>(7,7), Matrix<long double,7,9>(), 1e-13));
+
+  CALL_SUBTEST_7(testAliasing(Matrix<double,3,3,RowMajor>()));
+  CALL_SUBTEST_3(testAliasing(Matrix4cd()));
+  CALL_SUBTEST_4(testAliasing(MatrixXd(8,8)));
+  CALL_SUBTEST_5(testAliasing(Matrix3cf()));
+  CALL_SUBTEST_6(testAliasing(MatrixXf(8,8)));
 }