Матрица мощности/экспоненты с модулем?

Question

Возможно ли использовать linalg.matrix_power numpy с модулем, чтобы элементы не становились больше определенного значения?

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

Я приспособил его добавлением члена по модулю. HOWEVER, есть ошибка, при которой при переполнении возникает 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