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