def ix_(*args):
    """
    Construct an open mesh from multiple sequences.

    This function takes N 1-D sequences and returns N outputs with N
    dimensions each, such that the shape is 1 in all but one dimension
    and the dimension with the non-unit shape value cycles through all
    N dimensions.

    Using `ix_` one can quickly construct index arrays that will index
    the cross product. ``a[np.ix_([1,3],[2,5])]`` returns the array
    ``[[a[1,2] a[1,5]], [a[3,2] a[3,5]]]``.

    Parameters
    ----------
    args : 1-D sequences

    Returns
    -------
    out : tuple of ndarrays
        N arrays with N dimensions each, with N the number of input
        sequences. Together these arrays form an open mesh.

    See Also
    --------
    ogrid, mgrid, meshgrid

    Examples
    --------
    >>> a = np.arange(10).reshape(2, 5)
    >>> a
    array([[0, 1, 2, 3, 4],
           [5, 6, 7, 8, 9]])
    >>> ixgrid = np.ix_([0,1], [2,4])
    >>> ixgrid
    (array([[0],
           [1]]), array([[2, 4]]))
    >>> ixgrid[0].shape, ixgrid[1].shape
    ((2, 1), (1, 2))
    >>> a[ixgrid]
    array([[2, 4],
           [7, 9]])

    """
    out = []
    nd = len(args)
    for k, new in enumerate(args):
        new = asarray(new)
        if new.ndim != 1:
            raise ValueError("Cross index must be 1 dimensional")
        if new.size == 0:
            # Explicitly type empty arrays to avoid float default
            new = new.astype(_nx.intp)
        if issubdtype(new.dtype, _nx.bool_):
            new, = new.nonzero()
        new = new.reshape((1,)*k + (new.size,) + (1,)*(nd-k-1))
        out.append(new)
    return tuple(out)

 def compare(x, y):
     try:
         if npany(gisinf(x)) or npany( gisinf(y)):
             xinfid = gisinf(x)
             yinfid = gisinf(y)
             if not xinfid == yinfid:
                 return False
             # if one item, x and y is +- inf
             if x.size == y.size == 1:
                 return x == y
             x = x[~xinfid]
             y = y[~yinfid]
     except TypeError:
         pass
     z = abs(x-y)
     if not issubdtype(z.dtype, number):
         z = z.astype(float_) # handle object arrays
     return around(z, decimal) <= 10.0**(-decimal)

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 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.]])

    """
    M = asanyarray(M)
    if len(M.shape) != 2 or M.shape[0] != M.shape[1]:
        raise ValueError("input must be a square array")
    if not issubdtype(type(n), int):
        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=N.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 = N.dot(Z, Z)
        q += 1
#.........这里部分代码省略.........

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 type as M is returned. If n<0, the inverse is computed
    and raised to the exponent.

    Parameters
    ----------
    M : array_like
        Must be a square array (that is, of dimension two and with
        equal sizes).
    n : integer
        The exponent can be any integer or long integer, positive
        negative or zero.

    Returns
    -------
    M to the power n
        The return value is a an array the same shape and size as M;
        if the exponent was positive or zero then the type of the
        elements is the same as those of M. If the exponent was negative
        the elements are floating-point.

    Raises
    ------
    LinAlgException
        If the matrix is not numerically invertible, an exception is raised.

    See Also
    --------
    The matrix() class provides an equivalent function as the exponentiation
    operator.

    Examples
    --------
    >>> np.linalg.matrix_power(np.array([[0,1],[-1,0]]),10)
    array([[-1,  0],
           [ 0, -1]])

    """
    M = asanyarray(M)
    if len(M.shape) != 2 or M.shape[0] != M.shape[1]:
        raise ValueError("input must be a square array")
    if not issubdtype(type(n),int):
        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=N.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 = N.dot(Z,Z)
        q += 1
    result = Z
    for k in range(q+1,t):
        Z = N.dot(Z,Z)
        if beta[t-k-1] == '1':
            result = N.dot(result,Z)
    return result