# Build a sandpile

An abelian sandpile, for our purposes, is an infinite grid with integer coordinates, initially empty of sand. After each second, a grain of sand is placed at (0,0). Whenever a grid cell has 4 or more grains of sand, it spills one grain of sand to each of its four neighbors simultaneously. The neighbors of (x,y) are (x-1,y), (x+1,y), (x,y-1), and (x,y+1).

When a cell spills, it may cause its neighbors to spill. Some facts:

• This cascade will eventually stop.
• The order in which the cells spill is irrelevant; the result will be the same.

## Example

After 3 seconds, the grid looks like

.....
.....
..3..
.....
.....


After 4 seconds:

.....
..1..
.1.1.
..1..
.....


After 15 seconds:

.....
..3..
.333.
..3..
.....


And after 16 seconds:

..1..
.212.
11.11
.212.
..1..


## The challenge

In as few bytes as possible, write a function that takes a single positive integer t and outputs a picture of the sandpile after t seconds.

## Input

A single positive integer t, in any format you choose.

## Output

A picture of the sandpile after t seconds, using the characters

 . 1 2 3


Edit: Use any four distinct characters you like, or draw a picture. If you're not using ".123" or "0123", specify in your answer what the characters signify.

Unlike in the examples, your output should contain the minimal number of rows and columns necessary to show the nonzero part of the sandpile.

That is, for input 3, the output should be

 3


For 4, the output should be

 .1.
1.1
.1.


## Scoring

Standard golf scoring applies.

## Rules

No language functions or libraries that already know what a sandpile is are allowed.

Edit: The output section has been edited, the character set restriction has been completely lifted. Use any four distinct characters or colors you like.

• Can the input t be 0? What's the ouput then? – Luis Mendo Sep 4 '16 at 21:46
• Is it correct for at a given timestep multiple cascades can take place in a row? So in that timestep the cascades keep happening until every cell is again 3 or less? – flawr Sep 4 '16 at 21:55
• @flawr: Yes, that would be correct. Look at the difference between t=15 and t=16. – El'endia Starman Sep 4 '16 at 22:02
• @LuisMendo The input is specified as positive t, so zero isn't a valid input. – Eric Tressler Sep 4 '16 at 22:28
• Is it REALLY that necessary to have . for empty cells? Can we have 0 as a valid empty cell? – Andreï Kostyrka Sep 5 '16 at 0:31

## R, 378343297 291 bytes

As usually, the user supplies his/her input via scan() (I already used the variable t, so let us take z instead), so the second line should be launched separately, and then the rest:

e=numeric
a=1%*%scan()
x=1
o=a>3
n=1
while(any(o)){
v=which(o,T)
if(any(v==1)){a=rbind(e(n+2),cbind(e(n),a,e(n)),e(n+2));x=x+1;n=n+2;v=which(a>3,T)}
q=nrow(v)
u=cbind(e(q),1)
l=v-u[,1:2];r=v+u[,1:2];t=v-u[,2:1];b=v+u[,2:1]
a[l]=a[l]+1;a[r]=a[r]+1;a[t]=a[t]+1;a[b]=a[b]+1
a[v]=a[v]-4
o=a>3}
a


Outputs an array that contains the values of a at tth generation (0, 1, 2 or 3).

Test cases:

