Numpy矩阵幂/指数与模数？

Question

是否有可能使用numpy的linalg.matrix_power与模数,因此元素不会增长大于某个值？

Answer 1

为了防止溢出,如果您首先获取每个输入数字的模数,则可以使用相同结果的事实; 事实上:

(M**k) mod p = ([M mod p]**k) mod p,

对于矩阵 M.这来自以下两个基本身份,它们对整数有效,x并且y:

(x+y) mod p = ([x mod p]+[y mod p]) mod p  # All additions can be done on numbers *modulo p*
(x*y) mod p = ([x mod p]*[y mod p]) mod p  # All multiplications can be done on numbers *modulo p*

同样的标识也适用于矩阵,因为矩阵加法和乘法可以通过标量加法和乘法来表示.有了这个,你只对小数进行取幂(n mod p通常远小于n),并且不太可能出现溢出.在NumPy中,你只需这样做

((arr % p)**k) % p

为了得到(arr**k) mod p.

如果这仍然不够(即,如果存在[n mod p]**k导致溢出的风险,尽管n mod p很小),您可以将取幂分解为多个取幂.产生的基本身份

(n**[a+b]) mod p = ([{n mod p}**a mod p] * [{n mod p}**b mod p]) mod p

和

(n**[a*b]) mod p = ([n mod p]**a mod p)**b mod p.

因此,您可以将功率分解k为a+b+…或作为a*b*…其任何组合.上述标识允许您仅通过小数字执行小数的取幂,这大大降低了整数溢出的风险.

Answer 2

使用 Numpy 的实现：

https://github.com/numpy/numpy/blob/master/numpy/matrixlib/defmatrix.py#L98

我通过添加模项来对其进行调整。然而，存在一个错误，即如果发生溢出，则不会OverflowError引发或引发任何其他类型的异常。从那时起，解决方案将是错误的。这里有一个错误报告。

这是代码。小心使用：

from numpy.core.numeric import concatenate, isscalar, binary_repr, identity, asanyarray, dot
from numpy.core.numerictypes import issubdtype    
def matrix_power(M, n, mod_val):
    # Implementation shadows numpy's matrix_power, but with modulo included
    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 % mod_val
    if n <= 3:
        for _ in range(n-1):
            result = dot(result, M) % mod_val
        return result

    # binary decompositon 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 = dot(Z, Z) % mod_val
        q += 1
    result = Z
    for k in range(q+1, t):
        Z = dot(Z, Z) % mod_val
        if beta[t-k-1] == '1':
            result = dot(result, Z) % mod_val
    return result % mod_val