Find the identity sandpile

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 -> 6 0 6 -> 2 4 2 -> 3 0 3 -> 5 0 5 -> 1 4 1 -> 2 0 2 -> 4 0 4 -> 0 4 0 -> 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 -> 5 2 5 -> 1 6 1 -> 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.

• A word of warning to competitors: I described what the solution is, not how to compute it. Your code must be able to produce the 50 by 50 case in under a minute, so brute forcing is not a possibility. There are algorithms for solving this, and I'm not giving you them. That's intentional. I feel that too many challenges present you with pre-chewed food. I will give a +100 bounty to the first answer that doesn't trivialize the problem with a builtin (looking at you, Mathematica), at my discretion. – orlp Jan 15 '17 at 22:19
• I think the image of the 500x500 identity would benefit from saying what number each color corresponds to. – xnor Jan 15 '17 at 22:21
• What constitutes "sufficient contrast"? – Conor O'Brien Jan 15 '17 at 22:44
• @ConorO'Brien Any set of colors that are sufficiently distinguishable. I could make it 100% objective with some color measure, but I feel that's overkill. I don't care if you use red, yellow, green, grayscale, or whatever, just don't use 4 shades of grey that are within 1% of eachother (like #000000, #000001, #000002, #000003). – orlp Jan 15 '17 at 22:51
• ahem I believe this question is now eligible for bounty. Could I get the +100 bonus? :) – JungHwan Min Feb 1 '17 at 1:48

Octave, 120 113 bytes

function a=W(a);while nnz(b=a>3);a+=conv2(b,[t=[0 1 0];!t;t],'same')-b*4;end;end;@(n)imagesc(W((x=ones(n)*6)-W(x)))


Thanks to JungHwan Min for providing a link to the reference paper in his Mathematica answer.
Thanks to Stewie Griffin saved me 7 bytes [any(any(x)) -> nnz(x)]

Here two function is used:

1.f: for stabilization of a matrix
2. An anonymous function that takes n as input and shows the identity matrix.

Try It On rextester! for generation of a 50 * 50 matrix

Elapsed time for computation of the matrix: 0.0844409 seconds.

Explanation:

Consider a function f that stabilizes a matrix the task of finding the identity is simply

f(ones(n)*6 - f(ones(n)*6).

that ones(n)*6 means a n*n matrix of 6.

so for n=3:

M = [6 6 6
6 6 6
6 6 6];


The result will be f(M-f(M))

For stabilization function 2D convolution used to speed the task; In each iteration we make a binary matrix b with the same size of the input matrix and set it to 1 if corresponding element of the input matrix is >3. Then we apply a 2D convolution of the binary matrix with the following mask

0 1 0
1 0 1
0 1 0


representing four direct neighbors .
Result of convolution is added to the matrix and 4 times the binary matrix subtracted from it.

The loop continued until all element of the matrix are <= 3

Ungolfed version:

function a=stabilize(a)
while any(any(b=a>3))
end
end
n= 50;
M = ones(n)*6;
result = stabilize(M-stabilize(M));
imagesc(result);


Mathematica, 177157135 133 bytes

Colorize[f=BlockMap[⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&,#~ArrayPad~1,{3,3},1]&~FixedPoint~#&;k=Table[6,#,#];f[k-f@k]]&


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~#&


... 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&


... 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]&~FixedPoint~#&;k=Table[6,#,#];f[k-f@k]]&


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

• Considering that you followed the spirit of the time limit (implementing an actual algorithm rather than brute force), and it's very much plausible that a powerful desktop computer is 50% faster, I'll allow it. Since it implements an actual algorithm without a trivializing built-in, this qualifies for the +100 bonus. You'll have to wait for it though, as I can't start a bounty yet. – orlp Jan 16 '17 at 5:04
• That being said, implementing this quite trivially in Python (a notoriously slow language), only takes ~2s for n = 50. Maybe you can speed it up a bit? – orlp Jan 16 '17 at 5:15
• @orlp Done, but it's longer than the original code. Should I make the faster version my main answer or can I just put it at the end? – JungHwan Min Jan 16 '17 at 6:14
• Like this is fine I think. – orlp Jan 16 '17 at 6:30

Python 3 + Numpy + PIL, 385370 364 bytes

import numpy as q,PIL.Image as w
n=int(input())
z=n,n
def r(p):
while len(p[p>3]):
for x,y in q.ndindex(z):
if p[x,y]>3:
p[x,y]-=4;p[x-1,y]+=x>0;p[x,y-1]+=y>0
if~-n>x:p[x+1,y]+=1
if~-n>y:p[x,y+1]+=1
s=q.full(z,6)
t=s.copy()
r(t)
i=s-t
r(i)
w.fromarray(q.uint8(q.array(q.vectorize(lambda x:[x//1*65]*3,otypes=[object])(i).tolist()))).save('i.png')


Takes input on STDIN. Outputs the image as greyscale to i.png. Black corresponds to 0, dark grey to 1, light grey to 2, and white to 0.

Uses the formula I = R(S - R(S)), where I is the identity element, S is the matrix filled with sixes, and R is the reduction function.

I could probably save some bytes by switching to Python 2 and doing from numpy import*, but (1) I don't have Numpy installed on Python 2 and (2) the program wasn't terminating with from numpy import*.

Ungolfed:

import numpy as np
from PIL import Image

# Compute the identity element

n = int(input('Size of the sandpile: '))

def reduce_pile(sandpile):
while any(element >= 4 for element in np.nditer(sandpile)):
for x, y in np.ndindex((n, n)):
if sandpile[x, y] >= 4:
sandpile[x, y] -= 4
if x > 0: sandpile[x - 1, y] += 1
if y > 0: sandpile[x, y - 1] += 1
if x < n - 1: sandpile[x + 1, y] += 1
if y < n - 1: sandpile[x, y + 1] += 1

s = np.full((n, n), 6, dtype=np.int32)
s_prime = np.copy(s)

reduce_pile(s_prime)

identity = s - s_prime
reduce_pile(identity)

# Output it to identity.png as an image

colours = [[255, 255, 255], [255, 0, 0], [0, 255, 0], [0, 0, 255]]
img_array = np.vectorize(lambda x: colours[x], otypes=[object])(identity)
img_array = np.array(img_array.tolist(), dtype=np.uint8)

img = Image.fromarray(img_array)
img.save('identity.png')

• You may be able to save bytes by using scipy or matplotlib to display the data rather than generating an image explicitly with PIL. – orlp Jan 16 '17 at 20:18