z=3
[,1]
[1,]    3
z=4
[,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    1
[3,]    0    1    0
z=16
[,1] [,2] [,3] [,4] [,5]
[1,]    0    0    1    0    0
[2,]    0    2    1    2    0
[3,]    1    1    0    1    1
[4,]    0    2    1    2    0
[5,]    0    0    1    0    0


It helps us that this thing is symmetrical both vertcially and horizontally, which means that is the leftmost point gets a height of 4, this means that the topmost, rightmost and lowest points are also 4.

Oh, and did I say that you can make beautiful visualisations?

After 1000 drops:

After 50000 drops (≈4 seconds):

After 333333 drops (≈15 minutes):

You can draw it, too!

image(1:n,1:n,a,col=colorRampPalette(c("#FFFFFF","#000000"))(4), axes=F, xlab="", ylab="")


This thing took 4 seconds for 10000 iterations but slows down considerably for larger array sizes (e.g. a couple of minutes for 100000 iterations). This is why it gets so slow (I estimated the growth rate as in and obtained τ(i)≈689·i^1.08, so the average time per one additional grain until the whole sandpile settles after ith step is slightly larger than one), and the total time as a function of the number of grains grows a little bit slower than quadratically (T(i)≈0.028*i^1.74):

And now with a complete explanation:

e=numeric # Convenient abbreviation for further repeated use
a=1%*%scan() # Creates a 1×1 array with a user-supplied number
x=1 # The coordinate of the centre
o=a>3 # Remember which cells were overflown
n=1 # Array height that is going to change over time
while(any(o)){ # If there is still any overflow
v=which(o,T) # Get overflown cells' indices
if(any(v==1)){ # If overflow occurred at the border, grow the array
a=rbind(e(n+2),cbind(e(n),a,e(n)),e(n+2)) # Growing
x=x+1 # Move the centre
n=n+2 # Change the height
v=which(a>3,T) # Re-index the overflowed cells
}
q=nrow(v) # See how many indices are overflown
u=cbind(e(q),1) # Building block for neighbours' indices
l=v-u[,1:2];r=v+u[,1:2];t=v-u[,2:1];b=v+u[,2:1] # L, R, T, B neighbours
a[l]=a[l]+1;a[r]=a[r]+1;a[t]=a[t]+1;a[b]=a[b]+1 # Increment neighbours
a[v]=a[v]-4 # Remove 4 grains from the overflown indices
o=a>3} # See if still overflown indices remain
a # Output the matrix


This is the first time in my life when growing objects (like a <- c(a, 1)) works so much faster than pre-allocating a big empty matrix for values and filling it gradually with a ton of zeros unused.

Update. Golfed 18 bytes by removing arr.ind in which owing to Billywob and replacing rep(0,n) with e=numeric;e(n) in 5 instances owing to JDL, and 17 more bytes owing to JDL.

Update 2. Since the sandpile is Abelian, it may start with a stack of desired height, so I removed the redundant loop and gained a huge boost in productivity!

• I get your point about the extra column,row indices that you're outputting, but I think I want to restrict the output to just being "the answer", and nothing more. I'm glad you included pictures, though. – Eric Tressler Sep 5 '16 at 0:43
• Nice answer Andreï! You could definitely golf a few bytes though, for instance predefining rep(), given that you use it 6 times. Secondly, I don't think you need to write out the arr.ind=T option for the which() function. Simply use which(...,T). – Billywob Sep 5 '16 at 7:11
• It may be golfier to define n=numeric and use that instead, since n(k) is fewer characters than r(0,k). I like the pictures though. – JDL Sep 5 '16 at 7:21
• Another suggestion: 1%*%0 is fewer characters than array(0,c(1,1)). Also the second argument to u <- cbind can just be 1, cbind will extend it to the length of the first argument by default. – JDL Sep 5 '16 at 8:29
• @GregMartin Fixed this. Sorry for that; in my first language, we use the word “self’s” and never bother about the gender of the person in question (like “one small step for a man”); still, sometimes, at very rare occasions, I call a dog “she” or a “he”, whereas it should be “it”, unless, like, you are the owner and you really want to emphasize the sex of your anumal (in spite of the fact that telling a male from a female is not that hard). – Andreï Kostyrka Sep 7 '16 at 16:18

# MATL, 55534843 42 bytes

Graphical output:

0i:"Gto~+XytP*+t"t4=t1Y6Z+b+w~*]]tat3$)1YG  Try it at MATL Online!. It takes about 10 seconds for input 30. You may need to refresh the page and press "Run" again if it doesn't work. Here's an example result for input 100: ASCII output (43 bytes): 0i:"Gto~+XytP*+t"t4=t1Y6Z+b+w~*]]tat3$)48+c


Try it online!

### Explanation

