Replicating "Counting the Paths on a Grid" by Donald Knuth in 1976 (An example of importance sampling)

count_path_grid.py

'''
c style python because of the consideration of efficiency, please use pypy3 instead of python3!

This is a replication of "Counting the Paths on a Grid" introduced in the paper "Coping with Finiteness"
(Donald E. Knuth, 1976). This program estimates the number of self-avoiding paths on a 10×10 grid, 
starting from the point(0, 0) and ending at the point(10, 10).

~ pypy3 count_path_grid.py
1.5810562057001226e+24

which is very close to the true number: 1,568,758,030,464,750,013,214,106
'''
import random

N = 11
T = 1000000

grid = [[0 for _ in range(N)] for __ in range(N)]
grid_for_search = [[0 for _ in range(N)] for __ in range(N)]

directs = [[1, 0], [0, 1], [-1, 0], [0, -1]]

def init_grid():
    for i in range(N):
        for j in range(N):
            grid[i][j] = 0

def init_grid_for_search():
    for i in range(N):
        for j in range(N):
            grid_for_search[i][j] = grid[i][j]

def legal(grid, x, y):
    return x >= 0 and x < N and y >= 0 and y < N and grid[x][y] == 0

def flood_fill(x, y):
    if x == N - 1 and y == N - 1:
        return True

    for direct in directs:
        newx = x + direct[0]
        newy = y + direct[1]
        if legal(grid_for_search, newx, newy):
            grid_for_search[newx][newy] = 1
            if flood_fill(newx, newy) == True:
                return True

    return False

def connected(x, y):
    init_grid_for_search()
    grid_for_search[x][y] = 1
    return flood_fill(x, y)

def sampling():
    count = 1
    init_grid()
    cur_pos_x = 0
    cur_pos_y = 0
    grid[cur_pos_x][cur_pos_y] = 1

    while True:
        moves = []
        for direct in directs:
            newx = cur_pos_x + direct[0]
            newy = cur_pos_y + direct[1]
            if legal(grid, newx, newy) and connected(newx, newy):
                moves.append([newx, newy])
        count *= len(moves)
        random.shuffle(moves)
        cur_pos_x = moves[0][0]
        cur_pos_y = moves[0][1]
        grid[cur_pos_x][cur_pos_y] = 1
        if cur_pos_x == N-1 and cur_pos_y == N-1:
            break

    return count

if __name__ == "__main__":
    count = 0
    for _ in range(T):
        count += sampling()
    print(count / T)