Skip to content

Instantly share code, notes, and snippets.

@thejevans
Last active January 29, 2024 05:54
Show Gist options
  • Save thejevans/fb421e00773f5a0e8372620b6545653f to your computer and use it in GitHub Desktop.
Save thejevans/fb421e00773f5a0e8372620b6545653f to your computer and use it in GitHub Desktop.
euler213
import numpy as np
# build transition matrix
transition_matrix = np.zeros((900,900))
# corner cases
corner_transition_idx = [(0, 1), (0, 30), (29, 28), (29, 59), (870, 871), (870, 840), (899, 898), (899, 869)]
# edge cases
top_idx = [(x, x+30) for x in range(1, 29)] + [(x, x-1) for x in range(1, 29)] + [(x, x+1) for x in range(1, 29)]
left_idx = [(x, x+30) for x in range(30, 870, 30)] + [(x, x-30) for x in range(30, 870, 30)] + [(x, x+1) for x in range(30, 870, 30)]
right_idx = [(x, x+30) for x in range(59, 899, 30)] + [(x, x-30) for x in range(59, 899, 30)] + [(x, x-1) for x in range(59, 899, 30)]
bottom_idx = [(x, x-30) for x in range(871, 899)] + [(x, x-1) for x in range(871, 899)] + [(x, x+1) for x in range(871, 899)]
edge_transition_idx = top_idx + left_idx + right_idx + bottom_idx
# everything else
interior_transition_idx = []
for x in range(31,869):
if x % 30 in (0, 29):
continue
interior_transition_idx += [(x, x+1), (x, x-1), (x, x+30), (x, x-30)]
# combine
for x, y in corner_transition_idx:
transition_matrix[x, y] = 0.5
for x, y in edge_transition_idx:
transition_matrix[x, y] = 1/3
for x, y in interior_transition_idx:
transition_matrix[x, y] = 0.25
# compute A^50
after_fifty = np.linalg.matrix_power(transition_matrix, 50)
# get probability of no flea
no_flea_A = 1 - after_fifty
# multiply across rows
per_cell_no_flea = np.prod(no_flea_A, axis=0)
# sum all to final expectation value
np.sum(per_cell_no_flea)
import numpy as np
# build transition matrix
k1_diag = np.ones(899) * 0.25
k1_diag[np.arange(1, 30) * 30 - 1] = 0
k1_diag[np.arange(1, 30) * 30] = 1/3
k1_diag[[0, 870]] = 0.5
k1_diag[1:29] = 1/3
k1_diag[871:] = 1/3
k30_diag = np.ones(870) * 0.25
k30_diag[np.arange(1, 30) * 30 - 1] = 1/3
k30_diag[np.arange(1, 29) * 30] = 1/3
k30_diag[[0, 29]] = 0.5
k30_diag[1:29] = 1/3
tm = np.diag(k1_diag, k=1) + np.diag(k1_diag[::-1], k=-1) + np.diag(k30_diag, k=30) + np.diag(k30_diag[::-1], k=-30)
# compute expectation value
np.sum(np.prod(1 - np.linalg.matrix_power(tm, 50), axis=0))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment