MATH-435

Added "power" method in "RealMatrix" interface and default implementation in "AbstractRealMatrix".


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1083698 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/linear/AbstractRealMatrix.java

@@ -17,6 +17,8 @@
 
 package org.apache.commons.math.linear;
 
+import java.util.ArrayList;
+
 import org.apache.commons.math.exception.NoDataException;
 import org.apache.commons.math.exception.NotStrictlyPositiveException;
 import org.apache.commons.math.exception.DimensionMismatchException;
@@ -147,6 +149,66 @@ public abstract class AbstractRealMatrix implements RealMatrix {
     /** {@inheritDoc} */
     public RealMatrix preMultiply(final RealMatrix m) {
         return m.multiply(this);
+    }    
+
+    /** {@inheritDoc} */
+    @Override
+    public RealMatrix power(final int p) {
+        if (p < 0) {
+            throw new IllegalArgumentException("p must be >= 0");
+        }
+        
+        if (!isSquare()) {
+            throw new NonSquareMatrixException(getRowDimension(), getColumnDimension());
+        }
+        
+        if (p == 0) {
+            return MatrixUtils.createRealIdentityMatrix(this.getRowDimension());
+        }
+        
+        if (p == 1) {
+            return this.copy();
+        }
+
+        final int power = p - 1;
+        
+        /*
+         * Only log_2(p) operations is used by doing as follows:
+         * 5^214 = 5^128 * 5^64 * 5^16 * 5^4 * 5^2
+         * 
+         * In general, the same approach is used for A^p.
+         */        
+        
+        final char[] binaryRepresentation = Integer.toBinaryString(power).toCharArray();       
+        final ArrayList<Integer> nonZeroPositions = new ArrayList<Integer>();
+        int maxI = -1;
+        
+        for (int i = 0; i < binaryRepresentation.length; ++i) {            
+            if (binaryRepresentation[i] == '1') {
+                final int pos = binaryRepresentation.length - i - 1;
+                nonZeroPositions.add(pos);
+                
+                // The positions are taken in turn, so maxI is only changed once
+                if (maxI == -1) {
+                    maxI = pos;
+                }
+            }
+        }
+        
+        RealMatrix[] results = new RealMatrix[maxI + 1];
+        results[0] = this.copy();
+        
+        for (int i = 1; i <= maxI; ++i) {
+            results[i] = results[i-1].multiply(results[i-1]);
+        }
+        
+        RealMatrix result = this.copy();
+        
+        for (Integer i : nonZeroPositions) {
+            result = result.multiply(results[i]);
+        }        
+        
+        return result;
     }
 
     /** {@inheritDoc} */

src/main/java/org/apache/commons/math/linear/RealMatrix.java

@@ -98,6 +98,16 @@ public interface RealMatrix extends AnyMatrix {
      *             if rowDimension(this) != columnDimension(m)
      */
     RealMatrix preMultiply(RealMatrix m);
+    
+    /**
+     * Returns the result multiplying this with itself <code>p</code> times. 
+     * Depending on the underlying storage, instability for high powers might occur. 
+     * @param      p raise this to power p
+     * @return     this^p
+     * @throws     IllegalArgumentException if p < 0
+     *             NonSquareMatrixException if the matrix is not square
+     */
+    RealMatrix power(final int p);
 
     /**
      * Returns matrix entries as a two-dimensional array.

src/test/java/org/apache/commons/math/linear/Array2DRowRealMatrixTest.java

@@ -96,6 +96,7 @@ public final class Array2DRowRealMatrixTest {
     // tolerances
     protected double entryTolerance = 10E-16;
     protected double normTolerance = 10E-14;
+    protected double powerTolerance = 10E-16;
 
     /** test dimensions */
     @Test
@@ -222,6 +223,59 @@ public final class Array2DRowRealMatrixTest {
        TestUtils.assertEquals("m3*m4=m5", m3.multiply(m4), m5, entryTolerance);
    }
 
+    public void testPower() {
+        Array2DRowRealMatrix m = new Array2DRowRealMatrix(testData);
+        Array2DRowRealMatrix mInv = new Array2DRowRealMatrix(testDataInv);
+        Array2DRowRealMatrix mPlusInv = new Array2DRowRealMatrix(testDataPlusInv);
+        Array2DRowRealMatrix identity = new Array2DRowRealMatrix(id);
+        
+        TestUtils.assertEquals("m^0", m.power(0),
+            identity, entryTolerance);        
+        TestUtils.assertEquals("mInv^0", mInv.power(0),
+                identity, entryTolerance);        
+        TestUtils.assertEquals("mPlusInv^0", mPlusInv.power(0),
+                identity, entryTolerance);
+        
+        TestUtils.assertEquals("m^1", m.power(1),
+                m, entryTolerance);        
+        TestUtils.assertEquals("mInv^1", mInv.power(1),
+                mInv, entryTolerance);        
+        TestUtils.assertEquals("mPlusInv^1", mPlusInv.power(1),
+                mPlusInv, entryTolerance); 
+        
+        RealMatrix C1 = m.copy();
+        RealMatrix C2 = mInv.copy();
+        RealMatrix C3 = mPlusInv.copy();
+        
+        for (int i = 2; i <= 10; ++i) {
+            C1 = C1.multiply(m);
+            C2 = C2.multiply(mInv);
+            C3 = C3.multiply(mPlusInv);
+            
+            TestUtils.assertEquals("m^" + i, m.power(i),
+                    C1, entryTolerance);        
+            TestUtils.assertEquals("mInv^" + i, mInv.power(i),
+                    C2, entryTolerance);        
+            TestUtils.assertEquals("mPlusInv^" + i, mPlusInv.power(i),
+                    C3, entryTolerance);            
+        }
+        
+        try {
+            Array2DRowRealMatrix mNotSquare = new Array2DRowRealMatrix(testData2T);
+            mNotSquare.power(2);
+            Assert.fail("Expecting NonSquareMatrixException");
+        } catch (NonSquareMatrixException ex) {
+            // ignored
+        }
+        
+        try {
+            m.power(-1);
+            Assert.fail("Expecting IllegalArgumentException");
+        } catch (IllegalArgumentException ex) {
+            // ignored
+        }
+    }
+    
     /** test trace */
     @Test
     public void testTrace() {