如何在NumPy中进行模运算？矩阵幂取模的高效实现咨询

我需要执行如下计算：

import numpy as np
mat = np.array([
    [0, 0, 0, 1],
    [1, 0, 1, 0],
    [1, 1, 0, 0],
    [1, 1, 1, 1],
])
ret = np.linalg.matrix_power(mat, 1000)

计算结果会包含超出64位范围的超大数值，但我仅需对结果取10**9+7的模。目前我可以将矩阵幂展开为多次矩阵乘法，每次乘法后手动取模，示例代码如下：

import numpy as np
mod = 10**9 +7
ret = np.identity(4)
for _ in range(1000):
    ret @= mat
    ret %= mod

但这种方式会失去NumPy矩阵幂运算的优化特性，请问NumPy是否提供可直接指定结果取模的API？

NumPy原生的np.linalg.matrix_power并没有直接支持取模的参数，没办法在计算矩阵幂的过程中自动对结果取模。

不过你可以自己实现**快速幂（二分幂）**的版本，既能兼顾效率，又能满足取模需求——这个思路和NumPy矩阵幂的优化逻辑一致，都是通过二分法减少乘法次数，把时间复杂度从O(n)降到O(log n)：

import numpy as np

def matrix_power_mod(mat, power, mod):
    # 初始化结果为同维度单位矩阵，用int64避免中间溢出
    result = np.identity(mat.shape[0], dtype=np.int64)
    current_mat = mat.copy().astype(np.int64)
    while power > 0:
        # 若当前幂次为奇数，将当前矩阵乘入结果
        if power % 2 == 1:
            result = (result @ current_mat) % mod
        # 矩阵自乘，幂次折半
        current_mat = (current_mat @ current_mat) % mod
        power = power // 2
    return result

# 调用示例
mod = 10**9 + 7
mat = np.array([
    [0, 0, 0, 1],
    [1, 0, 1, 0],
    [1, 1, 0, 0],
    [1, 1, 1, 1],
])
ret = matrix_power_mod(mat, 1000, mod)

每次乘法后立即取模，能有效避免数值超出64位整数范围的问题，同时二分法的设计保证了运算效率接近NumPy原生的矩阵幂实现。

内容的提问来源于stack exchange，提问作者zzkluck