Python code to toy around with Sandpiles, as described on Numberphile: https://www.youtube.com/watch?v=1MtEUErz7Gg

sandpile

#!/usr/bin/env python

"""Computations on Sandpiles.

See: https://www.youtube.com/watch?v=1MtEUErz7Gg
"""


import numpy as np

class Sandpile:
    nums = np.zeros((3, 3), int)

    def __init__(self, nums: np.ndarray=None):
        if nums is not None:
            self.nums = np.array(nums, int)

    def __repr__(self) -> str:
        return 'Sandpile:\n%s' % self.nums

    def __add__(self, other: 'Sandpile') -> 'Sandpile':
        """Adds two sandpiles together."""

        intermed = Sandpile(self.nums + other.nums)
        intermed.topple()
        return intermed

    def topple(self):
        while True:
            to_topple = self.nums > 3
            if not to_topple.any():
                break

            for rowidx, row in enumerate(to_topple):
                for colidx, should_topple in enumerate(row):
                    if not should_topple: continue
                    self.nums[rowidx][colidx] -= 4
                    for drow, dcol in ((-1, 0), (0, -1), (1, 0), (0, 1)):
                        to_row = rowidx + drow
                        to_col = colidx + dcol
                        if not 0 <= to_row < self.nums.shape[0]: continue
                        if not 0 <= to_col < self.nums.shape[1]: continue

                        self.nums[to_row][to_col] += 1
            print('After toppling: %s' % self)

if __name__ == '__main__':
    # a = Sandpile([[2, 2, 2], [0, 0, 0], [0, 0, 0]])
    # b = Sandpile([[2, 2, 3], [1, 0, 1], [1, 1, 1]])
    # a = Sandpile([[2, 2, 2], [2, 2, 2], [2, 2, 2]])
    # a = Sandpile([[3, 3, 3], [3, 3, 3], [3, 3, 3]])
    # b = Sandpile([[3, 3, 3], [3, 1, 3], [3, 3, 3]])
    # b = Sandpile([[2, 1, 2], [1, 0, 1], [2, 1, 2]])
    a = Sandpile([[2, 2], [2, 2]])
    b = Sandpile([[1, 2], [0, 2]])

    import doctest
    doctest.testmod()
    print('a = %s' % a)
    print('b = %s' % b)
    print('a + b = %s' % (a + b))