CODEWITHSUNDEEP
python-3.xnumpy
Nick
Nick

Reputation: 3134

numpy: raising a matrix to a power produces strange result

My solution to the KnightDialer problem on 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)

It works correctly for values of N up to 51. E.g. for N = 1 it correctly returns 10, for N = 2 it correctly returns 20, for N = 3 it correctly returns 46 and so on. Than at N > 51 it stops producing an accurate result (for N = 52 it returns 107679695 instead of 690023703). I am not sure why but somehow raising the matrix to a power > 51 leads to an inaccurate result.

I tried replacing M**(N-1) with np.linalg.matrix_power(M, (N-1)) but the output is still not accurate. My hunch is that there is some numpy "magic" going behind the scenes but I am not sure what is it.

Upvotes: 1

Views: 210

Answers (1)

user555045
user555045

Reputation: 64904

Numpy struggles because it works with fixed-size integers such as int32 or int64 (which one is used here depends on your python installation). Exponentiating the matrix like this quickly makes the entries larger than the limit of that data type, resulting in truncation.

With regular Python integers, it is possible to use this kind of implementation (multiply the matrix first and apply the modular reduction afterwards), because regular Python integers can have much higher values. For example:

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)

Applying the modular reduction during the exponentiation lets numpy matrixes work without running out of bits:

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)

Using dtype=np.int64 is necessary in case for your installation the default integer type is int32. int32 is big enough to hold numbers below 10**9 + 7, but during the matrix multiplication there will be products between two such numbers (and then they are summed, too) and at that time overflow can occur if int32 was used.

Upvotes: 2

Related Questions