Numpy matrix power/exponent with modulo？

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

显而易见的方法有什么问题？

E.g.

import numpy as np 

x = np.arange(100).reshape(10,10) 
y = np.linalg.matrix_power(x, 2) % 50 

Answer 3

从numpy的使用实现：

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

我加入了模长期适应它。 但是，有一个错误，如果发生溢出，没有OverflowError或任何其他类型的异常引发。从那时起，解决方案将是错误的。有一个错误报告here。

这是代码。请谨慎使用：

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