Numpy matrix power/exponent with modulo?

Question

Czy można używać numinalu linalg.matrix_power z modulo, więc elementy nie rosną większe niż pewna wartość?

Answer 1

W celu uniknięcia przepełnienia, można wykorzystać fakt, że masz ten sam wynik, jeśli najpierw podjąć modulo każdego z numerów wejściowych; w rzeczywistości:

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

na matrycyM. Wynika to z następujących dwóch zasadniczych identyczności, które są ważne dla liczb x i 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* 

same identyfikatory przytrzymać matryc, a także, ponieważ dodawanie mnożenie macierzy i może być wyrażona przez skalarne dodawania i mnożenia. Dzięki temu tylko wykładniczo małe liczby (n mod p jest zwykle znacznie mniejsza niż n) i są znacznie mniej prawdopodobne, aby uzyskać przepełnienia. W NumPy, byś więc po prostu zrobić

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

w celu uzyskania (arr**k) mod p.

Jeśli to nadal za mało (tzn. Jeśli istnieje ryzyko, że [n mod p]**k spowoduje przepełnienie, mimo że n mod p jest małe), można podzielić potęgowanie na wiele potęgowań. Podstawowe tożsamości powyżej wydajnością

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

i

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

sposób można zerwać moc k jak a+b+… lub a*b*… lub ich dowolnej kombinacji. Powyższe tożsamości umożliwiają wykonywanie tylko potęgowań małych liczb przez małe liczby, co znacznie obniża ryzyko przekroczenia liczb całkowitych.

Answer 2

Co jest nie tak z oczywistym podejściem?

E.g.

import numpy as np 

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

Answer 3

Korzystanie realizację od NumPy:

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

I dostosowany dodając termin modulo. JEDNAK, jest błąd, ponieważ jeśli wystąpi przepełnienie, nie zostanie zgłoszony żaden inny wyjątek. Od tego momentu rozwiązanie będzie błędne. Istnieje zgłoszenie błędu here.

Oto kod. Używaj ostrożnie:

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