Build a sandpile

Question

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.

Answer 1

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:

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!

Answer 2

MATL, 55 53 48 43 42 bytes

Inspired by @flawr's answer.

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
  +        %   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

Answer 3

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.

Answer 4

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

Answer 5

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):
 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.

Answer 6

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$_}{

Answer 7

JavaScript, 418 416 406 400 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)">

Answer 8

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

Answer 9

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
• t=1000000 about 6 minutes

Answer 10

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

Answer 11

APL (Dyalog Unicode), 51 bytes^SBCS

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

Try it online!

Answer 12

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!

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