MAINT: Move matrix_power to linalg

The docstring already assumed it was in linalg, and this ensures
linalg becomes completely independent of matrixlib.

1 files changed, 111 insertions, 3 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 5ee230f92..dc1fa5dea 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -24,12 +24,12 @@ from numpy.core import (
     add, multiply, sqrt, maximum, fastCopyAndTranspose, sum, isfinite, size,
     finfo, errstate, geterrobj, longdouble, moveaxis, amin, amax, product, abs,
     broadcast, atleast_2d, intp, asanyarray, object_, ones, matmul,
-    swapaxes, divide, count_nonzero, ndarray, isnan
+    swapaxes, divide, count_nonzero, ndarray, isnan, identity, issubdtype,
+    binary_repr, integer
 )
 from numpy.core.multiarray import normalize_axis_index
-from numpy.lib import triu, asfarray
+from numpy.lib.twodim_base import triu
 from numpy.linalg import lapack_lite, _umath_linalg
-from numpy.matrixlib.defmatrix import matrix_power
 
 # For Python2/3 compatibility
 _N = b'N'
@@ -532,6 +532,114 @@ def inv(a):
     return wrap(ainv.astype(result_t, copy=False))
 
 
+def matrix_power(M, n):
+    """
+    Raise a square matrix to the (integer) power `n`.
+
+    For positive integers `n`, the power is computed by repeated matrix
+    squarings and matrix multiplications. If ``n == 0``, the identity matrix
+    of the same shape as M is returned. If ``n < 0``, the inverse
+    is computed and then raised to the ``abs(n)``.
+
+    Parameters
+    ----------
+    M : ndarray or matrix object
+        Matrix to be "powered."  Must be square, i.e. ``M.shape == (m, m)``,
+        with `m` a positive integer.
+    n : int
+        The exponent can be any integer or long integer, positive,
+        negative, or zero.
+
+    Returns
+    -------
+    M**n : ndarray or matrix object
+        The return value is the same shape and type as `M`;
+        if the exponent is positive or zero then the type of the
+        elements is the same as those of `M`. If the exponent is
+        negative the elements are floating-point.
+
+    Raises
+    ------
+    LinAlgError
+        If the matrix is not numerically invertible.
+
+    See Also
+    --------
+    matrix
+        Provides an equivalent function as the exponentiation operator
+        (``**``, not ``^``).
+
+    Examples
+    --------
+    >>> from numpy import linalg as LA
+    >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
+    >>> LA.matrix_power(i, 3) # should = -i
+    array([[ 0, -1],
+           [ 1,  0]])
+    >>> LA.matrix_power(np.matrix(i), 3) # matrix arg returns matrix
+    matrix([[ 0, -1],
+            [ 1,  0]])
+    >>> LA.matrix_power(i, 0)
+    array([[1, 0],
+           [0, 1]])
+    >>> LA.matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
+    array([[ 0.,  1.],
+           [-1.,  0.]])
+
+    Somewhat more sophisticated example
+
+    >>> q = np.zeros((4, 4))
+    >>> q[0:2, 0:2] = -i
+    >>> q[2:4, 2:4] = i
+    >>> q # one of the three quaternion units not equal to 1
+    array([[ 0., -1.,  0.,  0.],
+           [ 1.,  0.,  0.,  0.],
+           [ 0.,  0.,  0.,  1.],
+           [ 0.,  0., -1.,  0.]])
+    >>> LA.matrix_power(q, 2) # = -np.eye(4)
+    array([[-1.,  0.,  0.,  0.],
+           [ 0., -1.,  0.,  0.],
+           [ 0.,  0., -1.,  0.],
+           [ 0.,  0.,  0., -1.]])
+
+    """
+    M = asanyarray(M)
+    if M.ndim != 2 or M.shape[0] != M.shape[1]:
+        raise ValueError("input must be a square array")
+    if not issubdtype(type(n), integer):
+        raise TypeError("exponent must be an integer")
+
+    from numpy.linalg import inv
+
+    if n==0:
+        M = M.copy()
+        M[:] = identity(M.shape[0])
+        return M
+    elif n<0:
+        M = inv(M)
+        n *= -1
+
+    result = M
+    if n <= 3:
+        for _ in range(n-1):
+            result=dot(result, M)
+        return result
+
+    # binary decomposition to reduce the number of Matrix
+    # multiplications for n > 3.
+    beta = binary_repr(n)
+    Z, q, t = M, 0, len(beta)
+    while beta[t-q-1] == '0':
+        Z = dot(Z, Z)
+        q += 1
+    result = Z
+    for k in range(q+1, t):
+        Z = dot(Z, Z)
+        if beta[t-k-1] == '1':
+            result = dot(result, Z)
+    return result
+
+
 # Cholesky decomposition
 
 def cholesky(a):