0          % Push a 0. This is the initial array. Will be resized in first iteration
i:         % Take input n. Generate range [1 2 ... n]
"          % For each, i.e. repeat n times
Gto~+    %   Push input and add negate parity. This "rounds up" n to odd number
%   m = n or n+1
Xy       %   Identity matrix with that size
tP*      %   Multiply element-wise by vertically flipped copy. This produces a
%   matrix with a 1 in the center and the rest entries equal to 0
t        %   Duplicate
"        %   For each column (i.e. repeat m times)
t4=    %     Duplicate. Compare with 4 element-wise. This gives a 2D mask that
%     contains 1 for entries of the sandpile array that equal 4, and 0
%     for the rest
t      %     Duplicate
1Y6    %     Predefined literal: [0 1 0; 1 0 1; 0 1 0]
Z+     %     2D convolution, maintaining size
b      %     Bubble up to bring sandpile array to top
+      %     Element-wise addition. This adds 1 to the neighbours of a 4
w      %     Swap to bring copy of mask to top
~*     %     Multiply bu negated mask. This removes all previous 4
]        %  End
]          % End
t          % Duplicate the updated sandpile array
a          % 1D mask that contains 1 for columns that contain a 1. This will be
% used as a logical index to select columns
t          % Duplicate. This will be used as logical index to select rows (this
% can be done because of symmetry)
3$) % Keep only those rows and columns. This trims the outer zeros in the % sandpile array 1YG % Display as scaled image  • I'm jealous of 1Y6. – flawr Sep 5 '16 at 9:11 • @flawr But your ~mod(spiral(3),2) is much more clever :-) – Luis Mendo Sep 5 '16 at 9:22 # Matlab, 160 156 148 bytes n=input('');z=zeros(3*n);z(n+1,n+1)=n;for k=1:n;x=z>3;z=z+conv2(+x,1-mod(spiral(3),2),'s');z(x)=z(x)-4;end;v=find(sum(z));z=z(v,v);[z+48-(z<1)*2,'']  First a way too big matrix is created, with the n in the middle somewhere. Then the cascading is calculated with a very convenient 2d convolution. At the end the excess is trimmed off and the whole thing is converted to a string. Example output for t=100 ...121... ..32.23.. .3.323.3. 123.3.321 2.23.32.2 123.3.321 .3.323.3. ..32.23.. ...121...  As always: Convolution is the key to success. • v=any(z) instead of v=find(sum(z)) (I'm using that in my answer). Also, 2*~z instead of (z<1)*2 – Luis Mendo Sep 4 '16 at 23:18 • My computer froze on input n=500... It had been processing n=400 for several seconds. Am I doing something wrong? – Andreï Kostyrka Sep 5 '16 at 4:58 • @AndreïKostyrka It works for me (Matlab R2015b) – Luis Mendo Sep 5 '16 at 8:59 • @AndreïKostyrka For an input of n this program generates an 3*n x 3*n matrix, so it needs to store about 9*n^2 numbers. Also it is totally inefficient, because we also have a totally unnecessary long iteration from 1 upto n. But then again this is code-golf, making a program efficient is a different cup of tea. – flawr Sep 5 '16 at 9:13 • @AndreïKostyrka You can makie it more memory efficient by using sparse matrices (second line: z=sparse(zeros(2*n+1))) and changing the for loop to while any(z(:)>3). Then you can also perhaps compute the convolution kernel just once: kern = 1-mod(spiral(3),2). – flawr Sep 5 '16 at 9:15 # Python 2, 195 +1 +24 = 220 217 from pylab import* ifrom scipy.signal import convolve2d as c k=(arange(9)%2).reshape(3,3) def f(n):g=zeros((n,n),int);g[n/2,n/2]=n;exec"g=c(g/4,k,'same')+g%4;"*n;return g[any(g,0)].T[any(g,0)]  output for n=16 array([[0, 0, 1, 0, 0], [0, 2, 1, 2, 0], [1, 1, 0, 1, 1], [0, 2, 1, 2, 0], [0, 0, 1, 0, 0]])  there is a LOT of unnecessary padding and iterations, by using n as a "good-enough" upper bound, but n=200 still completed in a second and n=500 in about 12 seconds ungolfed from pylab import* from scipy.signal import convolve2d as c k=array([0,1,0], [1,0,1], [0,1,0]) def f(n): g=zeros((n,n)) # big grid of zeros, way bigger than necessary g[n/2,n/2]=n # put n grains in the middle exec"g=c(g/4,k,'same')+g%4;"*n # leave places with <4 grains as is, convolve the rest with the kernel k, repeat until convergence (and then some more) return g[any(g,0)].T[any(g,0)] # removes surrounding 0-rows and columns  replacing return x by imshow(x) adds one character and outputs an ugly interpolated image, adding imshow(x,'gray',None,1,'nearest') removes the blurry interpolation bringing the output up to specs • I get the following error when I run your code: ImportError: No module named convolve2d. Changing import scipy.signal.convolve2d as c to from scipy.signal import convolve2d as c resolves the issue. I am using scipy version 0.16.1, do I need an older or newer version? Or is the issue something else? – Andrew Epstein Sep 5 '16 at 20:31 • weird, now that you mention that it doesnt work for me anymore either. I probably did it right the first time in interactive mode, then shortened it and ignored the error, but the function stayed in memory – DenDenDo Sep 6 '16 at 12:02 # Perl, 157 147 bytes Includes +1 for -p Run with the count on STDIN, prints the map using 0123 to STDOUT: sandpile.pl <<< 16  sandpile.pl: #!/usr/bin/perl -p map{++substr$_,y///c/2-1,1;/4
/?$.+=s%^|\z%0 x$..$/%eg+!s/\b/0/g:s^.^$&%4+grep{3<substr$\,0|$_+"@+",1}-$.-2,-2,0,$.^eg while/[4-7]/}($\="0 ")x$_}{


## Python 3 2, 418385362342 330 bytes

w='[(i,j)for i in r(n)for j in r(n)if a[i][j]>3]'
def f(z):
a,x,r=[[z]],0,range
for _ in[0]*z:
n=len(a);v=eval(w)
if[1for b,c in v if(b==0)+(c==0)]:n+=2;a=[[0]*n]+[[0]+a[i]+[0]for i in r(n-2)]+[[0]*n];x+=1;v=eval(w)
for c,d in v:exec'a[c+%s][d+%s]+=1;'*4%(-1,0,1,0,0,-1,0,1);a[c][d]-=4
for i in a:print''.join(map(str,i))


Edit: saved 6 bytes thanks to @Qwerp-Derp

All credit to @Andreï Kostyrka, as this is a direct translation of his R code into Python.

• I think you can move the assignment of a,x,r into the function arguments. – Loovjo Sep 5 '16 at 5:25
• I've golfed down your code by a few bytes... it's not much, but it'll have to do. Do you mind if I put an edit on your answer and if I change the version of Python to Python 2? – Qwerp-Derp Sep 5 '16 at 9:00
• @Qwerp-Derp: Feel free! I'd love to see what you've done. – Andrew Epstein Sep 5 '16 at 18:32

# JavaScript, 418416406400 393 bytes

Creates an anonymous function that displays the output on the console.

var f =
t=>{a=(q,w)=>Math.max(q,w);c=_=>{x=a(p[0],x);y=a(p[1],y);m[p]=(g(p)+1)%4;if(!m[p]){s.push([p[0],p[1]]);}};x=y=0,m={};g=k=>{v=m[k];return!v?0:v;};m[o=[0,0]]=1;s=[];while(--t){m[o]=(m[o]+1)%4;if(!m[o]){s.push(o);}while(s.length){p=s.pop();p[0]++;c();p[0]-=2;c();p[0]++;p[1]++;c();p[1]-=2;c();p[1]++;}}s='';for(i=-x;i<=x;i++){for(j=-y;j<=y;j++){v=g([i,j]);s+=v==0?'.':v;}s+='\n';}console.log(s);}
<input id="i" type="number"><input type="button" value="Run" onclick="var v = +document.getElementById('i').value; if (v>0) f(v)">

• Warning: I pressed 'run' without input, and my screen crashed (infinite loop). Don't be as silly as I was. – Roberrrt Sep 5 '16 at 19:47
• @Roberrrt I updated my answer to prevent this. – hetzi Sep 6 '16 at 13:00

# Nim, 294 characters

import os,math,sequtils,strutils
var
z=parseFloat paramStr 1
y=z.sqrt.toInt+1
w=y/%2
b=y.newSeqWith newSeq[int] y
x=0
proc d(r,c:int)=
b[r][c]+=1;if b[r][c]>3:b[r][c]=0;d r-1,c;d r,c+1;d r+1,c;d r,c-1
for i in 1..z.toInt:d w,w
while b[w][x]<1:x+=1
for r in b[x..< ^x]:echo join r[x..< ^x]


Compile and Run:

nim c -r sandbox.nim 1000


Notes:

1. I was able to come up with a shorter version which uses a fixed table size, but I edited it out in favor of the dynamic one.
2. Once the sandbox is computed, x is calculated as the number of zero columns at the start of the middle row.
3. For display, the table is sliced down by excluding x rows and columns from each end.

## Performance

nim c --stackTrace:off --lineTrace:off --threads:off \
--checks:off --opt:speed sandbox.nim

time ./sandbox   10000       0.053s
time ./sandbox   20000       0.172s
time ./sandbox   30000       0.392s
time ./sandbox   40000       0.670s
time ./sandbox  100000       4.421s
time ./sandbox 1000000    6m59.047s


# Scala, 274 bytes

val t=args(0).toInt
val s=(Math.sqrt(t)+1).toInt
val (a,c)=(Array.ofDim[Int](s,s),s/2)
(1 to t).map{_=> ?(c,c)}
println(a.map{_.mkString}.mkString("\n"))
def?(b:Int,c:Int):Unit={
a(b)(c)+=1
if(a(b)(c)<4)return
a(b)(c)=0
?(b+1,c)
?(b-1,c)
?(b,c+1)
?(b,c-1)
}


Usage:

scala sandpile.scala <iterations>

I don't think there is a lot to explain about this one. Basically it just adds one grain of sand to the center. Then checks if it is larger than 4, if so it'll spill over and check for all neighbours larger than 4, spill over, etc. It's quite fast.

Performance:

• t=10000 72ms
• t=20000 167ms
• t=30000 419ms
• t=40000 659ms
• t=100000 3413ms
• My program suggests that, centered at (0,0), the sandpile first hits a radius of 15 at t=1552. That would require a 31x31 array to store (coordinates -15 to 15 inclusive). Are you sure this is correct through t=5000? – Eric Tressler Sep 6 '16 at 13:36
• I am not sure this is correct, though I think I got the logic right? I get an array index out of bounds exception on t>5593 – AmazingDreams Sep 6 '16 at 13:43
• When I increment and then immediately check for spillage it does go out of bounds at t=1552. I'd say that is the correct implementation. I updated the code. – AmazingDreams Sep 6 '16 at 14:05
• Your performance can be beaten only by direct array manipulation in C or Fortran with compiler optimisation. I envy you. – Andreï Kostyrka Sep 6 '16 at 19:40
• @AndreïKostyrka, yes, this is where scala shines! My output doesn't conform to the specs though so I'd have to work on that – AmazingDreams Sep 6 '16 at 20:00

# J, 76 bytes

p=:0,.~0,.0,~0,]
p(4&|+3 3([:+/@,*&(_3]\2|i.9));._3[:<.4%~p)@.([:*/4>{.)^:_


I define a verb p which pads a border of zeros around the input. The main verb takes an array as input. It then checks the first row for any sandpile that contains 4 or more grains. If one does, it outputs the same array except padded using p, and otherwise it performs 2d convolution to simulate falling grains. The main verb is repeated until convergence using the power operator ^:_.

## Usage

   p =: 0,.~0,.0,~0,]
f =: p(4&|+3 3([:+/@,*&(_3]\2|i.9));._3[:<.4%~p)@.([:*/4>{.)^:_
f 15
0 3 0
3 3 3
0 3 0
f 50
0 0 0 1 0 0 0
0 0 3 1 3 0 0
0 3 2 2 2 3 0
1 1 2 2 2 1 1
0 3 2 2 2 3 0
0 0 3 1 3 0 0
0 0 0 1 0 0 0
timex 'r =: f 50000'
46.3472
((256#.3&#)"0<.255*4%~i._4) viewmat r


It takes about 46 seconds to compute the result for n = 50000, and the result can be displayed using the viewmat addon with a monochrome color scheme.

# C 229 (with lots of warnings)

G[99][99],x,y,a=99,b=99,c,d;S(x,y){if(++G[y][x]>3)G[y][x]=0,S(x+1,y),S(x-1,y),S(x,y+1),S(x,y-1);a=x<a?x:a;b=y<b?y:b;c=x>c?x:c;d=y>d?y:d;}F(t){for(;t--;)S(49,49);for(y=b;y<=d;y++){for(x=a;x<=c;x++)printf("%d ",G[y][x]);puts("");}}

/* call it like this */
main(_,v)char**v;{F(atoi(v[1]));}

• Okay, I give up: why is your array 99 by 98? – Eric Tressler Nov 18 '16 at 14:26
• @EricTressler How did I not find that in testing?! – Jerry Jeremiah Nov 19 '16 at 19:13

# APL (Dyalog Unicode), 51 bytesSBCS

{m⌿r/⍨m←⌈/×r←(4∘|+{+/⊢/4 2⍴⌊⍵÷4}⌺3 3)0,∘⍉∘⌽⍣4⊢⍵}⍣≡⍪


Try it online!

# PHP, 213 bytes

function d($x,$y){global$p,$m;$m=max($m,$x);$q=&$p[$y][$x];if(++$q>3){$q=0;d($x+1,$y);d($x-1,$y);d($x,$y+1);d($x,$y-1);}}while($argv[1]--)d(0,0);for($y=-$m-1;$y++<$m;print"\n")for($x=-$m;$x<=$m;)echo+$p[$y][$x++];  recursively creates the pile in $p, remembering the size in \$m; then prints with nested loops.
Run with -r.