Find the identity sandpile This question is about abelian sandpiles. Read this previous challenge and watch this numberphile video to learn more.



An abelian sandpile of size n by n is a grid containing the number 0, 1, 2 and 3 (representing the number of grains of sand). Adding two sandpiles works by first adding element by element, and then toppling any element that goes above 3. The order in which you topple doesn't matter, the end result is the same. When a cell topples its number decreases by 4, and each of its direct neighbors increases by 1. This can cause a chain reaction. If a cell is on the edge of the grid, any grains that fall off the grid while toppling disappear.

For example, I'm adding two 3 by 3 sandpiles (giving a rather extreme chain reaction):

3 3 3   1 2 1   4 5 4    4 6 4    6 2 6    6 3 6    2 5 2    4 1 4    4 2 4    0 4 0    2 0 2    2 1 2
3 3 3 + 2 1 2 = 5 4 5 -&gt; 6 0 6 -&gt; 2 4 2 -&gt; 3 0 3 -&gt; 5 0 5 -&gt; 1 4 1 -&gt; 2 0 2 -&gt; 4 0 4 -&gt; 0 4 0 -&gt; 1 0 1
3 3 3   1 2 1   4 5 4    4 6 4    6 2 6    6 3 6    2 5 2    4 1 4    4 2 4    0 4 0    2 0 2    2 1 2


In this challenge we're interested in a subset of all possible n by n sandpiles. This subset contains any sandpile that you can get by adding an arbitrary sandpile to the all-3s n by n sandpile. For example, just above we saw that 212 | 101 | 212 is in the subset, because we got it by adding something to the all-3 sandpile.

Now this subset has an interesting element: the identity element. If you take this element and add it to any other element in the subset, the sum is unchanged. In other words, this sandpile acts like a zero of this subset. It just so happens that 212 | 101 | 212 is the zero element for the subset of 3 by 3. For example:

2 2 2   2 1 2   4 3 4    0 5 0    2 1 2    2 2 2
2 2 2 + 1 0 1 = 3 2 3 -&gt; 5 2 5 -&gt; 1 6 1 -&gt; 2 2 2
2 2 2   2 1 2   4 3 4    0 5 0    2 1 2    2 2 2


Now this is your challenge: given n, find the identity element of the subset of the n by n grid. Output it by assigning a unique color with sufficient contrast of your choice to each of 0, 1, 2, 3 and outputting an n by n image. Your code must be able to produce the 50 by 50 case in under a minute on a reasonable modern PC.



For example, the 500 by 500 identity element:



Here is blue = 3, green = 2, red = 1, white = 0. But you don't have to use this color scheme in your answer.
Mathematica, 177 157 135 133 bytes
Colorize[f=BlockMap[⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&amp;,#~ArrayPad~1,{3,3},1]&amp;~FixedPoint~#&amp;;k=Table[6,#,#];f[k-f@k]]&amp;

Takes a number n. The output is the identity sandpile. 0 is black, 1 is light gray, 2 is magenta, and 3 is blue-gray.
Sadly, Mathematica does not have a builtin for this...
Uses the algorithm stated in Scott Corry and David Perkinson's paper.
Takes 91.7 seconds on my 5-year-old laptop to calculate the 50x50 identity sandpile. I am confident that a reasonable modern desktop computer is more than 50% faster. (I also have a way faster code at the end).
Explanation
f= ...

Define function f (the input is a sandpile matrix): a function that ...
BlockMap[ ... ]~FixedPoint~#&amp;

... repeats the BlockMap operation until the output does not change. BlockMap operation: ...

#~ArrayPad~1

... pad the input array with one layer of 0 ...
{3,3},1

... partition it into 3x3 matrices, with offset 1 ...
⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&amp;

... and for each partition, add the number of sand grains toppled onto the center cell and the center cell value mod 4.
i.e. the output of f is the stabilized version of the input.

k=Table[6,#,#]

Define k as an n by n array of 6s.
f[k-f@k]]

Compute f(k - f(k)).
Colorize[ ... ]

Apply colors to the result.
Faster version (142 bytes)
Colorize[r=RotateLeft;a=ArrayPad;f=a[b=#~a~1;b+r[g=⌊b/4⌋,s={0,1}]+g~r~-s+r[g,1-s]+r[g,s-1]-4g,-1]&amp;~FixedPoint~#&amp;;k=Table[6,#,#];f[k-f@k]]&amp;

Same code, but uses built-in list rotation instead of BlockMap. Computes n=50 in 4.0 seconds on my laptop.