[Question]
  [
An [abelian sandpile](https://en.wikipedia.org/wiki/Abelian_sandpile_model), 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]
## 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 `t`th 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:
[![Abelian sandpile after 1000 steps](https://i.stack.imgur.com/6IGcD.png)](https://i.stack.imgur.com/6IGcD.png)
After 50000 drops (≈4 seconds):
[![Abelian sandpile after 50000 steps](https://i.stack.imgur.com/9awir.png)](https://i.stack.imgur.com/9awir.png)
After 333333 drops (≈15 minutes):
[![Abelian sandpile after 100000 steps](https://i.stack.imgur.com/gNEDM.png)](https://i.stack.imgur.com/gNEDM.png)
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](https://i.stack.imgur.com/5h0xt.png)](https://i.stack.imgur.com/5h0xt.png) and obtained τ(i)≈689·i^1.08, so the average time per one additional grain until the whole sandpile settles after `i`th 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](https://i.stack.imgur.com/dHfnE.png)](https://i.stack.imgur.com/dHfnE.png)
[![Approximate computation time](https://i.stack.imgur.com/cr3zX.png)](https://i.stack.imgur.com/cr3zX.png)
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](https://codegolf.stackexchange.com/users/56989/billywob) and replacing `rep(0,n)` with `e=numeric;e(n)` in 5 instances owing to [JDL](https://codegolf.stackexchange.com/users/55516/jdl), and 17 more bytes owing to [JDL](https://codegolf.stackexchange.com/users/55516/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]
# [MATL](http://github.com/lmendo/MATL), ~~55~~ ~~53~~ ~~48~~ ~~43~~ 42 bytes
*Inspired by [@flawr's answer](https://codegolf.stackexchange.com/a/92254/36398).*
**Graphical output**:
```
0i:"Gto~+XytP*+t"t4=t1Y6Z+b+w~*]]tat3$)1YG
```
[Try it at MATL Online!](https://matl.io/?code=0i%3A%22Gto%7E%2BXytP%2a%2Bt%22t4%3Dt1Y6Z%2Bb%2Bw%7E%2a%5D%5Dtat3%24%291YG&inputs=30&version=19.1.0). 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](https://i.stack.imgur.com/eS2xS.png)](https://i.stack.imgur.com/eS2xS.png)
**ASCII output (43 bytes)**:
```
0i:"Gto~+XytP*+t"t4=t1Y6Z+b+w~*]]tat3$)48+c
```
[Try it online!](http://matl.tryitonline.net/#code=MGk6Ikd0b34rWHl0UCordCJ0ND10MVk2WitiK3d-Kl1ddGF0MyQpNDgrYw&input=MTY)
### 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]
# 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]
# 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
[![n=100](https://i.stack.imgur.com/k9WTV.png)](https://i.stack.imgur.com/k9WTV.png)
[Answer]
## 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]
# 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]
# 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]
# 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]
# 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]
# 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.
![figure](https://i.stack.imgur.com/Y4Nzo.png)
[Answer]
# [APL (Dyalog Unicode)](https://www.dyalog.com/), 51 bytes[SBCS](https://github.com/abrudz/SBCS)
```
{m⌿r/⍨m←⌈/×r←(4∘|+{+/⊢/4 2⍴⌊⍵÷4}⌺3 3)0,∘⍉∘⌽⍣4⊢⍵}⍣≡⍪
```
[Try it online!](https://tio.run/##hY/BSsNAEIbv@xS5rdKW7O5ED75NwESKKZHQi9RcpJQm7hY9iF6tl@JZEQoi2DeZF4n/RogRCl5mZ77/n5md@CIbnV7GWX7W8Op@nPPiVjUp4mzC9rMI2W0mqNguw91Dgewg4uXj1WA2CLleh1Fg2L2yrdm97d6jku2WAjpUQ5jYVT7aD3bPEcywlEi5emL30jRTvwSMq7Wfm7Lbnsg0HmcSGXjBN3OZn8tSCJqig0SEZ/W14cWdVFrJQGqlEX0u9FFPJC8SkRfJi8c9EXYvK6NNOwNDdK/@0YVIAn92fQ2gtAEDTIYdI6MM/WUESO2@jukWkNE91jb69n98@@bt3fv7P8Tu0OQb "APL (Dyalog Unicode) – Try It Online")
[Answer]
# 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!](https://tio.run/##NY7LasMwEAB/pRgK2mhFIhoK8kb2LfmAHI0pepDUhybBUYyEybe7Ek0PC7Mze1gnzs4ty6FTqi@DERMatOjQ05HlDWbOD13qu9g3H@2L9AZz5DLnAuIfMHH5AiGh3pDRcWfaWBuyOu1sm2pLTsfGZefI69T47Dw99yzAfLqOzGirlaIgBMGRbRVuFVAJSVvy60SJ69sj3FlVARQftSG3jgS3cbiEE6ve/VuFf49y3gM9lx8zXNgXTuC@zbhaTTTvmQnXgU2d7KFcLPLzFw "C (gcc) – Try It Online")
[Answer]
# 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`.
]