numpy: повышение матрицы до мощности дает странный результат

Question

Мое решение проблемы KnightDialer на Leetcode :

import numpy as np
class Solution:
  def knightDialer(self, N: int) -> int:
    M = np.matrix([[0, 0, 0, 0, 1, 0, 1, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
                   [0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
                   [0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
                   [1, 0, 0, 1, 0, 0, 0, 0, 0, 1],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [1, 1, 0, 0, 0, 0, 0, 1, 0, 0],
                   [0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
                   [0, 1, 0, 1, 0, 0, 0, 0, 0, 0],
                   [0, 0, 1, 0, 1, 0, 0, 0, 0, 0]])
    return np.sum(M**(N-1)) % (10**9 + 7)

Он работает правильно для значений от N до 51. Например, для N = 1 он правильно возвращает 10, для N = 2 корректно возвращает 20, для N = 3 корректно возвращает 46 и т. Д. Чем при N> 51 он перестает давать точный результат (для N = 52 он возвращает 107679695 вместо 690023703). Я не уверен, почему, но как-то повышение матрицы до степени> 51 приводит к неточному результату.

Я пытался заменить M**(N-1) на np.linalg.matrix_power(M, (N-1)), но вывод все еще не точен. Я догадываюсь, что за кулисами идут какие-то numpy "маги" c, но я не уверен, что это такое.

Answer 1

Numpy борется, потому что он работает с целыми числами фиксированного размера, такими как int32 или int64 (какое из них используется здесь, зависит от вашей установки python). Такое экспонирование матрицы быстро делает записи больше предела этого типа данных, что приводит к усечению.

С обычными Python целыми числами можно использовать этот вид реализации (сначала умножьте матрицу и применять модульное сокращение впоследствии), поскольку обычные целые числа Python могут иметь гораздо более высокие значения. Например:

def mm(A, B):
    return [[sum([x*y for (x, y) in zip(row, col)]) for col in zip(*B)] for row in A]

def knightDialer(N: int) -> int:
    M = [[0, 0, 0, 0, 1, 0, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
        [0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
        [0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
        [1, 0, 0, 1, 0, 0, 0, 0, 0, 1],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
        [1, 1, 0, 0, 0, 0, 0, 1, 0, 0],
        [0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
        [0, 1, 0, 1, 0, 0, 0, 0, 0, 0],
        [0, 0, 1, 0, 1, 0, 0, 0, 0, 0]]
    N = N - 1
    res = [[int(i == j) for j in range(len(M))] for i in range(len(M))]
    while N:
        if N & 1: res = mm(res, M)
        M = mm(M, M)
        N >>= 1
    print(M)
    return sum([sum(i) for i in zip(*res)]) % (10**9 + 7)

Применение модульного сокращения во время возведения в степень позволяет numpy матрицам работать без битов:

def knightDialer(N: int) -> int:
    M = np.matrix([[0, 0, 0, 0, 1, 0, 1, 0, 0, 0],
                   [0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
                   [0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
                   [0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
                   [1, 0, 0, 1, 0, 0, 0, 0, 0, 1],
                   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                   [1, 1, 0, 0, 0, 0, 0, 1, 0, 0],
                   [0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
                   [0, 1, 0, 1, 0, 0, 0, 0, 0, 0],
                   [0, 0, 1, 0, 1, 0, 0, 0, 0, 0]], dtype=np.int64)
    N = N - 1
    res = np.eye(M.shape[0], dtype=np.int64)
    while N:
        if N & 1: res = res * M % (10**9 + 7)
        M = M * M % (10**9 + 7)
        N >>= 1
    return np.sum(res) % (10**9 + 7)

Использование dtype=np.int64 необходимо в случае, если ваш установка целочисленный тип по умолчанию int32. int32 достаточно велик, чтобы содержать числа ниже 10**9 + 7, но во время умножения матриц между двумя такими числами (а затем они также суммируются) будут произведения, и в это время может произойти переполнение, если использовалось int32.