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.


After 3 seconds, the grid looks like


After 4 seconds:


After 15 seconds:


And after 16 seconds:


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.


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


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


For 4, the output should be



Standard golf scoring applies.


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.

  • 2
    \$\begingroup\$ Can the input t be 0? What's the ouput then? \$\endgroup\$
    – Luis Mendo
    Commented Sep 4, 2016 at 21:46
  • 1
    \$\begingroup\$ 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? \$\endgroup\$
    – flawr
    Commented Sep 4, 2016 at 21:55
  • 2
    \$\begingroup\$ @flawr: Yes, that would be correct. Look at the difference between t=15 and t=16. \$\endgroup\$ Commented Sep 4, 2016 at 22:02
  • \$\begingroup\$ @LuisMendo The input is specified as positive t, so zero isn't a valid input. \$\endgroup\$ Commented Sep 4, 2016 at 22:28
  • 1
    \$\begingroup\$ Is it REALLY that necessary to have . for empty cells? Can we have 0 as a valid empty cell? \$\endgroup\$ Commented Sep 5, 2016 at 0:31

13 Answers 13


R, 378 343 297 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:


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

Test cases:

[1,]    3
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    1
[3,]    0    1    0
     [,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:

Abelian sandpile after 1000 steps

After 50000 drops (≈4 seconds):

Abelian sandpile after 50000 steps

After 333333 drops (≈15 minutes):

Abelian sandpile after 100000 steps

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 Growth rate 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):

Average iteration until the pile settles

Approximate computation time

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!

  • 1
    \$\begingroup\$ 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. \$\endgroup\$ Commented Sep 5, 2016 at 0:43
  • 1
    \$\begingroup\$ 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). \$\endgroup\$
    – Billywob
    Commented Sep 5, 2016 at 7:11
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – JDL
    Commented Sep 5, 2016 at 7:21
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – JDL
    Commented Sep 5, 2016 at 8:29
  • 1
    \$\begingroup\$ @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). \$\endgroup\$ Commented Sep 7, 2016 at 16:18

MATL, 55 53 48 43 42 bytes

Inspired by @flawr's answer.

Graphical output:


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:

enter image description here

ASCII output (43 bytes):


Try it online!


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
  +        %   Add to previous array. This updates the sandpile array
  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
  • 3
    \$\begingroup\$ I'm jealous of 1Y6. \$\endgroup\$
    – flawr
    Commented Sep 5, 2016 at 9:11
  • 1
    \$\begingroup\$ @flawr But your ~mod(spiral(3),2) is much more clever :-) \$\endgroup\$
    – Luis Mendo
    Commented Sep 5, 2016 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


As always:

Convolution is the key to success.

  • \$\begingroup\$ v=any(z) instead of v=find(sum(z)) (I'm using that in my answer). Also, 2*~z instead of (z<1)*2 \$\endgroup\$
    – Luis Mendo
    Commented Sep 4, 2016 at 23:18
  • \$\begingroup\$ My computer froze on input n=500... It had been processing n=400 for several seconds. Am I doing something wrong? \$\endgroup\$ Commented Sep 5, 2016 at 4:58
  • \$\begingroup\$ @AndreïKostyrka It works for me (Matlab R2015b) \$\endgroup\$
    – Luis Mendo
    Commented Sep 5, 2016 at 8:59
  • 1
    \$\begingroup\$ @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. \$\endgroup\$
    – flawr
    Commented Sep 5, 2016 at 9:13
  • \$\begingroup\$ @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). \$\endgroup\$
    – flawr
    Commented Sep 5, 2016 at 9:15

Python 2, 195 +1 +24 = 220 217

from pylab import*
ifrom scipy.signal import convolve2d as c
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


from pylab import*
from scipy.signal import convolve2d as c
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


  • \$\begingroup\$ 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? \$\endgroup\$ Commented Sep 5, 2016 at 20:31
  • \$\begingroup\$ 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 \$\endgroup\$
    – DenDenDo
    Commented Sep 6, 2016 at 12:02

Python 3 2, 418 385 362 342 330 bytes

w='[(i,j)for i in r(n)for j in r(n)if a[i][j]>3]'
def f(z):
 for _ in[0]*z:
  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.

  • \$\begingroup\$ I think you can move the assignment of a,x,r into the function arguments. \$\endgroup\$
    – xenia
    Commented Sep 5, 2016 at 5:25
  • 1
    \$\begingroup\$ 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? \$\endgroup\$
    – clismique
    Commented Sep 5, 2016 at 9:00
  • \$\begingroup\$ @Qwerp-Derp: Feel free! I'd love to see what you've done. \$\endgroup\$ Commented Sep 5, 2016 at 18:32

Perl, 157 147 bytes

Includes +1 for -p

Run with the count on STDIN, prints the map using 0123 to STDOUT:

sandpile.pl <<< 16


#!/usr/bin/perl -p
/?$.+=s%^|\z%0 x$..$/%eg+!s/\b/0/g:s^.^$&%4+grep{3<substr$\,0|$_+"@+",1}-$.-2,-2,0,$.^eg while/[4-7]/}($\="0

JavaScript, 418 416 406 400 393 bytes

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

var f =
<input id="i" type="number"><input type="button" value="Run" onclick="var v = +document.getElementById('i').value; if (v>0) f(v)">

  • 1
    \$\begingroup\$ Warning: I pressed 'run' without input, and my screen crashed (infinite loop). Don't be as silly as I was. \$\endgroup\$
    – roberrrt-s
    Commented Sep 5, 2016 at 19:47
  • 1
    \$\begingroup\$ @Roberrrt I updated my answer to prevent this. \$\endgroup\$
    – hetzi
    Commented Sep 6, 2016 at 13:00

Nim, 294 characters

import os,math,sequtils,strutils
 z=parseFloat paramStr 1
 b=y.newSeqWith newSeq[int] y
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


  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.


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)}


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.


  • t=10000 72ms
  • t=20000 167ms
  • t=30000 419ms
  • t=40000 659ms
  • t=100000 3413ms
  • t=1000000 about 6 minutes
  • \$\begingroup\$ 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? \$\endgroup\$ Commented Sep 6, 2016 at 13:36
  • \$\begingroup\$ 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 \$\endgroup\$ Commented Sep 6, 2016 at 13:43
  • \$\begingroup\$ 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. \$\endgroup\$ Commented Sep 6, 2016 at 14:05
  • \$\begingroup\$ Your performance can be beaten only by direct array manipulation in C or Fortran with compiler optimisation. I envy you. \$\endgroup\$ Commented Sep 6, 2016 at 19:40
  • \$\begingroup\$ @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 \$\endgroup\$ Commented Sep 6, 2016 at 20:00

J, 76 bytes

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 ^:_.


   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'
   load 'viewmat'
   ((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.



APL (Dyalog Unicode), 51 bytesSBCS

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

Try it online!


C 222

Thanks to @ceilingcat for some very nice pieces of golfing - now even shorter

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

Try it online!

  • \$\begingroup\$ Okay, I give up: why is your array 99 by 98? \$\endgroup\$ Commented Nov 18, 2016 at 14:26
  • \$\begingroup\$ @EricTressler How did I not find that in testing?! \$\endgroup\$ Commented Nov 19, 2016 at 19:13

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.