モジュロを持つナンビ行列パワー/指数？

Question

numpyのlinalg.matrix_powerをモジュロで使用して、要素が特定の値よりも大きくならないようにすることはできますか？

Answer 1

オーバーフローを防止するために、あなたが最初にあなたの入力番号のそれぞれの剰余を取る場合は、同じ結果を得るという事実を使用することができます。実際には：行列Mため

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

。行列の加算と乗算はスカラー加算と乗算で表現することができますので、

(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* 

同じアイデンティティ、同様の行列のために保持する：これは、整数xとyのための有効な次の二つの基本的なアイデンティティは、から来ています。これにより、小さな数値だけが指数化され（n mod pは一般にnよりもずっと小さくなります）、オーバーフローを起こす可能性は非常に低くなります。 NumPyでは、を得るには、

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

を取得するだけです。

まだ十分でない場合（n mod pが小さいにもかかわらず[n mod p]**kがオーバーフローするリスクがある場合）、累乗を複数の累乗に分割することができます。収量

(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. 

上記の基本的なアイデンティティがこのように、あなたはa+b+…またはa*b*…またはそれらの任意の組み合わせとしてパワーkを破ることができます。上記のアイデンティティを使用すると、小さい数値の累乗だけを小数で実行できるため、整数のオーバーフローのリスクが大幅に低下します。 numpyのの実装を使用して

Answer 2

明白なアプローチで何が問題になっていますか？

など。

import numpy as np 

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

Answer 3

：

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