def matrix_power(a, n)
    a = NArray.asarray(a)
    m,k = a.shape[-2..-1]
    unless m==k
      raise NArray::ShapeError, "input must be a square array"
    end
    unless Integer===n
      raise ArgumentError, "exponent must be an integer"
    end
    if n == 0
      return a.class.eye(m)
    elsif n < 0
      a = inv(a)
      n = n.abs
    end
    if n <= 3
      r = a
      (n-1).times do
        r = matmul(r,a)
      end
    else
      while (n & 1) == 0
        a = matmul(a,a)
        n >>= 1
      end
      r = a
      while n != 0
        a = matmul(a,a)
        n >>= 1
        if (n & 1) != 0
          r = matmul(r,a)
        end
      end
    end
    r
  end
def matrix_power(a, n)
    a = NArray.asarray(a)
    m,k = a.shape[-2..-1]
    unless m==k
      raise NArray::ShapeError, "input must be a square array"
    end
    unless Integer===n
      raise ArgumentError, "exponent must be an integer"
    end
    if n == 0
      return a.class.eye(m)
    elsif n < 0
      a = inv(a)
      n = n.abs
    end
    if n <= 3
      r = a
      (n-1).times do
        r = matmul(r,a)
      end
    else
      while (n & 1) == 0
        a = matmul(a,a)
        n >>= 1
      end
      r = a
      while n != 0
        a = matmul(a,a)
        n >>= 1
        if (n & 1) != 0
          r = matmul(r,a)
        end
      end
    end
    r
  end
Compute a square matrix `a` to the power `n`.

  * If n > 0: return `a**n`.
  * If n == 0: return identity matrix.
  * If n < 0: return `(a*\*-1)*\*n.abs`.

 @param a [Numo::NArray] square matrix (>= 2-dimensinal NArray).
 @param n [Integer] the exponent.
 @example
     i = Numo::DFloat[[0, 1], [-1, 0]]
     => Numo::DFloat#shape=[2,2]
     [[0, 1],
      [-1, 0]]
     Numo::Linalg.matrix_power(i,3)
     => Numo::DFloat#shape=[2,2]
     [[0, -1],
      [1, 0]]
     Numo::Linalg.matrix_power(i,0)
     => Numo::DFloat#shape=[2,2]
     [[1, 0],
      [0, 1]]
     Numo::Linalg.matrix_power(i,-3)
     => Numo::DFloat#shape=[2,2]
     [[0, 1],
      [-1, 0]]

     q = Numo::DFloat.zeros(4,4)
     q[0..1,0..1] = -i
     q[2..3,2..3] = i
     q
     => Numo::DFloat#shape=[4,4]
     [[-0, -1, 0, 0],
      [1, -0, 0, 0],
      [0, 0, 0, 1],
      [0, 0, -1, 0]]
     Numo::Linalg.matrix_power(q,2)
     => Numo::DFloat#shape=[4,4]
     [[-1, 0, 0, 0],
      [0, -1, 0, 0],
      [0, 0, -1, 0],
      [0, 0, 0, -1]]