# Topple the sandpile

(There are related questions about infinite sandpiles, and finding identity elements of sandpiles.)

Given a matrix of non-negative integers, return a matrix of the same dimensions, but toppled:

1. If the matrix doesn't contain any values larger than 4, return it.
2. Every "cell" that is larger than 3 gets reduced by 4, and all directly neighboring cells (above, below, left, and right) are incremented, if they exist.
3. GOTO 1.

### Examples:

0 1 0        0 2 0
2 4 0   ->   3 0 1
0 0 3        0 1 3

1 2 3    2 3 4    2 5 1    4 1 2    0 3 3    0 3 3    0 3 3
4 5 6 -> 2 4 4 -> 4 2 3 -> 0 5 4 -> 3 2 1 -> 3 3 1 -> 3 3 2
7 8 9    5 7 7    2 6 5    4 3 2    0 5 3    1 1 4    1 2 0


(You only need to return the final result. The path on which you reach it may differ from the one shown here: it doesn't matter in which order you perform the toppling operations, they all lead to the same result.)

For a deeper explanation and some motivation see this Numberphile video or the Wikipedia article on the Abelian sandpile model.

## Rules:

• You may take input and output in any of the standard ways
• Loopholes are forbidden
• Input and output may be:
• a nested list: [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
• a simple list: [1, 2, 3, 4, 5, 6, 7, 8, 9] and the shape
• some kind of native matrix type
• a string, e.g. 1 2 3\n4 5 6\n7 8 9
• or whatever else works in your language.
• Input and output must be in the same form
• The input may contain larger numbers than the ones shown here, but the size may be bound by the limits of your language (MAXINT equivalents, if applicable)
• The matrix may have any shape (e.g. 1x1, 2x2, 3x3, 4x4, 2x7, 11x3, ...)
• You don't need to handle the case where the shape is 0xN or Nx0.

## Testcases

[[2, 5, 4], [8, 6, 4], [1, 2, 3]] -> [[3, 3, 0], [1, 2, 2], [1, 3, 2]]
[[0, 0, 2], [1, 3, 3], [0, 0, 0]] -> [[0, 0, 2], [1, 3, 3], [0, 0, 0]]
[[9, 9, 9], [9, 9, 9], [9, 9, 9]] -> [[1, 3, 1], [3, 1, 3], [1, 3, 1]]
[[4, 5], [2, 3]] -> [[2, 3], [0, 1]]
[[2, 3, 5], [2, 2, 0]] -> [[3, 0, 2], [2, 3, 1]]
[[7]] -> [[3]]


This is , the shortest code (per language) wins.

• Is it OK to display all the intermediate results? Mar 30 '17 at 13:29
• @feersum I guess so, as long as it's clear what the final result is. Mar 30 '17 at 13:38

# MATL, 17 bytes

tss:"t3>t1Y6Z+w4*-+


Try it at MATL Online! Or verify all test cases.

### Explanation

The program iterates for as many times as the sum of the input. This is a loose upper bound on the required number of iterations.

For each iteration, entries in the sandpile matrix exceeding 3 are detected, giving a matrix of 1 and 0, which is convolved with the 4-neighbour mask. The entries exceeding 3 in the sandpile matrix are reduced by 4, and the result of the convolution is added.

For the last iterations, in which the sandpile matrix does not have any numbers exceeding 3, zeros are subtracted from and added to it, so it is unaffected.

t       % Implicit input (matrix). Duplicate
ss      % Sum of matrix entries
:"      % Repeat that many times
t     %   Duplicate
3>    %   True for matrix entries that exceed 3
t     %   Duplicate
1Y6   %   Push predefined literal [0, 1, 0; 1, 0, 1; 0, 1, 0]
Z+    %   2D convolution, keeping size
w     %   Swap
4*    %   Multiply by 4
-     %   Subtract
% Implicit end. Implicit display

• Convolution high five. Mar 30 '17 at 13:41
• @MartinEnder Ah, you also used that :-) Nice to see a fellow convoluter! I'm sure flawr will join us soon Mar 30 '17 at 13:43
• @LuisMendo Convolutionista Mar 30 '17 at 14:43

## Mathematica, 65 bytes

#//.s_:>s+ListConvolve[{v={0,1,0},1-v,v},x=UnitStep[s-4],2,0]-4x&


### Explanation

#//.s_:>...&


Repeatedly transform the input by toppling all piles greater than 3. This process stops automatically when the transformation fails to change the matrix (i.e. when no large piles exist any more). In the following expression the matrix is called s.

...x=UnitStep[s-4]...


Create a matrix that has a 1 whenever the current matrix has a 4 or greater, and a zero otherwise. This is essentially a mask that indicates which piles need to be toppled. Call the mask x.

ListConvolve[{v={0,1,0},1-v,v},x=UnitStep[s-4],2,0]


First we compute the number of sand that is added to each pile due to toppled neighbouring piles. This is done with a convolution of the following matrix over x:

0 1 0
1 0 1
0 1 0


Essentially, it adds one to the current cell for each of its von-Neumann neighbours in the mask.

s+...-4x


We add the previous result to s and then we subtract four times the mask from it to reduce the toppled piles.

# Octave, 65 bytes

This doesn't seem very good, I must be missing some tricks...

m=input(0);do;m+=conv2(m>3,[0 1 0;1 -4 1;0 1 0],"same")
until m<4

• What version of Octave are you using that allows input(0)? Mar 30 '17 at 14:44
• @Suever >> version ans = 4.0.1 Mar 30 '17 at 14:46

## JavaScript (ES6), 101 95 bytes

Takes the width of the matrix w and an array of values a in currying syntax (w)(a). Returns an array of values.

w=>g=a=>(b=a.map((n,i)=>n%4+(F=d=>~m|i%w&&a[i+d]>>2)(m=w)+F(-w)+F(m=-1)+F(!++i)))+0==a+0?a:g(b)


### Formatted and commented

w =>                      // main function: takes w as input, returns g
g = a =>                // recursive function g: takes a as input
(                     //
b = a.map((n, i) => // for each element n at position i in a:
n % 4 + (         //   keep only n MOD 4
F = d =>        //   define F(): function that takes d as input
~m |          //     if m is not equal to -1
i % w &&      //     or i MOD w is not null:
a[i + d] >> 2 //       return a fourth of the value of the cell at i + d
)(m = w) +        //   test the cell below the current cell
F(-w) +           //   test the cell above
F(m = -1) +       //   test the cell on the left
F(!++i)           //   test the cell on the right
)                   // end of map(): assign the result to b
) + 0 == a + 0 ?      // if b is equal to a:
a                   //   stop recursion and return a
:                     // else:
g(b)                //   do a recursive call with b


### Test cases

let f =

w=>g=a=>(b=a.map((n,i)=>n%4+(F=d=>~m|i%w&&a[i+d]>>2)(m=w)+F(-w)+F(m=-1)+F(!++i)))+0==a+0?a:g(b)

console.log(JSON.stringify(f(3)([2, 5, 4, 8, 6, 4, 1, 2, 3]))); // -> [3, 3, 0, 1, 2, 2, 1, 3, 2]
console.log(JSON.stringify(f(3)([0, 0, 2, 1, 3, 3, 0, 0, 0]))); // -> [0, 0, 2, 1, 3, 3, 0, 0, 0]
console.log(JSON.stringify(f(3)([9, 9, 9, 9, 9, 9, 9, 9, 9]))); // -> [1, 3, 1, 3, 1, 3, 1, 3, 1]
console.log(JSON.stringify(f(2)([4, 5, 2, 3]))); // -> [2, 3, 0, 1]
console.log(JSON.stringify(f(3)([2, 3, 5, 2, 2, 0]))); // -> [3, 0, 2, 2, 3, 1]
console.log(JSON.stringify(f(1)([7]))); // -> [3]

## JavaScript (ES6), 118114 104 bytes

Saved 2 bytes thanks to @Neil

f=a=>a.find(b=>++y&&b.find(c=>++x&&c>3,x=0),y=0)?f(a.map(b=>b.map(c=>c+=--i|y?i*i+y*y==1:-4,i=x,--y))):a

• Does (i-=x)|y-j?i*i+ help?
– Neil
Mar 30 '17 at 14:20
• @Neil It does indeed, thanks! Mar 30 '17 at 14:23
• ...I was on the phone but I was also considering a.find(...b.find(...c>3&&a.map(...)))&&f(a).
– Neil
Mar 30 '17 at 15:11
• @Neil I don't think that would work, since .map doesn't mutate... Mar 30 '17 at 15:34
• It seems that making it mutate costs slightly less than moving the map inside the find saves: f=a=>a.find((b,x)=>b.find((c,y)=>c>3&&a.map(b=>b.map((_,j)=>b[j]+=x|(j-=y)?x*x+j*j==1:-4)&x--)))&&f(a)
– Neil
Mar 30 '17 at 20:24

# C++, 261258 250 bytes

#import<vector>
#define S size()
void f(std::vector<std::vector<int>>&m){s:int i,j,r;for(i=r=0;i<m.S;++i)for(j=0;j<m[i].S;++j){if(m[i][j]>3){r=1;m[i][j]-=4;j>0&&m[i][j-1]++;i>0&&m[i-1][j]++;j<m[i].S-1&&m[i][j+1]++;i<m.S-1&&m[i+1][j]++;}}if(r)goto s;}


Takes input as a reference to a vector of vectors and modifies it directly.

Try it online!

• 237 bytes Jul 29 '20 at 1:35