As with the last problem that I solved, I sought out a problem that could be solved using linear algebra.
This problem is another random walks problem. This time there are two objects moving through the system and there are absorbing states, locations at which the object cannot move out of.
The first challenge I needed to overcome was how to represent the two objects moving through a single system. In the end, I settled on each entry in the transition matrix representing the single state where object one is in location $a$ and object two is in location $b$.
I did some research on Markov chains with absorbing states. The Wikipedia article on Absorbing Markov chains was very helpful. It described how to find the fundamental matrix $N$ by first finding matrix $Q$. The latter is the normal transition matrix for the system, minus the rows and columns representing absorbing states. Therefore, the matrix $Q^n$ gives the probability of transitioning from one non-absorbing state to another after $n$ turns. Our fundamental matrix $N$ is the sum of all powers of $Q$.
\[N=\displaystyle\sum_{k=0}^\infty Q^k\]Notice how similar this sum is to the geometric sequence
\[n=\displaystyle\sum_{k=0}^\infty q^k=\frac{1}{1-q}=(1-q)^{-1}\]It follows that $N=(I-Q)^{-1}$. The final solution is the row sum of the row that represents the starting state where the dice are on opposite sides of the table.
I wish my program was faster. Currently runs in about 20 seconds with Python
v3.9.6. Most of the time is taken up in the numpy.array allocation on line
8.
Coming back to this problem now that I have more numpy experience, I see that
I could avoid the large allocation on line 8 by constructing the matrix using
numpy.zeros. I am now using a Macbook Air with an M1 chip which reduces the
running time to about 18 seconds with Python 3.9.10. The update to use
numpy.zeros further reduces the time to about 14 seconds.
import numpy
def q(n):
def place(to_i, to_j, val):
to_i, to_j = to_i%n, to_j%n
row[to_i*n+to_j] += val
P = numpy.array([[0] * n**2 for _ in range(n**2)])
for i in range(n):
for j in range(n):
if i == j:
continue
row = P[i*n+j]
place(i, j, 16)
place(i+1, j, 4)
place(i-1, j, 4)
place(i, j+1, 4)
place(i, j-1, 4)
place(i+1, j+1, 1)
place(i+1, j-1, 1)
place(i-1, j+1, 1)
place(i-1, j-1, 1)
mask = [i*n+i for i in range(n)]
P = numpy.delete(P, mask, 0)
P = numpy.delete(P, mask, 1)
return P / 36
def solve(n):
I = numpy.identity(n**2-n)
Q = q(n)
N = numpy.linalg.inv(I - Q)
starting_row = N[n//2-1]
expected = sum(starting_row)
return round(expected, 6)
if __name__ == '__main__':
import sys
n = eval(sys.argv[1])
print(solve(n))
import pytest
from problem import solve, q
_test_q = (
(2, [[16, 4],
[4, 16]]),
(3, [[16, 4, 1, 1, 1, 4],
[4, 16, 1, 4, 1, 1],
[1, 1, 16, 4, 4, 1],
[1, 4, 4, 16, 1, 1],
[1, 1, 4, 1, 16, 4],
[4, 1, 1, 1, 4 ,16]]),
(4, [[16, 4, 0, 1, 1, 0, 0, 0, 0, 1, 4, 1],
[4, 16, 4, 0, 4, 1, 0, 0, 0, 0, 1, 4],
[0, 4, 16, 1, 1, 4, 0, 0, 0, 1, 0, 1],
[1, 0, 1, 16, 0, 4, 4, 1, 1, 0, 0, 0],
[1, 4, 1, 0, 16, 4, 0, 1, 1, 0, 0, 0],
[0, 1, 4, 4, 4, 16, 1, 0, 4, 0, 0, 0],
[0, 0, 0, 4, 0, 1, 16, 4, 4, 4, 1, 0],
[0, 0, 0, 1, 1, 0, 4, 16, 0, 1, 4, 1],
[0, 0, 0, 1, 1, 4, 4, 0, 16, 1, 0, 1],
[1, 0, 1, 0, 0, 0, 4, 1, 1, 16, 4, 0],
[4, 1, 0, 0, 0, 0, 1, 4, 0, 4, 16, 4],
[1, 4, 1, 0, 0, 0, 0, 1, 1, 0, 4, 16]]),
)
@pytest.mark.parametrize('n,expect', _test_q)
def test_q(n, expect):
assert expect == (q(n) * 36).tolist()
_test_solve = (
(2, 2.25),
(4, 7.2),
(6, 15.340909),
(8, 26.44898),
(10, 40.561927),
(12, 57.674227),
(14, 77.786608),
(16, 100.898979),
(18, 127.011352),
(20, 156.123724),
(50, 952.809311),
(55, 1150.840242),
(60, 1368.371173),
)
@pytest.mark.parametrize('n,expect', _test_solve)
def test_solve(n, expect):
assert expect == solve(n)