没有 Numpy 的 3x3 矩阵的矩阵功率

Question

我试图在不使用 Numpy 的 3x3 矩阵（不使用任何库函数）的情况下计算 Matrix 的 n 次幂

这是我目前写的代码：

   def matmul(M1, M2):
        a1 = (M1[0][0]*M2[0][0]) + (M1[0][1] * M2[1][0]) + (M1[0][2] * M2[2][0])
        a2 = (M1[0][0]*M2[0][1]) + (M1[0][1] * M2[1][1]) + (M1[0][2] * M2[2][1])
        a3 = (M1[0][0]*M2[0][2]) + (M1[0][1] * M2[1][2]) + (M1[0][2] * M2[2][2])
        a4 = (M1[1][0]*M2[0][0]) + (M1[1][1] * M2[1][0]) + (M1[1][2] * M2[2][0])
        a5 = (M1[1][0]*M2[0][1]) + (M1[1][1] * M2[1][1]) + (M1[1][2] * M2[2][1])
        a6 = (M1[1][0]*M2[0][2]) + (M1[1][1] * M2[1][2]) + (M1[1][2] * M2[2][2])
        a7 = (M1[2][0]*M2[0][0]) + (M1[2][1] * M2[1][0]) + (M1[2][2] * M2[2][0])
        a8 = (M1[2][0]*M2[0][1]) + (M1[2][1] * M2[1][1]) + (M1[2][2] * M2[2][1])
        a9 = (M1[2][0]*M2[0][2]) + (M1[2][1] * M2[1][2]) + (M1[2][2] * M2[2][2])
        return [[a1, a2, a3], [a4, a5, a6], [a7, a8, a9]]
    def matpow(mat, p):
        if p == 1:
            return mat
        m2 = matpow(mat, p//2)
        if p%2 == 0:
            return matmul(m2, m2)
        else:
            return matmul(matmul(m2, m2), mat)

这适用于较小数量的“n”，但在应用于较大数量时会失败并出现以下错误：

 File "/workspace/default/solution.py", line 17, in matpow
    m2 = matpow(mat, p//2)
  File "/workspace/default/solution.py", line 17, in matpow
    m2 = matpow(mat, p//2)
  [Previous line repeated 990 more times]
  File "/workspace/default/solution.py", line 15, in matpow
    if p == 1:
RecursionError: maximum recursion depth exceeded in comparison

似乎在尝试大数时达到了递归堆栈深度，有没有一种方法可以在不使用库的情况下以优化的方式实现功能？

Answer 1

如果幂参数确实巨大，则使用迭代解决方案，通过迭代幂的二进制位：

def matpow(mat, p):
    # Start with the identity matrix. This way it will also work for p==0
    # If you are using floats, then make this also with floats:
    result = [[1,0,0],[0,1,0],[0,0,1]]
    while p > 0:
        if p & 1:
            result = matmul(result, mat)
        mat = matmul(mat, mat)
        p >>= 1
    return result

Answer 2

使用递归时，堆在内存中新建一条记录来存放返回寄存器，这样每次调用函数matpow时，程序就知道条件满足时返回哪个寄存器。因此，我不使用递归，而是使用迭代。

使用 while/for 循环进行迭代，直到找到正确的值。

这应该可以工作（经过测试），

def matpow_iteratively(mat, p):
    new_mat = mat
    arr = []

    while p > 1:
        arr.append(p)
        p = p // 2

    index = len(arr) - 1
    while index >= 0:
        if arr[index] % 2 == 0:
            new_mat = matmul(new_mat, new_mat)
        else:
            new_mat = matmul(matmul(new_mat, new_mat), mat)
        index = index - 1

    return new_mat