[Question]
  [
This question is about [abelian sandpiles](https://en.wikipedia.org/wiki/Abelian_sandpile_model). Read [this previous challenge](https://codegolf.stackexchange.com/questions/92251/build-a-sandpile) and [watch this numberphile video](https://www.youtube.com/watch?v=1MtEUErz7Gg) to learn more.
---
An abelian sandpile of size *n* by *n* is a grid containing the number 0, 1, 2 and 3 (representing the number of grains of sand). Adding two sandpiles works by first adding element by element, and then *toppling* any element that goes above 3. The order in which you topple doesn't matter, the end result is the same. When a cell topples its number decreases by 4, and each of its direct neighbors increases by 1. This can cause a chain reaction. If a cell is on the edge of the grid, any grains that fall off the grid while toppling disappear.
For example, I'm adding two 3 by 3 sandpiles (giving a rather extreme chain reaction):
```
3 3 3   1 2 1   4 5 4    4 6 4    6 2 6    6 3 6    2 5 2    4 1 4    4 2 4    0 4 0    2 0 2    2 1 2
3 3 3 + 2 1 2 = 5 4 5 -> 6 0 6 -> 2 4 2 -> 3 0 3 -> 5 0 5 -> 1 4 1 -> 2 0 2 -> 4 0 4 -> 0 4 0 -> 1 0 1
3 3 3   1 2 1   4 5 4    4 6 4    6 2 6    6 3 6    2 5 2    4 1 4    4 2 4    0 4 0    2 0 2    2 1 2
```
**In this challenge we're interested in a subset of all possible *n* by *n* sandpiles.** This subset contains any sandpile that you can get by adding an arbitrary sandpile to the all-3s *n* by *n* sandpile. For example, just above we saw that `212 | 101 | 212` is in the subset, because we got it by adding something to the all-3 sandpile.
Now this subset has an interesting element: the *identity* element. If you take this element and add it to *any other element in the **subset***, the sum is unchanged. In other words, this sandpile acts like a *zero* of this subset. It just so happens that `212 | 101 | 212` is the zero element for the subset of 3 by 3. For example:
```
2 2 2   2 1 2   4 3 4    0 5 0    2 1 2    2 2 2
2 2 2 + 1 0 1 = 3 2 3 -> 5 2 5 -> 1 6 1 -> 2 2 2
2 2 2   2 1 2   4 3 4    0 5 0    2 1 2    2 2 2
```
Now this is your challenge: **given *n*, find the identity element of the subset of the *n* by *n* grid**. Output it by assigning a unique color with sufficient contrast of your choice to each of `0, 1, 2, 3` and outputting an n by n image. **Your code must be able to produce the 50 by 50 case in under a minute on a reasonable modern PC.**
---
For example, the 500 by 500 identity element:
[![500 by 500 identity element](https://i.stack.imgur.com/q8Sqm.png)](https://i.stack.imgur.com/q8Sqm.png)
Here is blue = 3, green = 2, red = 1, white = 0. But you don't have to use this color scheme in your answer.
  
[Answer]
# Mathematica, ~~177~~ ~~157~~ ~~135~~ 133 bytes
```
Colorize[f=BlockMap[⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&,#~ArrayPad~1,{3,3},1]&~FixedPoint~#&;k=Table[6,#,#];f[k-f@k]]&
```
Takes a number `n`. The output is the identity sandpile. 0 is black, 1 is light gray, 2 is magenta, and 3 is blue-gray.
Sadly, Mathematica does not have a builtin for this...
Uses the algorithm stated in [Scott Corry and David Perkinson's paper](http://people.reed.edu/%7Edavidp/divisors_and_sandpiles/draft-11.20.2016.pdf).
Takes 91.7 seconds on my 5-year-old laptop to calculate the 50x50 identity sandpile. I am confident that a reasonable modern desktop computer is more than 50% faster. (I also have a way faster code at the end).
**Explanation**
```
f= ...
```
Define function `f` (the input is a sandpile matrix): a function that ...
```
BlockMap[ ... ]~FixedPoint~#&
```
... repeats the `BlockMap` operation until the output does not change. `BlockMap` operation: ...
---
```
#~ArrayPad~1
```
... pad the input array with one layer of 0 ...
```
{3,3},1
```
... partition it into 3x3 matrices, with offset 1 ...
```
⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&
```
... and for each partition, add the number of sand grains toppled onto the center cell and the center cell value mod 4.
i.e. the output of `f` is the stabilized version of the input.
---
```
k=Table[6,#,#]
```
Define `k` as an *n* by *n* array of 6s.
```
f[k-f@k]]
```
Compute f(k - f(k)).
```
Colorize[ ... ]
```
Apply colors to the result.
**Faster version (142 bytes)**
```
Colorize[r=RotateLeft;a=ArrayPad;f=a[b=#~a~1;b+r[g=⌊b/4⌋,s={0,1}]+g~r~-s+r[g,1-s]+r[g,s-1]-4g,-1]&~FixedPoint~#&;k=Table[6,#,#];f[k-f@k]]&
```
Same code, but uses built-in list rotation instead of `BlockMap`. Computes n=50 in 4.0 seconds on my laptop.
[Answer]
# Octave, ~~120~~ 113 bytes
```
function a=W(a);while nnz(b=a>3);a+=conv2(b,[t=[0 1 0];!t;t],'same')-b*4;end;end;@(n)imagesc(W((x=ones(n)*6)-W(x)))
```
Thanks to [JungHwan Min](https://codegolf.stackexchange.com/users/60043/junghwan-min) for providing a link to the reference paper in his Mathematica answer.  
Thanks to [Stewie Griffin](https://codegolf.stackexchange.com/users/31516/stewie-griffin) saved me 7 bytes `[any(any(x)) -> nnz(x)]`
Here two function is used:
1.`f`: for stabilization of a matrix  
2. An anonymous function that takes `n` as input and shows the identity matrix.
[Try It On rextester!](http://rextester.com/JVF23775) for generation of a 50 \* 50 matrix
Elapsed time for computation of the matrix: `0.0844409 seconds`.
**Explanation:**
Consider a function `f` that stabilizes a matrix the task of finding the identity is simply
`f(ones(n)*6 - f(ones(n)*6)`.
that `ones(n)*6` means a n\*n matrix of 6.
so for `n=3`:
```
M = [6 6 6
     6 6 6
     6 6 6];
```
The result will be `f(M-f(M))`
For stabilization function 2D convolution used to speed the task;
In each iteration we make a binary matrix `b` with the same size of the input matrix and set it to 1 if corresponding element of the input matrix is >3.
Then we apply a 2D convolution of the binary matrix with the following mask
```
0 1 0
1 0 1
0 1 0
```
representing four direct neighbors .  
Result of convolution is added to the matrix and 4 times the binary matrix subtracted from it.
The loop continued until all element of the matrix are <= 3
**Ungolfed version**:
```
function a=stabilize(a)
    mask = [t=[0 1 0];!t;t];
    while any(any(b=a>3))
        a+=conv2(b,mask,'same')-b*4;
    end
end
n= 50;
M = ones(n)*6;
result = stabilize(M-stabilize(M));
imagesc(result);
```
[Answer]
# Python 3 + Numpy + PIL, ~~385~~ ~~370~~ 364 bytes
```
import numpy as q,PIL.Image as w
n=int(input())
z=n,n
def r(p):
 while len(p[p>3]):
  for x,y in q.ndindex(z):
   if p[x,y]>3:
    p[x,y]-=4;p[x-1,y]+=x>0;p[x,y-1]+=y>0
    if~-n>x:p[x+1,y]+=1
    if~-n>y:p[x,y+1]+=1
s=q.full(z,6)
t=s.copy()
r(t)
i=s-t
r(i)
w.fromarray(q.uint8(q.array(q.vectorize(lambda x:[x//1*65]*3,otypes=[object])(i).tolist()))).save('i.png')
```
Takes input on STDIN. Outputs the image as greyscale to `i.png`. Black corresponds to 0, dark grey to 1, light grey to 2, and white to 0.
Uses the formula `I = R(S - R(S))`, where `I` is the identity element, `S` is the matrix filled with sixes, and `R` is the reduction function.
I could probably save some bytes by switching to Python 2 and doing `from numpy import*`, but (1) I don't have Numpy installed on Python 2 and (2) the program wasn't terminating with `from numpy import*`.
Ungolfed:
```
import numpy as np
from PIL import Image
# Compute the identity element
n = int(input('Size of the sandpile: '))
def reduce_pile(sandpile):
  while any(element >= 4 for element in np.nditer(sandpile)):
    for x, y in np.ndindex((n, n)):
      if sandpile[x, y] >= 4:
        sandpile[x, y] -= 4
        if x > 0: sandpile[x - 1, y] += 1
        if y > 0: sandpile[x, y - 1] += 1
        if x < n - 1: sandpile[x + 1, y] += 1
        if y < n - 1: sandpile[x, y + 1] += 1
s = np.full((n, n), 6, dtype=np.int32)
s_prime = np.copy(s)
reduce_pile(s_prime)
identity = s - s_prime
reduce_pile(identity)
# Output it to identity.png as an image
colours = [[255, 255, 255], [255, 0, 0], [0, 255, 0], [0, 0, 255]]
img_array = np.vectorize(lambda x: colours[x], otypes=[object])(identity)
img_array = np.array(img_array.tolist(), dtype=np.uint8)
img = Image.fromarray(img_array)
img.save('identity.png')
```
]