def matrix_power(M, n):
M = asanyarray(M)<EOL>if len(M.shape) != <NUM_LIT:2> or M.shape[<NUM_LIT:0>] != M.shape[<NUM_LIT:1>]:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>if not issubdtype(type(n), int):<EOL><INDENT>raise TypeError("<STR_LIT>")<EOL><DEDENT>from numpy.linalg import inv<EOL>if n==<NUM_LIT:0>:<EOL><INDENT>M = M.copy()<EOL>M[:] = identity(M.shape[<NUM_LIT:0>])<EOL>return M<EOL><DEDENT>elif n<<NUM_LIT:0>:<EOL><INDENT>M = inv(M)<EOL>n *= -<NUM_LIT:1><EOL><DEDENT>result = M<EOL>if n <= <NUM_LIT:3>:<EOL><INDENT>for _ in range(n-<NUM_LIT:1>):<EOL><INDENT>result=N.dot(result, M)<EOL><DEDENT>return result<EOL><DEDENT>beta = binary_repr(n)<EOL>Z, q, t = M, <NUM_LIT:0>, len(beta)<EOL>while beta[t-q-<NUM_LIT:1>] == '<STR_LIT:0>':<EOL><INDENT>Z = N.dot(Z, Z)<EOL>q += <NUM_LIT:1><EOL><DEDENT>result = Z<EOL>for k in range(q+<NUM_LIT:1>, t):<EOL><INDENT>Z = N.dot(Z, Z)<EOL>if beta[t-k-<NUM_LIT:1>] == '<STR_LIT:1>':<EOL><INDENT>result = N.dot(result, Z)<EOL><DEDENT><DEDENT>return result<EOL>
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 quarternion 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.]])