Find the identity sandpile
Mathematica, 177 157 135 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.
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)
mask = [t=[0 1 0];!t;t];
while any(any(b=a>3))
a+=conv2(b,mask,'same')-b*4;
end
end
n= 50;
M = ones(n)*6;
result = stabilize(M-stabilize(M));
imagesc(result);