GoKawiilPublished 2026-3-8
When I wrote Beautiful Abelian Sandpiles I wanted to show off some nice images of large identity sandpiles. But the simple algorithm I used was horrendously slow. Showing a sandpile identity that was larger than 100 by 100 took multiple seconds! That's not good enough! I became obsessed with trying to find a faster way, all in an effort to compute bigger and bigger sandpile identities, bigger than anything anyone had seen before. In the end, I did exactly that.
Before I discuss anything else, I want to show off this gigantic sandpile identity that is 16384 by 16384. Here it is in all its glory! Feel free to zoom in and explore all the hidden detail. Explore in fullscreen to get the best experience
This is coloured with the matplotlib inferno colourmap, so black, purple, orange, yellow correspond to 0, 1, 2, 3, grains of sand respectively. The sandpile identity was computed on an AMD Ryzen 7 4800H, and it took just under an hour, slightly beating the 10 days required to compute this 10 000 by 10 000 example which, prior to this, was the largest sandpile identity that I knew of. But how do we compute sandpiles quickly? Let's first get an overview of the current algorithms that exist and understand why they work!
Current Methods
There are essentially two different methods for computing an identity. The first is the method that I used in Beautiful Abelian Sandpiles which is given in pseudocode below:
def identity(width, height): sixes = [[6 for _ in range(width)] for _ in range(height) ] sixes_stabilized = stabilize(sixes) difference = [[6 - sixes_stabilized[j][i] for i in range(width)] for j in range(height) ] return stabilize(difference)
I will refer to this as the Difference method. The stabilize function topples all cells until all cells are below 4 . Note that creating a 2D array of 6s and subtracting two 2D arrays is fast. So the time to compute stabilize twice determines the speed of this method.
The second method iteratively builds up the identity by adding sand around the edges in a special configuration called the Burning Configuration. Again, the pseudocode for this method is given below:
def add_burning_config(grid, height, width): for i in range(height): grid[i][0] += 1 grid[i][-1] += 1 for i in range(width): grid[0][i] += 1 grid[-1][i] += 1 def identity(width, height): grid = [[0 for _ in range(width)] for _ in range(height) ] while True: old_grid = clone(grid) add_burning_config(grid, height, width) stabilize(grid) if grids_equal(grid, old_grid): return grid
... continue reading