Матрица с коэффициентом мощности / экспонента с модулем? - PullRequest
Новая
       37

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

10 голосов
/ 15 декабря 2011

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

Ответы [ 3 ]

10 голосов
/ 15 декабря 2011

Чтобы предотвратить переполнение, вы можете использовать тот факт, что вы получите тот же результат, если вы сначала берете по модулю каждого из ваших входных чисел; на самом деле: 

(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*… или любую их комбинацию. Приведенные выше идентификаторы позволяют выполнять возведения только небольших чисел в маленькие числа, что значительно снижает риск переполнения целых чисел.

4 голосов
/ 09 октября 2014

Использование реализации от Numpy:

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

Я адаптировал его, добавив по модулю термин. ОДНАКО , есть ошибка, в которой при переполнении не возникает OverflowError или любое другое исключение. С этого момента решение будет неправильным. Отчет об ошибке здесь .

Вот код. Используйте с осторожностью: 

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
1 голос
/ 15 декабря 2011

Что не так с очевидным подходом?

Например 

import numpy as np

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