Abelian sandpile model

From Rosetta Code
Task
Abelian sandpile model
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Abelian sandpile model. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)


Implement the Abelian sandpile model also known as Bak–Tang–Wiesenfeld model. Its history, mathematical definition and properties can be found under its wikipedia article.

The task requires the creation of a 2D grid of arbitrary size on which "piles of sand" can be placed. Any "pile" that has 4 or more sand particles on it collapses, resulting in four particles being subtracted from the pile and distributed among its neighbors.

It is recommended to display the output in some kind of image format, as terminal emulators are usually too small to display images larger than a few dozen characters tall. As an example of how to accomplish this, see the Bitmap/Write a PPM file task.
Examples up to 2^30, wow!
javascript running on web
Examples:

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 4 0 0 -> 0 1 0 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0 0 0 0 0    0 0 0 0 0
0 0 0 0 0    0 0 1 0 0
0 0 6 0 0 -> 0 1 2 1 0
0 0 0 0 0    0 0 1 0 0
0 0 0 0 0    0 0 0 0 0

0  0 0  0  0    0 0 1 0 0
0  0 0  0  0    0 2 1 2 0
0  0 16 0  0 -> 1 1 0 1 1
0  0 0  0  0    0 2 1 2 0
0  0 0  0  0    0 0 1 0 0

11l

V grid = [[0] * 10] * 10
grid[5][5] = 64

print(‘Before:’)
L(row) grid
   print(row.map(c -> ‘#3’.format(c)).join(‘’))

F simulate(&grid)
   L
      V changed = 0B
      L(arr) grid
         V ii = L.index
         L(val) arr
            V jj = L.index
            I val > 3
               grid[ii][jj] -= 4
               I ii > 0
                  grid[ii - 1][jj]++
               I ii < grid.len - 1
                  grid[ii + 1][jj]++
               I jj > 0
                  grid[ii][jj - 1]++
               I jj < grid.len - 1
                  grid[ii][jj + 1]++
               changed = 1B
      I !changed
         L.break

simulate(&grid)

print("\nAfter:")
L(row) grid
   print(row.map(c -> ‘#3’.format(c)).join(‘’))

grid = [[0] * 65] * 65
grid[32][32] = 64 * 64

simulate(&grid)

V ppm = File(‘sand_pile.ppm’, ‘w’)
ppm.write_bytes(("P6\n#. #.\n255\n".format(grid.len, grid.len)).encode())
V colors = [[Byte(0),   0, 0],
            [Byte(255), 0, 0],
            [Byte(0), 255, 0],
            [Byte(0), 0, 255]]
L(row) grid
   L(c) row
      ppm.write_bytes(colors[c])
Output:
Before:
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0 64  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0

After:
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  1  2  1  0  0  0
  0  0  0  2  2  2  2  2  0  0
  0  0  1  2  2  2  2  2  1  0
  0  0  2  2  2  0  2  2  2  0
  0  0  1  2  2  2  2  2  1  0
  0  0  0  2  2  2  2  2  0  0
  0  0  0  0  1  2  1  0  0  0
  0  0  0  0  0  0  0  0  0  0

AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits or android 64 bits with application Termux
/* ARM assembly AARCH64 Raspberry PI 3B or android 64 bits */
/*  program abelian64.s   */ 

/* run : abelian 256 12 12  */

/*******************************************/
/* Constantes file                         */
/*******************************************/
/* for this file see task include a file in language AArch64 assembly*/
.include "../includeConstantesARM64.inc"
.equ MAXI, 25

/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessValue:        .asciz "@ "
szMessErrParam:     .asciz "error : command line = abelian size posx posy  \n"
szMessFin:          .asciz "End display :\n"
szCarriageReturn:   .asciz "\n"
 
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss
sZoneConv:        .skip 24
iSandPile:        .skip 8 * MAXI * MAXI
/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:                            // entry of program 
    mov fp,sp
    ldr x4,[fp]                  // load number of parameters command line
    cmp x4,#3                    // < 4 -> error
    ble 99f
    add x0,fp,#32                // load address param 4 = pos y
    ldr x0,[x0]
    bl conversionAtoD            // conversion ascii -> numeric
    mov x3,x0
    add x0,fp,#24                // load address param 3 = pos x
    ldr x0,[x0]
    bl conversionAtoD
    mov x2,x0
    add x0,fp,#16                 // load address param 2 = size begin pile
    ldr x0,[x0]
    bl conversionAtoD
    ldr x4,qAdriSandPile
    mov x5,#MAXI
    madd x5,x3,x5,x2              // compute offset = maxi * y + x
    str x0,[x4,x5,lsl #3]         // and store size in pos x,y
    //mov x0,x4                   // display start position
    //bl displaySandPile
    
    mov x0,x4                     // sandpile address
    mov x1,x2                     // pos x to start
    mov x2,x3                     // pos y to start
    bl addSand
    
    ldr x0,qAdrszMessFin
    bl affichageMess
    mov x0,x4
    bl displaySandPile
    b 100f
99:                               // line command error
   ldr x0,qAdrszMessErrParam
   bl affichageMess
100:                              // standard end of the program 
    mov x0,0                      // return code
    mov x8,EXIT                   // request to exit program
    svc 0                         // perform the system call
 
qAdrszCarriageReturn:     .quad szCarriageReturn
qAdrsZoneConv:            .quad sZoneConv
qAdrszMessErrParam:       .quad szMessErrParam
qAdrszMessFin:            .quad szMessFin
qAdriSandPile:            .quad iSandPile
/***************************************************/
/*     display  sandpile               */
/***************************************************/
// x0 contains address to sandpile
displaySandPile:
    stp x1,lr,[sp,-16]!         // save  registres
    stp x2,x3,[sp,-16]!         // save  registres
    stp x4,x5,[sp,-16]!         // save  registres
    stp x6,x7,[sp,-16]!         // save  registres
    mov x6,x0
    mov x3,#0                   // indice y
    mov x4,#MAXI
1:
    mov x2,#0                   // indice x
2:
    madd x5,x3,x4,x2            // compute offset
    ldr x0,[x6,x5,lsl #3]       // load value at pos x,y
    ldr x1,qAdrsZoneConv
    bl conversion10             // call decimal conversion
    add x1,x1,1
    mov x7,#0
    strb w7,[x1,x0]
    ldr x0,qAdrszMessValue
    ldr x1,qAdrsZoneConv        // insert value conversion in message
    bl strInsertAtCharInc
    bl affichageMess
    add x2,x2,1
    cmp x2,MAXI
    blt 2b
    ldr x0,qAdrszCarriageReturn
    bl affichageMess
    add x3,x3,1
    cmp x3,MAXI
    blt 1b

100:
    ldp x6,x7,[sp],16         // restaur des  2 registres
    ldp x4,x5,[sp],16         // restaur des  2 registres
    ldp x2,x3,[sp],16         // restaur des  2 registres
    ldp x1,lr,[sp],16         // restaur des  2 registres
    ret
qAdrszMessValue:       .quad szMessValue
/***************************************************/
/*     display  sandpile               */
/***************************************************/
// x0 contains address to sanspile
// x1 contains position x
// x2 contains position y
addSand:
    stp x1,lr,[sp,-16]!         // save  registres
    stp x2,x3,[sp,-16]!         // save  registres
    stp x4,x5,[sp,-16]!         // save  registres
    mov x3,#MAXI
    madd x4,x3,x2,x1            // compute offset
    ldr x5,[x0,x4,lsl #3]
1:
    cmp x5,#4                   // 4 grains ?
    blt 100f
    sub x5,x5,4                 // yes sustract
    str x5,[x0,x4,lsl #3]
    cmp x1,MAXI-1               // right position ok ?
    beq 2f
    add x1,x1,1                 // yes
    bl add1Sand                 // add 1 grain
    bl addSand                  // and compute new pile
    sub x1,x1,1
2:
    cmp x1,0                    // left position ok ?
    beq 3f
    sub x1,x1,1
    bl add1Sand
    bl addSand
    add x1,x1,1
3:
    cmp x2,0                    // higt position ok ?
    beq 4f
    sub x2,x2,1
    bl add1Sand
    bl addSand
    add x2,x2,1
4:
    cmp x2,MAXI-1               // low position ok ?
    beq 5f
    add x2,x2,1
    bl add1Sand
    bl addSand
    sub x2,x2,1
5:
   ldr x5,[x0,x4,lsl #3]       // reload value
   b 1b                        // and loop
100:
    ldp x4,x5,[sp],16         // restaur des  2 registres
    ldp x2,x3,[sp],16         // restaur des  2 registres
    ldp x1,lr,[sp],16         // restaur des  2 registres
    ret
/***************************************************/
/*     add 1 grain of sand              */
/***************************************************/
// x0 contains address to sanspile
// x1 contains position x
// x2 contains position y
add1Sand:
    stp x3,lr,[sp,-16]!       // save  registres
    stp x4,x5,[sp,-16]!       // save  registres
    mov x3,#MAXI
    madd x4,x3,x2,x1          // compute offset
    ldr x5,[x0,x4,lsl #3]     // load value at pos x,y
    add x5,x5,1
    str x5,[x0,x4,lsl #3]     // and store 
100:
    ldp x4,x5,[sp],16         // restaur des  2 registres
    ldp x3,lr,[sp],16         // restaur des  2 registres
    ret
/********************************************************/
/*        File Include fonctions                        */
/********************************************************/
/* for this file see task include a file in language AArch64 assembly */
.include "../includeARM64.inc"
Output:
~/.../rosetta/asm1 $ abelian64 64 12 12
End display :
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 2 2 2 0 2 2 2 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

ALGOL 68

Using code from the Abelian sandpile model/Identity task.

BEGIN # model Abelian sandpiles #
    # returns TRUE if the sandpile s is stable, FALSE otherwise #
    OP STABLE = ( [,]INT s )BOOL:
       BEGIN
           BOOL result := TRUE;
           FOR i FROM 1 LWB s TO 1 UPB s WHILE result DO
               FOR j FROM 2 LWB s TO 2 UPB s WHILE result := s[ i, j ] < 4 DO SKIP OD
           OD;
           result
       END # STABLE # ;
    # returns the sandpile s after avalanches #
    OP AVALANCHE = ( [,]INT s )[,]INT:
       BEGIN
           [ 1 : 1 UPB s, 1 : 2 UPB s ]INT result := s[ AT 1, AT 1 ];
           WHILE BOOL had avalanche := FALSE;
                 FOR i TO 1 UPB s DO
                     FOR j TO 2 UPB s DO
                         IF result[ i, j ] >= 4 THEN
                             # unstable pile #
                             had avalanche := TRUE;
                             result[ i, j ] -:= 4;
                             IF i > 1       THEN result[ i - 1, j     ] +:= 1 FI;
                             IF i < 1 UPB s THEN result[ i + 1, j     ] +:= 1 FI;
                             IF j > 1       THEN result[ i,     j - 1 ] +:= 1 FI;
                             IF j < 2 UPB s THEN result[ i,     j + 1 ] +:= 1 FI
                         FI
                     OD
                 OD;
                 had avalanche
           DO SKIP OD;
           result
       END # AVALANCHE # ;
    # returns the maximum element of s #
    OP   MAX = ( [,]INT s )INT:
         BEGIN
            INT result := s[ 1 LWB s, 2 LWB s ];
            FOR i FROM 1 LWB s TO 2 UPB s DO
                FOR j FROM 2 LWB s TO 2 UPB s DO
                    IF s[ i, j ] > result THEN result := s[ i, j ] FI
                OD
            OD;
            result
         END # MAX # ;
    # prints the sandpile s #
    PROC show sandpile = ( STRING title, [,]INT s )VOID:
    BEGIN
        print( ( title, newline ) );
        IF 1 UPB s >= 1 LWB s AND 2 UPB s >= 2 LWB s THEN
            # non-empty sandpile #
            INT width := 1; # find tthe width needed for each element #
            INT v     := MAX s;
            WHILE v > 9 DO
                v OVERAB 10;
                width +:= 1
            OD;
            FOR i TO 1 UPB s DO
                FOR j TO 2 UPB s DO
                    print( ( " ", whole( s[ i, j ], - width ) ) )
                OD;
                print( ( newline ) )
            OD
        FI
    END # show sandpile # ;
    # printys a sandpile before and after the avalanches #
    PROC show sandpile before and after = ( [,]INT s )VOID:
         BEGIN
             [ 1 LWB s : 1 UPB s, 2 LWB s : 2 UPB s ]INT t := s;
             show sandpile( "before: ", t );
             WHILE NOT STABLE t DO
                 t := AVALANCHE t
             OD;
             show sandpile( "after: ", t );
             print( ( newline ) )
         END # show sandpile before and after # ;
    # task test case #
    [,]INT s1 = ( (  0,  0,  0,  0,  0 )
                , (  0,  0,  0,  0,  0 )
                , (  0,  0, 16,  0,  0 )
                , (  0,  0,  0,  0,  0 )
                , (  0,  0,  0,  0,  0 )
                );
    show sandpile before and after( s1 );
    # test case from 11l, C, etc. #
    [ 1 : 10, 1 : 10 ]INT s2;
    FOR i FROM 1 LWB s2 TO 1 UPB s2 DO
        FOR j FROM 2 LWB s2 TO 2 UPB s2 DO
            s2[ i, j ] := 0
        OD
    OD;
    s2[ 6, 6 ] := 64;
    show sandpile before and after( s2 )
END
Output:
before:
  0  0  0  0  0
  0  0  0  0  0
  0  0 16  0  0
  0  0  0  0  0
  0  0  0  0  0
after:
 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

before:
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0 64  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
after:
 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 1 2 1 0 0 0
 0 0 0 2 2 2 2 2 0 0
 0 0 1 2 2 2 2 2 1 0
 0 0 2 2 2 0 2 2 2 0
 0 0 1 2 2 2 2 2 1 0
 0 0 0 2 2 2 2 2 0 0
 0 0 0 0 1 2 1 0 0 0
 0 0 0 0 0 0 0 0 0 0

ARM Assembly

Works with: as version Raspberry Pi
/* ARM assembly Raspberry PI  or android 32 bits */
/*  program abelian.s   */ 

/* run : abelian 256 12 12  */

/* REMARK 1 : this program use routines in a include file 
   see task Include a file language arm assembly 
   for the routine affichageMess conversion10 
   see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */
/************************************/
/* Constantes                       */
/************************************/
.include "../constantes.inc"
.equ MAXI, 25

/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessValue:        .asciz "@ "
szMessErrParam:     .asciz "error : command line = abelian size posx posy  \n"
szMessFin:          .asciz "End display :\n"
szCarriageReturn:   .asciz "\n"
 
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss
sZoneConv:        .skip 24
iSandPile:        .skip 4 * MAXI * MAXI
/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:                            @ entry of program 
    mov fp,sp
    ldr r4,[fp]                  @ load number of parameters commend line
    cmp r4,#3                    @ < 4 -> error
    ble 99f
    add r0,fp,#16                @ load address param 4 = pos y
    ldr r0,[r0]
    bl conversionAtoD            @ conversion ascii -> numeric
    mov r3,r0
    add r0,fp,#12                @ load address param 3 = pos x
    ldr r0,[r0]
    bl conversionAtoD
    mov r2,r0
    add r0,fp,#8                 @ load address param 2 = size begin pile
    ldr r0,[r0]
    bl conversionAtoD
    ldr r4,iAdriSandPile
    mov r5,#MAXI
    mul r5,r3,r5                 @ compute offset = maxi * y
    add r5,r2                    @ + x
    str r0,[r4,r5,lsl #2]        @ and store size in pos x,y
    //mov r0,r4                    @ display start position
    //bl displaySandPile
    
    mov r0,r4                   @ sandpile address
    mov r1,r2                   @ pos x to start
    mov r2,r3                   @ pos y to start
    bl addSand
    
    ldr r0,iAdrszMessFin
    bl affichageMess
    mov r0,r4
    bl displaySandPile
    b 100f
99:                                  @ line command error
   ldr r0,iAdrszMessErrParam
   bl affichageMess
100:                                  @ standard end of the program 
    mov r0, #0                        @ return code
    mov r7, #EXIT                     @ request to exit program
    svc #0                            @ perform the system call
 
iAdrszCarriageReturn:     .int szCarriageReturn
iAdrsZoneConv:            .int sZoneConv
iAdrszMessErrParam:       .int szMessErrParam
iAdrszMessFin:            .int szMessFin
iAdriSandPile:            .int iSandPile
/***************************************************/
/*     display  sandpile               */
/***************************************************/
// r0 contains address to sandpile
displaySandPile:
    push {r1-r6,lr}             @ save  registers 
    mov r6,r0
    mov r3,#0                   @ indice y
    mov r4,#MAXI
1:
    mov r2,#0                   @ indice x
2:
    mul r5,r3,r4
    add r5,r2                   @ compute offset
    ldr r0,[r6,r5,lsl #2]       @ load value at pos x,y
    ldr r1,iAdrsZoneConv
    bl conversion10             @ call decimal conversion
    add r1,#1
    mov r7,#0
    strb r7,[r1,r0]
    ldr r0,iAdrszMessValue
    ldr r1,iAdrsZoneConv        @ insert value conversion in message
    bl strInsertAtCharInc
    bl affichageMess
    add r2,#1
    cmp r2,#MAXI
    blt 2b
    ldr r0,iAdrszCarriageReturn
    bl affichageMess
    add r3,#1
    cmp r3,#MAXI
    blt 1b

100:
    pop {r1-r6,lr}             @ restaur registers
    bx lr                      @ return
iAdrszMessValue:       .int szMessValue
/***************************************************/
/*     display  sandpile               */
/***************************************************/
// r0 contains address to sanspile
// r1 contains position x
// r2 contains position y
addSand:
    push {r1-r5,lr}             @ save  registers 
    mov r3,#MAXI
    mul r4,r3,r2
    add r4,r1
    ldr r5,[r0,r4,lsl #2]
1:
    cmp r5,#4                   @ 4 grains ?
    blt 100f
    sub r5,#4                   @ yes sustract
    str r5,[r0,r4,lsl #2]
    cmp r1,#MAXI-1              @ right position ok ?
    beq 2f
    add r1,#1                   @ yes
    bl add1Sand                 @ add 1 grain
    bl addSand                  @ and compute new pile
    sub r1,#1
2:
    cmp r1,#0                   @ left position ok ?
    beq 3f
    sub r1,#1
    bl add1Sand
    bl addSand
    add r1,#1
3:
    cmp r2,#0                   @ higt position ok ?
    beq 4f
    sub r2,#1
    bl add1Sand
    bl addSand
    add r2,#1
4:
    cmp r2,#MAXI-1               @ low position ok ?
    beq 5f
    add r2,#1
    bl add1Sand
    bl addSand
    sub r2,#1
5:
   ldr r5,[r0,r4,lsl #2]       @ reload value
   b 1b                        @ and loop
100:
    pop {r1-r5,lr}             @ restaur registers
    bx lr                      @ return
/***************************************************/
/*     add 1 grain of sand              */
/***************************************************/
// r0 contains address to sanspile
// r1 contains position x
// r2 contains position y
add1Sand:
    push {r3-r5,lr}           @ save  registers 
    mov r3,#MAXI
    mul r4,r3,r2
    add r4,r1                 @ compute offset
    ldr r5,[r0,r4,lsl #2]     @ load value at pos x,y
    add r5,#1
    str r5,[r0,r4,lsl #2]     @ and store 
100:
    pop {r3-r5,lr}            @ restaur registers
    bx lr                     @ return
/***************************************************/
/*      ROUTINES INCLUDE                           */
/***************************************************/
.include "../affichage.inc"
Output:
pi@raspberrypi:~/asm32/rosetta32/ass10 $ abelian 512 12 12
End display :
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  1  3  0  2  2  2  0  3  1  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  2  1  2  3  2  3  2  3  2  1  2  0  0  0  0  0  0  0
0  0  0  0  0  0  2  3  3  0  2  0  3  0  2  0  3  3  2  0  0  0  0  0  0
0  0  0  0  0  2  3  0  3  3  3  3  3  3  3  3  3  0  3  2  0  0  0  0  0
0  0  0  0  1  1  3  3  0  3  1  3  3  3  1  3  0  3  3  1  1  0  0  0  0
0  0  0  0  3  2  0  3  3  2  2  0  3  0  2  2  3  3  0  2  3  0  0  0  0
0  0  0  1  0  3  2  3  1  2  2  2  3  2  2  2  1  3  2  3  0  1  0  0  0
0  0  0  1  2  2  0  3  3  0  2  0  3  0  2  0  3  3  0  2  2  1  0  0  0
0  0  0  1  2  3  3  3  3  3  3  3  0  3  3  3  3  3  3  3  2  1  0  0  0
0  0  0  1  2  2  0  3  3  0  2  0  3  0  2  0  3  3  0  2  2  1  0  0  0
0  0  0  1  0  3  2  3  1  2  2  2  3  2  2  2  1  3  2  3  0  1  0  0  0
0  0  0  0  3  2  0  3  3  2  2  0  3  0  2  2  3  3  0  2  3  0  0  0  0
0  0  0  0  1  1  3  3  0  3  1  3  3  3  1  3  0  3  3  1  1  0  0  0  0
0  0  0  0  0  2  3  0  3  3  3  3  3  3  3  3  3  0  3  2  0  0  0  0  0
0  0  0  0  0  0  2  3  3  0  2  0  3  0  2  0  3  3  2  0  0  0  0  0  0
0  0  0  0  0  0  0  2  1  2  3  2  3  2  3  2  1  2  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  1  3  0  2  2  2  0  3  1  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

AutoHotkey

Size := Size2 := 10
celula := [], Deltas := ["-1,0","1,1","0,-1","0,1"], Width := Size * 2.5
Gui, font, S%Size% 
Gui, add, text, y1
loop, 19
{
	Row := A_Index
	loop, 19
	{
		Col := A_Index
		Gui, add, button, % (Col=1 ? "xs y+1" : "x+1 yp") " v" Col "_" Row " w" Width " -TabStop"
		celula[Col,Row] := 0
		GuiControl, hide, %Col%_%Row%
	}
}
gui, color, black
Gui, show,, Abelian SandPile
InputBox, areia,, How much particles?,, 140, 140
celula[10,10] := areia
inicio:
loop 
{
GuiControl,, 10_10, % celula[10,10]
fim := true
loop 19
	{
		linha := A_Index
		loop 19
		{
			coluna := A_Index
			if (celula[coluna,linha] >= 4)
				{
					celula[coluna,linha] -= 4
					celula[coluna-1,linha] += 1
					celula[coluna+1,linha] += 1
					celula[coluna,linha-1] += 1
					celula[coluna,linha+1] += 1
					fim := false
				}
		}
	}
if (fim = true)
	break
}
loop 19
	{
		line := A_Index
		loop 19
			{
				column := A_Index
				GuiControlGet, posicao, Pos, % column "_" line
				switch celula[column,line]
					{
						case 0:
							gui, add, progress, Backgroundblue w24 h24 x%posicaoX% y%posicaoY%
						case 1:
							gui, add, progress, backgroundred w24 h24 x%posicaoX% y%posicaoY%
						case 2:
							gui, add, progress, backgroundgreen w24 h24 x%posicaoX% y%posicaoY%
						case 3:
							gui, add, progress, backgroundyellow w24 h24 x%posicaoX% y%posicaoY%
					}
			}
	}
msgbox sair
ExitApp
Image for 1000 particles grid 19x19

BQN

#!/usr/bin/env BQN
# Takes the size of grid (1 side) and the size of the starting pile as command line arguments
sizepile•BQN¨ 2↑•args
_Fp{𝔽𝔽_𝕣𝔽𝕩} # Fixed point function
Sand{p𝕩4  (𝕩+¯4×p)++´«,»,«˘,»˘{𝕎𝕩}¨<p}_Fp # Calculates given sand grid until it stops changing

lastSand pile˙((˜size÷2)) sizesize0 # Calculate the last state of the sand grid

# PPM output, code taken from https://rosettacode.org/wiki/Bitmap/Write_a_PPM_file#BQN
header_ppm"P6"(@+10)(•Repr size)" "(•Repr size)(@+10)"255"@+10
colors4255×=⌜˜3 # Black, Red, Green, Blue
image_ppm@+⥊colors˜4|last
"sandpile.ppm" •file.Bytes header_ppmimage_ppm

C

Writes out the initial and final sand piles to the console and the final sand pile to a PPM file.

#include<stdlib.h>
#include<string.h>
#include<stdio.h>

int main(int argc, char** argv)
{
	int i,j,sandPileEdge, centerPileHeight, processAgain = 1,top,down,left,right;	
	int** sandPile;
	char* fileName;
	static unsigned char colour[3];

	if(argc!=3){
		printf("Usage: %s <Sand pile side> <Center pile height>",argv[0]);
		return 0;
	}

	sandPileEdge = atoi(argv[1]);
	centerPileHeight = atoi(argv[2]);

	if(sandPileEdge<=0 || centerPileHeight<=0){
		printf("Sand pile and center pile dimensions must be positive integers.");
		return 0;
	}

	sandPile = (int**)malloc(sandPileEdge * sizeof(int*));

	for(i=0;i<sandPileEdge;i++){
		sandPile[i] = (int*)calloc(sandPileEdge,sizeof(int));
	}

	sandPile[sandPileEdge/2][sandPileEdge/2] = centerPileHeight;

	printf("Initial sand pile :\n\n");

	for(i=0;i<sandPileEdge;i++){
		for(j=0;j<sandPileEdge;j++){
			printf("%3d",sandPile[i][j]);
		}
		printf("\n");
	}

	while(processAgain == 1){

		processAgain = 0;
		top = 0;
		down = 0;
		left = 0;
		right = 0;

		for(i=0;i<sandPileEdge;i++){
			for(j=0;j<sandPileEdge;j++){
				if(sandPile[i][j]>=4){				
					if(i-1>=0){
						top = 1;
						sandPile[i-1][j]+=1;
						if(sandPile[i-1][j]>=4)
							processAgain = 1;
					}
					if(i+1<sandPileEdge){
						down = 1;
						sandPile[i+1][j]+=1;
						if(sandPile[i+1][j]>=4)
							processAgain = 1;
					}
					if(j-1>=0){
						left = 1;
						sandPile[i][j-1]+=1;
						if(sandPile[i][j-1]>=4)
							processAgain = 1;
					}
					if(j+1<sandPileEdge){
						right = 1;
						sandPile[i][j+1]+=1;
						if(sandPile[i][j+1]>=4)
							processAgain = 1;
					}
				sandPile[i][j] -= (top + down + left + right);
				if(sandPile[i][j]>=4)
					processAgain = 1;
				}
			}
		}
	}

	printf("Final sand pile : \n\n");

	for(i=0;i<sandPileEdge;i++){
		for(j=0;j<sandPileEdge;j++){
			printf("%3d",sandPile[i][j]);
		}
		printf("\n");
	}

	fileName = (char*)malloc((strlen(argv[1]) + strlen(argv[2]) + 23)*sizeof(char));

	strcpy(fileName,"Final_Sand_Pile_");
	strcat(fileName,argv[1]);
	strcat(fileName,"_");
	strcat(fileName,argv[2]);
	strcat(fileName,".ppm");
	
	FILE *fp = fopen(fileName,"wb");

	fprintf(fp,"P6\n%d %d\n255\n",sandPileEdge,sandPileEdge);

	for(i=0;i<sandPileEdge;i++){
		for(j=0;j<sandPileEdge;j++){
			colour[0] = (sandPile[i][j] + i)%256;
			colour[1] = (sandPile[i][j] + j)%256;
			colour[2] = (sandPile[i][j] + i*j)%256;
			fwrite(colour,1,3,fp);
		}
	}
	
	fclose(fp);

	printf("\nImage file written to %s\n",fileName);

	return 0;
}

Console output :

abhishek_ghosh@Azure:~/doodles$ ./a.out 10 64
Initial sand pile :

  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0 64  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
Final sand pile :

  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0  0
  0  0  0  0  1  2  1  0  0  0
  0  0  0  2  2  2  2  2  0  0
  0  0  1  2  2  2  2  2  1  0
  0  0  2  2  2  0  2  2  2  0
  0  0  1  2  2  2  2  2  1  0
  0  0  0  2  2  2  2  2  0  0
  0  0  0  0  1  2  1  0  0  0
  0  0  0  0  0  0  0  0  0  0

Image file written to Final_Sand_Pile_10_64.ppm

C++

Works with: g++ version 9.2.0 20061115
Library: xtensor
Library: xtensor-io


#include <iostream>
#include "xtensor/xarray.hpp"
#include "xtensor/xio.hpp"
#include "xtensor-io/ximage.hpp"

xt::xarray<int> init_grid (unsigned long x_dim, unsigned long y_dim)
{
    xt::xarray<int>::shape_type shape = { x_dim, y_dim };
    xt::xarray<int> grid(shape);

    grid(x_dim/2, y_dim/2) = 64000;
    
    return grid;
}

int print_grid(xt::xarray<int>& grid)
{
    // for output to the terminal uncomment next line
    // only makes sense for small grid < 32x32;
    // std::cout << grid << std::endl << std::endl;

    // output result to an image
    xt::dump_image("grid.jpg", grid);

    return 0;
}

bool iterate_grid(xt::xarray<int>& grid, const unsigned long& x_dim, const unsigned long& y_dim)
{
    bool changed = false;

    for (unsigned long i=0; i < x_dim; ++i)
    {
        for (unsigned long j=0; j < y_dim; ++j)
        {
            if ( grid(i, j) >= 4 )
            {
                grid(i, j) -= 4;
                changed = true;
                try
                {
                    grid.at(i-1, j) += 1;
                    grid.at(i+1, j) += 1;
                    grid.at(i, j-1) += 1;
                    grid.at(i, j+1) += 1;
                }
                catch (const std::out_of_range& oor)
                {
                }
            }
        }
    }

    return changed;
}

int main(int argc, char* argv[])
{
    const unsigned long x_dim { 200 };
    const unsigned long y_dim { 200 };

    xt::xarray<int> grid = init_grid(x_dim, y_dim);
    bool changed { true };

    iterate_grid(grid, x_dim, y_dim);

    while (changed == true)
    {
        changed = iterate_grid(grid, x_dim, y_dim);
    }
    print_grid(grid);

    return 0;
}

Compile with following CMakeList.txt:

cmake_minimum_required(VERSION 3.1)
project(abelian_sandpile)

find_package(xtl REQUIRED)
find_package(xtensor REQUIRED)
# if xtensor was built with xsimd support:
# find_package(xsimd REQUIRED)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
include_directories(/usr/include/OpenImageIO)
find_library(OIIO "OpenImageIO")

add_executable(abelian_sandpile src/abelian_sandpile.cpp)

target_compile_options(abelian_sandpile PRIVATE -march=native -std=c++14)
target_link_libraries(abelian_sandpile xtensor ${OIIO})

Delphi

Translation of: Python
program Abelian_sandpile_model;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  System.SysUtils,
  Vcl.Graphics,
  System.Classes;

type
  TGrid = array of array of Integer;

function Iterate(var Grid: TGrid): Boolean;
var
  changed: Boolean;
  i: Integer;
  j: Integer;
  val: Integer;
  Alength: Integer;
begin
  Alength := length(Grid);
  changed := False;

  for i := 0 to High(Grid) do
    for j := 0 to High(Grid[0]) do
    begin
      val := Grid[i, j];
      if val > 3 then
      begin
        Grid[i, j] := Grid[i, j] - 4;

        if i > 0 then
          Grid[i - 1, j] := Grid[i - 1, j] + 1;

        if i < Alength - 1 then
          Grid[i + 1, j] := Grid[i + 1, j] + 1;

        if j > 0 then
          Grid[i, j - 1] := Grid[i, j - 1] + 1;

        if j < Alength - 1 then
          Grid[i, j + 1] := Grid[i, j + 1] + 1;
        changed := True;
      end;
    end;
  Result := changed;
end;

procedure Simulate(var Grid: TGrid);
var
  changed: Boolean;
begin
  while Iterate(Grid) do
    ;
end;

procedure Zeros(var Grid: TGrid; Size: Integer);
var
  i, j: Integer;
begin
  SetLength(Grid, Size, Size);
  for i := 0 to Size - 1 do
    for j := 0 to Size - 1 do
      Grid[i, j] := 0;
end;

procedure Println(Grid: TGrid);
var
  i, j: Integer;
begin
  for i := 0 to High(Grid) do
  begin
    Writeln;
    for j := 0 to High(Grid[0]) do
      Write(Format('%3d', [Grid[i, j]]));
  end;
  Writeln;
end;

function Grid2Bmp(Grid: TGrid): TBitmap;
const
  Colors: array[0..2] of TColor = (clRed, clLime, clBlue);
var
  Alength: Integer;
  i: Integer;
  j: Integer;
begin
  Alength := Length(Grid);

  Result := TBitmap.Create;
  Result.SetSize(Alength, Alength);

  for i := 0 to Alength - 1 do
    for j := 0 to Alength - 1 do
    begin
      Result.Canvas.Pixels[i, j] := Colors[Grid[i, j]];
    end;
end;

procedure Grid2P6(Grid: TGrid; FileName: TFileName);
var
  f: text;
  i, j, Alength: Integer;
  ppm: TFileStream;
  Header: AnsiString;
const
  COLORS: array[0..3] of array[0..2] of byte =
 //  R,   G,    B
((0   ,   0,    0),
 (255 ,   0,    0),
 (0   , 255,   0),
 (0   ,   0, 255));
begin
  Alength := Length(Grid);
  ppm := TFileStream.Create(FileName, fmCreate);
  Header := Format('P6'#10'%d %d'#10'255'#10, [Alength, Alength]);
  writeln(Header);
  ppm.Write(Tbytes(Header), Length(Header));

  for i := 0 to Alength - 1 do
    for j := 0 to Alength - 1 do
    begin
      ppm.Write(COLORS[Grid[i, j]], 3);
    end;
  ppm.Free;
end;

const
  DIMENSION = 10;

var
  Grid: TGrid;
  bmp: TBitmap;

begin
  Zeros(Grid, DIMENSION);
  Grid[4, 4] := 64;
  Writeln('Before:');
  Println(Grid);

  Simulate(Grid);

  Writeln(#10'After:');
  Println(Grid);

  // Output bmp
  with Grid2Bmp(Grid) do
  begin
    SaveToFile('output.bmp');
    free;
  end;

  // Output ppm
  Grid2P6(Grid, 'output.ppm');

  Readln;
end.


Forth

Works with: gforth version 0.7.3


#! /usr/bin/gforth -d 20M
\ Abelian Sandpile Model

0 assert-level !

\ command-line

: parse-number  s>number? invert throw drop ;
: parse-size    ." size  : " next-arg parse-number dup . cr ;
: parse-height  ." height: " next-arg parse-number dup . cr ;
: parse-args    cr parse-size parse-height ;

parse-args constant HEIGHT constant SIZE

: allot-erase   create here >r dup allot r> swap erase ;
: size^2        SIZE dup * cells ;
: 2cells        [ 2 cells ] literal ;
: -2cells       [ 2cells negate ] literal ;

size^2 allot-erase arr

\ array processing
: ix            swap SIZE * + cells arr + ;
: center        SIZE 2/ dup ;
: write-cell    ix @ u. ;
: write-row     SIZE 0 ?do dup i write-cell loop drop cr ;
: arr.          SIZE 0 ?do i write-row loop ;

\ stack processing

: stack-empty?  dup -1 = ;
: stack-full?   stack-empty? invert ;

\ pgm-handling

: concat        { a1 l1 a2 l2 } l1 l2 + allocate throw dup dup a1 swap l1 cmove a2 swap l1 + l2 cmove l1 l2 + ;
: write-pgm     ." P2" cr SIZE u. SIZE u. cr ." 3" cr arr. ;
: u>s           0 <# #s #> ;
: filename      s" sandpile-" SIZE u>s concat s" -" concat HEIGHT u>s concat s" .pgm" concat ;
: to-pgm        filename w/o create-file throw ['] write-pgm over outfile-execute close-file throw ;

\ sandpile

: prep-arr      HEIGHT center ix ! ;
: prep-stack    -1 HEIGHT 4 u>= if center then ;
: prepare       prep-arr prep-stack ;
: ensure        if else 2drop 0 2rdrop exit then ;
: col>=0        dup 0>= ensure ;
: col<SIZE      dup SIZE < ensure ;
: row>=0        over 0>= ensure ;
: row<SIZE      over SIZE < ensure ;
: legal?        col>=0 col<SIZE row>=0 row<SIZE 2drop true ;
: north         1. d- ;
: east          1+ ;
: south         1. d+ ;
: west          1- ;
: reduce        2dup ix dup -4 swap +! @ 4 < if 2drop then ;
: increase      2dup legal? if 2dup ix dup 1 swap +! @ 4 = if 2swap else 2drop then else 2drop then ; 
: inc-north     2dup north increase ;
: inc-east      2dup east increase ;
: inc-south     2dup south increase ;
: inc-west      2dup west increase ;
: inc-all       inc-north inc-east inc-south inc-west 2drop ;
: simulate      prepare begin stack-full? while 2dup 2>r reduce 2r> inc-all repeat drop to-pgm ." written to " filename type cr ;

simulate bye
Output:

sandpile with 5000 grains of sand: ./sandpile.fs 61 5000: [1]
sandpile with 50000 grains of sand: ./sandpile.fs 201 50000: [2]
sandpile with 500000 grains of sand: ./sandpile.fs 601 500000: [3]

Fortran

Works with: gfortran version 9.2.0

The Abelian sandpile operations are defined here.

module abelian_sandpile_m

  implicit none

  private
  public :: pile

  type :: pile
    !! usage:
    !!    1) init
    !!    2) run

    integer, allocatable :: grid(:,:)
    integer              :: n(2)

  contains
    procedure :: init
    procedure :: run

    procedure, private :: process_node
    procedure, private :: inside
  end type

contains

  logical function inside(this, i)
    class(pile), intent(in) :: this
    integer,     intent(in) :: i(2)

    inside = ((i(1) > 0) .and. (i(1) <= this%n(1)) .and. (i(2) > 0) .and. (i(2) <= this%n(2)) )
  end function

  recursive subroutine process_node(this, i)
    !! start process

    class(pile), intent(inout) :: this
    integer,     intent(in)    :: i(2)
      !! node coordinates to process

    integer :: i0(2,2), j(2), d, k

    ! if node has more than 4 grains -> redistribute
    if (this%grid(i(1),i(2)) >= 4) then
      ! unit vectors: help shift only one dimension (see below)
      i0 = reshape([1,0,0,1], [2,2])

      ! subtract 4 grains
      this%grid(i(1),i(2)) = this%grid(i(1),i(2))-4

      ! add one grain to neighbor if not out of bound
      do d = 1, 2               ! loop dimensions
        do k = -1, 1, 2         ! loop +-1 step in direction d
          j = i+k*i0(:,d)       ! j = i, but one element is shifted by +-1
          if (this%inside(j)) this%grid(j(1),j(2)) = this%grid(j(1),j(2)) + 1
        end do
      end do

      ! check neighbor nodes
      do d = 1, 2               ! loop dimensions
        do k = -1, 1, 2         ! loop +-1 step in direction d
          j = i+k*i0(:,d)       ! j = i, but one element is shifted by +-1
          if (this%inside(j)) call this%process_node(j)
        end do
      end do

      ! check itself
      call this%process_node(i)
    end if
  end subroutine

  subroutine run(this)
    !! start process

    class(pile), intent(inout) :: this

    ! only node that could be unstable is inital node
    call this%process_node(this%n/2)
  end subroutine

  subroutine init(this, nx, ny, h)
    class(pile), intent(out) :: this
    integer,     intent(in)  :: nx, ny
      !! grid dimensions
    integer,     intent(in)  :: h
      !! height of and grains in middle of grid

    this%n = [nx, ny]
    allocate (this%grid(nx,ny), source=0)
    this%grid(nx/2, ny/2) = h
  end subroutine

end module

The main program calls the abelian_sandpile_m and creates an ppm bitmap file by loading rgbimage_m module, which is defined here.

program main

  use rgbimage_m
  use abelian_sandpile_m

  implicit none

  integer :: nx, ny, i, j

  integer :: colors(0:3,3)

  type(rgbimage) :: im
  type(pile) :: p

  colors(0,:) = [255,255,255]
  colors(1,:) = [0,0,90]
  colors(2,:) = [0,0,170]
  colors(3,:) = [0,0,255]

  nx = 200
  ny = 100

  call p%init(nx, ny, 2000)
  call p%run

  call im%init(nx, ny)

  do i = 1, nx
    do j = 1, ny
      call im%set_pixel(i, j, colors(p%grid(i,j),:))
    end do
  end do

  call im%write('fig.ppm')

end program

Fōrmulæ

Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.

Programs in Fōrmulæ are created/edited online in its website.

In this page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.

F#

// Abelian sandpile model. Nigel Galloway: July 20th., 2020
type Sandpile(x,y,N:int[])=
  member private this.x=x
  member private this.y=y
  member private this.i=let rec topple n=match Array.tryFindIndex(fun n->n>3)n with
                                          None->n
                                         |Some g->let i=n.[g]/4
                                                  n.[g]<-n.[g]%4
                                                  match g%x,g/x with
                                                   (0,0)->n.[x]<-n.[x]+i;n.[1]<-n.[1]+i;topple n
                                                  |(α,0) when α=x-1->n.[g+x]<-n.[g+x]+i;n.[g-1]<-n.[g-1]+i;topple n
                                                  |(_,0)->n.[g-1]<-n.[g-1]+i;n.[g+1]<-n.[g+1]+i;n.[g+x]<-n.[g+x]+i;topple n
                                                  |(0,β) when β=y-1->n.[g-x]<-n.[g-x]+i;n.[g+1]<-n.[g+1]+i;topple n
                                                  |(0,β)->n.[g-x]<-n.[g-x]+i;n.[g+1]<-n.[g+1]+i;n.[g+x]<-n.[g+x]+i;topple n
                                                  |(α,β) when α=x-1 && β=y-1->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;topple n
                                                  |(α,_) when α=x-1->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;n.[g+x]<-n.[g+x]+i;topple n
                                                  |(_,β) when β=y-1->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;n.[g+1]<-n.[g+1]+i;topple n
                                                  |_->n.[g-1]<-n.[g-1]+i;n.[g-x]<-n.[g-x]+i;n.[g+x]<-n.[g+x]+i;n.[g+1]<-n.[g+1]+i;topple n
                        topple N
  static member (+) (n:Sandpile, g:Sandpile)=Sandpile(n.x,n.y,Array.map2(fun n g->n+g) n.i g.i)
  member this.toS=sprintf "%A" (this.i|>Array.chunkBySize x|>array2D)

printfn "%s\n" (Sandpile(3,3,[|4;3;3;3;1;2;0;2;3|])).toS
let e1=Array.zeroCreate<int> 25 in e1.[12]<-4; printfn "%s\n" (Sandpile(5,5,e1)).toS
let e1=Array.zeroCreate<int> 25 in e1.[12]<-6; printfn "%s\n" (Sandpile(5,5,e1)).toS
let e1=Array.zeroCreate<int> 25 in e1.[12]<-16; printfn "%s\n" (Sandpile(5,5,e1)).toS
Output:
[[2; 1; 0]
 [0; 3; 3]
 [1; 2; 3]]

[[0; 0; 0; 0; 0]
 [0; 0; 1; 0; 0]
 [0; 1; 0; 1; 0]
 [0; 0; 1; 0; 0]
 [0; 0; 0; 0; 0]]

[[0; 0; 0; 0; 0]
 [0; 0; 1; 0; 0]
 [0; 1; 2; 1; 0]
 [0; 0; 1; 0; 0]
 [0; 0; 0; 0; 0]]

[[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]]


FreeBASIC

ScreenRes 320, 200, 8
WindowTitle "Abelian sandpile model"

Dim As Long dimen = 220
Dim As Long pila1(dimen*dimen), pila2(dimen*dimen)
Dim As Long i, x, y
Dim As Long partic = 400000
Dim As Double t0 = Timer
Do 
    i = 0
    For y = 0 To dimen-1
        For x = 0 To dimen-1
            If x = dimen/2 And y = dimen/2 Then
                partic -= 4
                pila1(i) += 4
            End If
            If pila1(i) >= 4 Then
                pila1(i) -= 4
                pila2(i-1) += 1
                pila2(i+1) += 1
                pila2(i-dimen) += 1
                pila2(i+dimen) += 1
            End If
            Pset(x, y), (pila1(i)*2)
            i += 1
            Swap pila1(i), pila2(i)
        Next x
    Next y
Loop Until partic < 4
Bsave "abelian_sandpile.bmp",0
Sleep
Output:
[https://www.dropbox.com/s/lky2ji9n15gn5u6/abelian_sandpile.bmp?dl=0 Abelian sandpile model FreeBasic image]

FutureBasic

_mFile = 1
begin enum
  _iClose
end enum

_window = 1
begin enum 1
  _gridView
  _gridSizeLabel
  _gridSizeFld
  _centerNumLabel
  _centerNumFld
  _colorRadio
  _monoRadio
  _avalancheBtn
end enum


void local fn BuildMenu
  menu _mFile,,, @"File"
  menu _mFile, _iClose,, @"Close", @"w"
  MenuItemSetAction( _mFile, _iClose, @"performClose:" )
end fn


void local fn BuildWindow
  window _window, @"Abelian Sandpile Model", (0,0,513,360), NSWindowStyleMaskTitled + NSWindowStyleMaskClosable + NSWindowStyleMaskMiniaturizable
  
  subclass view _gridView, (20,20,320,320)
  
  textlabel _gridSizeLabel, @"Grid size:", (385,322,61,16)
  textfield _gridSizeFld,, @"5", (452,319,41,21)
  ControlSetFormat( _gridSizeFld, @"0123456789", YES, 4, 0 )
  
  textlabel _centerNumLabel, @"Center number:", (347,285,99,16)
  textfield _centerNumFld,, @"32", (452,282,41,21)
  ControlSetFormat( _centerNumFld, @"0123456789", YES, 4, 0 )
  
  radiobutton _colorRadio,, NSControlStateValueOn, @"Color", (367,249,59,18)
  radiobutton _monoRadio,, NSControlStateValueOff, @"Mono", (432,249,61,18)
  
  button _avalancheBtn,,, @"Avalanche", (375,13,96,32)
  
  WindowMakeFirstResponder( _window, _gridSizeFld )
end fn


void local fn ViewDrawRect
  long gridSize = fn ControlIntegerValue(_gridSizeFld)
  CGRect bounds = fn ViewBounds( _gridView )
  CGFloat cellSize = bounds.size.width/gridSize
  ColorRef col0 = fn ColorWhite, col1, col2, col3
  long r, c, value
  CGFloat x = 0, y = 0
  ColorRef color
  
  if ( fn ButtonState( _colorRadio ) == NSControlStateValueOn )
    col1 = fn ColorRed
    col2 = fn ColorGreen
    col3 = fn ColorBlue
  else
    col1 = fn ColorWithRGB( 0.25, 0.25, 0.25, 1.0 )
    col2 = fn ColorWithRGB( 0.5, 0.5, 0.5, 1.0 )
    col3 = fn ColorWithRGB( 0.75, 0.75, 0.75, 1.0 )
  end if
  
  for r = 0 to gridSize-1
    for c = 0 to gridSize-1
      value = mda_integer(r,c)
      select ( value )
        case 1    : color = col1
        case 2    : color = col2
        case 3    : color = col3
        case else : color = col0
      end select
      
      BezierPathFillRect( fn CGRectMake( x, y, cellSize, cellSize ), color )
      x += cellSize
    next
    x = 0
    y += cellSize
  next
end fn


void local fn AvalancheAction
  long r, c, gridSize = fn ControlIntegerValue(_gridSizeFld)
  long centerNum = fn ControlIntegerValue(_centerNumFld)
  long midNum = gridSize/2
  long limit = gridSize-1
  BOOL stable = NO
  long value
  
  // initialize array
  mda_kill
  for r = 0 to gridSize-1
    for c = 0 to gridSize-1
      mda(r,c) = 0
    next
  next
  mda(midNum,midNum) = centerNum
  
  // collapse
  while ( stable == NO )
    stable = YES
    for r = 0 to gridSize-1
      for c = 0 to gridSize-1
        value = mda_integer(r,c)
        if ( value > 3 )
          mda(r,c) = @(mda_integer(r,c)-4)
          if ( r > 0 ) then mda(r-1,c) = @(mda_integer(r-1,c) + 1)
          if ( r < limit ) then mda(r+1,c) = @(mda_integer(r+1,c) + 1)
          if ( c > 0 ) then mda(r,c-1) = @(mda_integer(r,c-1) + 1)
          if ( c < limit ) then mda(r,c+1) = @(mda_integer(r,c+1) + 1)
          stable = NO : break
        end if
      next
      if ( stable == NO ) then break
    next
  wend
  
  ViewSetNeedsDisplay( _gridView )
end fn


void local fn DoAppEvent( ev as long )
  select ( ev )
    case _appWillFinishLaunching
      fn BuildMenu
      fn BuildWindow
      fn AvalancheAction
    case _appShouldTerminateAfterLastWindowClosed : AppEventSetBool(YES)
  end select
end fn

void local fn DoDialog( ev as long, tag as long, wnd as long, obj as CFTypeRef )
  select ( ev )
    case _btnClick
      select ( tag )
        case _avalancheBtn : fn AvalancheAction
        case _gridSizeFld, _centerNumFld  : fn AvalancheAction
        case _colorRadio, _monoRadio : ViewSetNeedsDisplay( _gridView )
      end select
    case _viewDrawRect : fn ViewDrawRect
  end select
end fn

on appevent fn DoAppEvent
on dialog fn DoDialog

HandleEvents

Go

Translation of: Rust


Stack management in Go is automatic, starting very small (2KB) for each goroutine and expanding as necessary until the maximum allowed size is reached.

package main

import (
    "fmt"
    "log"
    "os"
    "strings"
)

const dim = 16 // image size

func check(err error) {
    if err != nil {
        log.Fatal(err)
    }
}

// Outputs the result to the terminal using UTF-8 block characters.
func drawPile(pile [][]uint) {
    chars:= []rune(" ░▓█")
    for _, row := range pile {
        line := make([]rune, len(row))
        for i, elem := range row {
            if elem > 3 { // only possible when algorithm not yet completed.
                elem = 3
            }
            line[i] = chars[elem]
        }
        fmt.Println(string(line))
    }
}

// Creates a .ppm file in the current directory, which contains
// a colored image of the pile.
func writePile(pile [][]uint) {
    file, err := os.Create("output.ppm")
    check(err)
    defer file.Close()
    // Write the signature, image dimensions and maximum color value to the file.
    fmt.Fprintf(file, "P3\n%d %d\n255\n", dim, dim)
    bcolors := []string{"125 0 25 ", "125 80 0 ", "186 118 0 ", "224 142 0 "}
    var line strings.Builder
    for _, row := range pile {        
        for _, elem := range row {
            line.WriteString(bcolors[elem])
        }
        file.WriteString(line.String() + "\n")
        line.Reset() 
    }
}

// Main part of the algorithm, a simple, recursive implementation of the model.
func handlePile(x, y uint, pile [][]uint) {
    if pile[y][x] >= 4 {
        pile[y][x] -= 4
        // Check each neighbor, whether they have enough "sand" to collapse and if they do,
        // recursively call handlePile on them.
        if y > 0 {
            pile[y-1][x]++
            if pile[y-1][x] >= 4 {
                handlePile(x, y-1, pile)
            }
        }
        if x > 0 {
            pile[y][x-1]++
            if pile[y][x-1] >= 4 {
                handlePile(x-1, y, pile)
            }
        }
        if y < dim-1 {
            pile[y+1][x]++
            if pile[y+1][x] >= 4 {
                handlePile(x, y+1, pile)
            }
        }
        if x < dim-1 {
            pile[y][x+1]++
            if pile[y][x+1] >= 4 {
                handlePile(x+1, y, pile)
            }
        }

        // Uncomment this line to show every iteration of the program.
        // Not recommended with large input values.
        // drawPile(pile)

        // Finally call the function on the current cell again,
        // in case it had more than 4 particles.
        handlePile(x, y, pile)
    }
}

func main() {
    // Create 2D grid and set size using the 'dim' constant.
    pile := make([][]uint, dim)
    for i := 0; i < dim; i++ {
        pile[i] = make([]uint, dim)
    }

    // Place some sand particles in the center of the grid and start the algorithm.
    hdim := uint(dim/2 - 1)
    pile[hdim][hdim] = 16
    handlePile(hdim, hdim, pile)
    drawPile(pile)

    // Uncomment this to save the final image to a file
    // after the recursive algorithm has ended.
    // writePile(pile)
}
Output:
                
                
                
                
                
       ░        
      ▓░▓       
     ░░ ░░      
      ▓░▓       
       ░        
                
                
                
                
                
                
       

Haskell

Works with: GHC version 8.8.1
Library: base version 4.13.0.0
Library: array version 0.5.4.0
Library: mtl version 2.2.2


Using a custom monad to make the code cleaner.

{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables        #-}

module Rosetta.AbelianSandpileModel.ST 
    ( simulate
    , test
    , toPGM
    ) where

import Control.Monad.Reader (asks, MonadReader (..), ReaderT, runReaderT)
import Control.Monad.ST (runST, ST)
import Control.Monad.State (evalStateT, forM_, lift, MonadState (..), StateT, modify, when)
import Data.Array.ST (freeze, readArray, STUArray, thaw, writeArray)
import Data.Array.Unboxed (array, assocs, bounds, UArray, (!))
import Data.Word (Word32)
import System.IO (hPutStr, hPutStrLn, IOMode (WriteMode), withFile)
import Text.Printf (printf)

type Point     = (Int, Int)
type ArrayST s = STUArray s Point Word32
type ArrayU    = UArray Point Word32

newtype M s a = M (ReaderT (S s) (StateT [Point] (ST s)) a)
    deriving (Functor, Applicative, Monad, MonadReader (S s), MonadState [Point])

data S s = S 
    { bMin :: !Point
    , bMax :: !Point
    , arr  :: !(ArrayST s)
    }

runM :: M s a -> S s -> [Point]-> ST s a
runM (M m) = evalStateT . runReaderT m

liftST :: ST s a -> M s a
liftST = M . lift . lift

simulate :: ArrayU -> ArrayU
simulate a = runST $ simulateST a

simulateST :: forall s. ArrayU -> ST s ArrayU
simulateST a = do
    let (p1, p2) = bounds a
        s = [p | (p, c) <- assocs a, c >= 4]
    b <- thaw a :: ST s (ArrayST s)
    let st = S { bMin = p1
               , bMax = p2
               , arr  = b
               }
    runM simulateM st s

simulateM :: forall s. M s ArrayU
simulateM = do
    ps <- get
    case ps of
        []      -> asks arr >>= liftST . freeze
        p : ps' -> do
            c <- changeArr p $ \x -> x - 4
            when (c < 4) $ put ps'
            forM_ [north, east, south, west] $ inc . ($ p)
            simulateM

changeArr :: Point -> (Word32 -> Word32) -> M s Word32
changeArr p f = do
    a    <- asks arr
    oldC <- liftST $ readArray a p
    let newC = f oldC
    liftST $ writeArray a p newC
    return newC

inc :: Point -> M s ()
inc p = do
    b <- inBounds p
    when b $ do
        c <- changeArr p succ
        when (c == 4) $ modify $ (p :)

inBounds :: Point -> M s Bool
inBounds p = do
    st <- ask
    return $ p >= bMin st && p <= bMax st

north, east, south, west :: Point -> Point
north (x, y) = (x, y + 1)
east  (x, y) = (x + 1, y)
south (x, y) = (x, y - 1)
west  (x, y) = (x - 1, y)

toPGM :: ArrayU -> FilePath -> IO ()
toPGM a fp = withFile fp WriteMode $ \h -> do
    let ((x1, y1), (x2, y2)) = bounds a
        width  = x2 - x1 + 1
        height = y2 - y1 + 1
    hPutStrLn h "P2"
    hPutStrLn h $ show width ++ " " ++ show height
    hPutStrLn h "3"
    forM_ [y1 .. y2] $ \y -> do
        forM_ [x1 .. x2] $ \x -> do
            let c = min 3 $ a ! (x, y)
            hPutStr h $ show c ++ " "
        hPutStrLn h ""

initArray :: Int -> Word32 -> ArrayU
initArray size height = array 
    ((-size, -size), (size, size))
    [((x, y), if x == 0 && y == 0 then height else 0) | x <- [-size .. size], y <- [-size .. size]]

test :: Int -> Word32 -> IO ()
test size height = do
    printf "size = %d, height = %d\n" size height
    let a  = initArray size height
        b  = simulate a
        fp = printf "sandpile_%d_%d.pgm" size height
    toPGM b fp
    putStrLn $ "wrote image to " ++ fp
Output:

sandpile with 1000 grains of sand: test 15 1000: [4]
sandpile with 10000 grains of sand: test 40 10000: [5]
sandpile with 100000 grains of sand: test 150 100000: [6]
sandpile with 1000000 grains of sand: test 400 1000000: [7]

J

grid=: 4 : 'x (<<.-:2$y)} (2$y)$0'         NB. y by y grid with x grains in middle
ab=: - [: +/@(-"2 ((,-)=/~i.2)|.!.0]) 3&<  NB. abelian sand pile for grid graph
require 'viewmat'                          NB. viewmat utility
viewmat ab ^: _ (1024 grid 25)             NB. visual

Java

This is based on the JavaScript implementation linked to in the task description.

import java.awt.*;
import java.awt.event.*;
import javax.swing.*;

public class AbelianSandpile {
    public static void main(String[] args) {
        SwingUtilities.invokeLater(new Runnable() {
            public void run() {
                Frame frame = new Frame();
                frame.setVisible(true);
            }
        });
    }

    private static class Frame extends JFrame {
        private Frame() {
            super("Abelian Sandpile Model");
            setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
            Container contentPane = getContentPane();
            JPanel controlPanel = new JPanel(new FlowLayout(FlowLayout.LEFT));
            JButton start = new JButton("Restart Simulation");
            start.addActionListener(e -> restartSimulation());
            JButton stop = new JButton("Stop Simulation");
            stop.addActionListener(e -> stopSimulation());
            controlPanel.add(start);
            controlPanel.add(stop);
            contentPane.add(controlPanel, BorderLayout.NORTH);
            contentPane.add(canvas = new Canvas(), BorderLayout.CENTER);
            timer = new Timer(100, e -> canvas.runAndDraw());
            timer.start();
            pack();
        }

        private void restartSimulation() {
            timer.stop();
            canvas.initGrid();
            timer.start();
        }

        private void stopSimulation() {
            timer.stop();
        }

        private Timer timer;
        private Canvas canvas;
    }

    private static class Canvas extends JComponent {
        private Canvas() {
            setBorder(BorderFactory.createEtchedBorder());
            setPreferredSize(new Dimension(600, 600));
        }

        public void paintComponent(Graphics g) {
            int width = getWidth();
            int height = getHeight();
            g.setColor(Color.WHITE);
            g.fillRect(0, 0, width, height);
            int cellWidth = width/GRID_LENGTH;
            int cellHeight = height/GRID_LENGTH;
            for (int i = 0; i < GRID_LENGTH; ++i) {
                for (int j = 0; j < GRID_LENGTH; ++j) {
                    if (grid[i][j] > 0) {
                        g.setColor(COLORS[grid[i][j]]);
                        g.fillRect(i * cellWidth, j * cellHeight, cellWidth, cellHeight);
                    }
                }
            }
        }

        private void initGrid() {
            for (int i = 0; i < GRID_LENGTH; ++i) {
                for (int j = 0; j < GRID_LENGTH; ++j) {
                    grid[i][j] = 0;
                }
            }
        }

        private void runAndDraw() {
            for (int i = 0; i < 100; ++i)
                addSand(GRID_LENGTH/2, GRID_LENGTH/2);
            repaint();
        }

        private void addSand(int i, int j) {
            int grains = grid[i][j];
            if (grains < 3) {
                grid[i][j]++;
            }
            else {
                grid[i][j] = grains - 3;
                if (i > 0)
                    addSand(i - 1, j);
                if (i < GRID_LENGTH - 1)
                    addSand(i + 1, j);
                if (j > 0)
                    addSand(i, j - 1);
                if (j < GRID_LENGTH - 1)
                    addSand(i, j + 1);
            }
        }

        private int[][] grid = new int[GRID_LENGTH][GRID_LENGTH];
    }

    private static final Color[] COLORS = {
        Color.WHITE,
        new Color(0x00, 0xbf, 0xff),
        new Color(0xff, 0xd7, 0x00),
        new Color(0xb0, 0x30, 0x60)
    };
    private static final int GRID_LENGTH = 300;
}
Output:

Media:Abelian sandpile java.png

jq

Adapted from Wren

Works with jq and gojq, the C and Go implementations of jq

For consistency with Abelian sandpile model/Identity#jq, the function for reducing a sandpile to its equilibrium position is called `avalanche` here.

# Generic functions
def array($n): . as $in | [range(0;$n)|$in];

def when(filter; action): if filter // null then action else . end;

def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;

# module  Sandpile 

# 'a' is a list of integers in row order
def new($a):
  ($a|length) as $length
  | ($length|sqrt|floor) as $rows
  | if ($rows * $rows != $length) then "The matrix of values must be square." | error
    else
    {$a,
     $rows,
     neighbors:
        (reduce range(0; $length) as $i (null;
         .[$i] = []
         | when($i % $rows > 0;       .[$i] += [$i-1] )
         | when(($i + 1) % $rows > 0; .[$i] += [$i+1] )
         | when($i - $rows >= 0;      .[$i] += [$i-$rows] )
         | when($i + $rows < $length; .[$i] += [$i+$rows] ) ) )
    }
    end;

def isStable:
  all(.a[]; . <= 3);

def tos:
  . as $in
  | .rows as $rows
  | reduce range(0; $rows) as $i ("";
      reduce range(0; $rows) as $j (.;
        . + " \($in.a[$rows*$i + $j] | lpad(2))" )
      | . +"\n" );

# just topple once so we can observe intermediate results
def topple:
  last(
    label $out
    | foreach range(0; .a|length) as $i (.;
        if .a[$i] > 3
        then .a[$i] += -4
        | reduce .neighbors[$i][] as $j (.; .a[$j] += 1)
        | ., break $out
       else .
       end ) );

def avalanche:
  until(isStable; topple);

# str1 and str2 should be strings representing a sandpile (i.e. .a)
def printAcross(str1; str2):
  (str1|split("\n")) as $r1
  | (str2|split("\n")) as $r2
  | ($r1|length - 1) as $rows
  | ($rows/2|floor) as $cr
  | reduce range(0; $rows) as $i ("";
        (if $i == $cr then "->" else "  " end) as $symbol
        | . + "\($r1[$i]) \($symbol) \($r2[$i])\n" ) ;

{ a1: (0|array(25))}
| .a2 = .a1
| .a3 = .a1
| .a1[12] = 4
| .a2[12] = 6
| .a3[12] = 16
| .a4 = (0|array(100))
| .a4[55] = 64

| (.a1, .a2, .a3, .a4) as $a
| .s = new($a)
| (.s|tos) as $str1
| .s |= avalanche
| (.s|tos) as $str2
| printAcross($str1; $str2)
Output:

As for Wren


Julia

Modified from code by Hayk Aleksanyan, viewable at github.com/hayk314/Sandpiles, license viewable there.

module AbelSand

# supports output functionality for the results of the sandpile simulations
# outputs the final grid in CSV format, as well as an image file

using CSV, DataFrames, Images

function TrimZeros(A)
    # given an array A trims any zero rows/columns from its borders
    # returns a 4 tuple of integers, i1, i2, j1, j2, where the trimmed array corresponds to A[i1:i2, j1:j2]
    # A can be either numeric or a boolean array

    i1, j1 = 1, 1
    i2, j2 = size(A)

    zz = typeof(A[1, 1])(0)    # comparison of a value takes into account the type as well

    # i1 is the first row which has non zero element
    for i = 1:size(A, 1)
        q = false
        for k = 1:size(A, 2)
            if A[i, k] != zz
                q = true
                i1 = i
                break
            end
        end

        if q == true
            break
        end
    end

    # i2 is the first from below row with non zero element
    for i in size(A, 1):-1:1
        q = false
        for k = 1:size(A, 2)
            if A[i, k] != zz
                q = true
                i2 = i
                break
            end
        end

        if q == true
            break
        end
    end

    # j1 is the first column with non zero element

    for j = 1:size(A, 2)
        q = false
        for k = 1:size(A, 1)
            if A[k, j] != zz
                j1 = j
                q = true
                break
            end
        end

        if q == true
            break
        end
    end

    # j2 is the last column with non zero element

    for j in size(A, 2):-1:1
        q=false
        for k=1:size(A,1)
            if A[k, j] != zz
                j2 = j
                q=true
                break
            end
        end

        if q==true
            break
        end
    end

    return i1, i2, j1, j2
end

function addLayerofZeros(A, extraLayer)
    # adds layer of zeros from all corners to the given array A

    if extraLayer <= 0
        return A
    end

    N, M = size(A)


    Z = zeros( typeof(A[1,1]), N + 2*extraLayer, M + 2*extraLayer)
    Z[(extraLayer+1):(N + extraLayer ), (extraLayer+1):(M+extraLayer)] = A

    return Z

end

function printIntoFile(A, extraLayer, strFileName, TrimSmallValues = false)
    # exports a 2d matrix A into a csv file
    # @extraLayer is an integers adding layer of 0-s sorrounding the output matrix

    # trimming off very small values; tiny values affect the performance of CSV export
    if TrimSmallValues == true
        A = map(x -> if (abs(x - floor(x)) < 0.01) floor(x) else x end, A) 
    end

    i1, i2, j1, j2  = TrimZeros( A )
    A = A[i1:i2, j1:j2]

    A = addLayerofZeros(A, extraLayer)

    CSV.write(string(strFileName,".csv"), DataFrame(A), writeheader = false)

    return A

end

function Array_magnifier(A, cell_mag, border_mag)
    # A is the main array; @cell_mag is the magnifying size of the cell,
    # @border_mag is the magnifying size of the border between lattice cells

    # creates a new array where each cell of the original array A appears magnified by size = cell_mag


    total_factor = cell_mag + border_mag

    A1 = zeros(typeof(A[1, 1]), total_factor*size(A, 1), total_factor*size(A, 2))

    for i = 1:size(A,1), j = 1:size(A,2), u = ((i-1)*total_factor+1):(i*total_factor),
                                          v = ((j-1)*total_factor+1):(j*total_factor)
        if(( u - (i - 1) * total_factor <= cell_mag) && (v - (j - 1) * total_factor <= cell_mag))
            A1[u, v] = A[i, j] 
        end
    end

    return A1

end

function saveAsGrayImage(A, fileName, cell_mag, border_mag, TrimSmallValues = false)
    # given a 2d matrix A, we save it as a gray image after magnifying by the given factors
    A1 = Array_magnifier(A, cell_mag, border_mag)
    A1 = A1/maximum(maximum(A1))

    # trimming very small values from A1 to improve performance
    if TrimSmallValues == true
        A1 = map(x -> if ( x < 0.01) 0.0 else round(x, digits = 2) end, A1) 
    end

    save(string(fileName, ".png") , colorview(Gray, A1)) 
end

function saveAsRGBImage(A, fileName, color_codes, cell_mag, border_mag)
    # color_codes is a dictionary, where key is a value in A and value is an RGB triplet
    # given a 2d array A, and color codes (mapping from values in A to RGB triples), save A
    # into fileName as png image after applying the magnifying factors

    A1 = Array_magnifier(A, cell_mag, border_mag)
    color_mat = zeros(UInt8, (3, size(A1, 1), size(A1, 2)))

    for i = 1:size(A1,1)
        for j = 1:size(A1,2)
            color_mat[:, i, j]  = get(color_codes, A1[i, j] , [0, 0, 0])
        end
    end

    save(string(fileName, ".png") , colorview(RGB, color_mat/255)) 
end

const N_size = 700       # the radius of the lattice Z^2, the actual size becomes (2*N+1)x(2*N+1)
const dx = [1, 0, -1, 0] # for a given (x,y) in Z^2, (x + dx, y + dy) for all (dx,dy) covers the neighborhood of (x,y)
const dy = [0, 1, 0, -1]

struct L_coord
    # represents a lattice coordinate
    x::Int
    y::Int
end

function FindCoordinate(Z::Array{L_coord,1}, a::Int, b::Int)
    # in the given array Z of coordinates finds the (first) index of the tuple (a,b)
    # if no match, returns -1

    for i=1:length(Z)
        if (Z[i].x == a) && (Z[i].y == b)
            return i
        end
    end

    return -1
end

function move(N)
    # the main function moving the pile sand grains of size N at the origin of Z^2 until the sandpile becomes stable

    Z_lat = zeros(UInt8, 2 * N_size + 1, 2 * N_size + 1)     # models the integer lattice Z^2, we will have at most 4 sands on each vertex
    V_sites = falses(2 * N_size + 1, 2 * N_size + 1)         # all sites which are visited by the sandpile process, are being marked here
    Odometer = zeros(UInt64, 2 * N_size + 1, 2 * N_size + 1) # stores the values of the odometer function


    walking = L_coord[]    # the coordinates of sites which need to move

    V_sites[N_size + 1, N_size + 1] = true

    # i1, ... j2  -> show the boundaries of the box which is visited by the sandpile process
    i1, i2, j1, j2 = N_size + 1, N_size + 1, N_size + 1, N_size + 1 
    n = N

    t1 = time_ns()
    
    while n > 0
        n -= 1

        Z_lat[N_size + 1, N_size + 1] += 1
        if (Z_lat[N_size + 1, N_size + 1] >= 4)
            push!(walking, L_coord(N_size + 1, N_size + 1))
        end

        while(length(walking) > 0)
            w = pop!(walking)
            x = w.x
            y = w.y

            Z_lat[x, y] -= 4
            Odometer[x, y] += 4

            for k = 1:4
                Z_lat[x + dx[k], y + dy[k]] += 1
                V_sites[x + dx[k], y + dy[k]] = true
                if Z_lat[x + dx[k], y + dy[k]] >= 4
                    if FindCoordinate(walking, x + dx[k] , y + dy[k]) == -1
                        push!(walking, L_coord( x + dx[k], y + dy[k]))
                    end
                end
            end

            i1 = min(i1, x - 1)
            i2 = max(i2, x + 1)
            j1 = min(j1, y - 1)
            j2 = max(j2, y + 1)
        end


    end #end of the main while
    t2 = time_ns()

    println("The final boundaries are:: ", (i2 - i1 + 1),"x",(j2 - j1 + 1), "\n")
    print("time elapsed: " , (t2 - t1) / 1.0e9, "\n")

    Z_lat = printIntoFile(Z_lat, 0, string("Abel_Z_", N))
    Odometer = printIntoFile(Odometer, 1, string("Abel_OD_", N))

    saveAsGrayImage(Z_lat, string("Abel_Z_", N), 20, 0)
    color_code = Dict(1=>[255, 128, 255], 2=>[255, 0, 0],3 => [0, 128, 255])
    saveAsRGBImage(Z_lat, string("Abel_Z_color_", N), color_code, 20, 0)

    # for the total elapsed time, it's better to use the @time macros on the main call

    return Z_lat, Odometer # these are trimmed in output module

end # end of function move


end # module


using .AbelSand

Z_lat, Odometer = AbelSand.move(100000)
Output:

Link to PNG output file for N=100000 ie. AbelSand.move(100000)
Link to PNG output file (run time >90 min) for N=1000000 (move(1000000))

Lua

local sandpile = {
  init = function(self, dim, val)
    self.cell, self.dim = {}, dim
    for r = 1, dim do
      self.cell[r] = {}
      for c = 1, dim do
        self.cell[r][c] = 0
      end
    end
    self.cell[math.floor(dim/2)+1][math.floor(dim/2)+1] = val
  end,
  iter = function(self)
    local dim, cel, more = self.dim, self.cell
    repeat
      more = false
      for r = 1, dim do
        for c = 1, dim do
          if cel[r][c] >= 4 then
            cel[r][c] = cel[r][c] - 4
            if c > 1 then cel[r][c-1], more = cel[r][c-1]+1, more or cel[r][c-1]>=3 end
            if c < dim then cel[r][c+1], more = cel[r][c+1]+1, more or cel[r][c+1]>=3 end
            if r > 1 then cel[r-1][c], more = cel[r-1][c]+1, more or cel[r-1][c]>=3 end
            if r < dim then cel[r+1][c], more = cel[r+1][c]+1, more or cel[r+1][c]>=3 end
          end
          more = more or cel[r][c] >= 4
        end
      end
    until not more
  end,
  draw = function(self)
    for r = 1, self.dim do
      print(table.concat(self.cell[r]," "))
    end
  end,
}
sandpile:init(15, 256)
sandpile:iter()
sandpile:draw()
Output:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 2 2 2 1 0 0 0 0 0
0 0 0 0 2 1 0 2 0 1 2 0 0 0 0
0 0 0 2 3 3 3 2 3 3 3 2 0 0 0
0 0 2 3 2 2 2 3 2 2 2 3 2 0 0
0 1 1 3 2 2 3 0 3 2 2 3 1 1 0
0 2 0 3 2 3 2 3 2 3 2 3 0 2 0
0 2 2 2 3 0 3 0 3 0 3 2 2 2 0
0 2 0 3 2 3 2 3 2 3 2 3 0 2 0
0 1 1 3 2 2 3 0 3 2 2 3 1 1 0
0 0 2 3 2 2 2 3 2 2 2 3 2 0 0
0 0 0 2 3 3 3 2 3 3 3 2 0 0 0
0 0 0 0 2 1 0 2 0 1 2 0 0 0 0
0 0 0 0 0 1 2 2 2 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Mathematica/Wolfram Language

ClearAll[sp]
sp[s_List] + sp[n_Integer] ^:= sp[s] + sp[ConstantArray[n, Dimensions[s]]]
sp[s_List] + sp[t_List] ^:= Module[{dim, r, tmp, neighbours}, dim = Dimensions[s];
  r = s + t;
  While[Max[r] > 3, r = ArrayPad[r, 1, 0];
   tmp = Quotient[r, 4];
   r -= 4 tmp;
   r += RotateLeft[tmp, {0, 1}] + RotateLeft[tmp, {1, 0}] + 
     RotateLeft[tmp, {0, -1}] + RotateLeft[tmp, {-1, 0}];
   r = ArrayPad[r, -1];];
  sp[r]
  ]
u = sp[CenterArray[250, {15, 15}]];
u += sp[0];
StringRiffle[StringJoin /@ Map[ToString, u[[1]], {2}], "\n"]
Output:
000000000000000
000001222100000
000021020120000
000233323332000
002322222223200
011322232223110
020322030223020
022223323322220
020322030223020
011322232223110
002322222223200
000233323332000
000021020120000
000001222100000
000000000000000

MiniScript

For use with the Mini Micro.

colors = [color.black, color.yellow, color.orange, 
          color.brown, color.red, color.fuchsia, 
          color.purple, color.blue, color.navy]

rows = 48; rowRange = range(0, rows-1)
cols = 72; colRange = range(0, cols-1)
particlesOfSand = rows * cols
divBase = particlesOfSand / (colors.len - 4)
deltas = [[0,-1],[-1, 0], [1, 0],[0, 1]]

displayGrid = function(grid, td)
	for y in globals.rowRange
		for x in globals.colRange
			colorIx = grid[y][x]
			// determine the rest of the colors if > 3 by division 
			if colorIx > 3 then colorIx = (colorIx - 3) / divBase + 4
			td.setCell x,y, colorIx
		end for
	end for
end function

clear

// Prepare a tile display
// Generate image used for the tiles from the defined above.
// The colors are to indicate height of a sand pile.
img = Image.create(colors.len, 1); 
for i in range(0, colors.len - 1)
	img.setPixel(i, 0, colors[i])
end for

grid = []
for y in rowRange
	row = []
	for x in colRange
		row.push(0)
	end for
	grid.push(row)
end for

grid[rows/2][cols/2] = particlesOfSand

display(4).mode = displayMode.tile
td = display(4)
td.cellSize = 640/48   // size of cells on screen
td.extent = [cols, rows]
td.overlap = 0   // adds a small gap between cells
td.tileSet = img; td.tileSetTileSize = 1
td.clear 0

toTopple = []
for y in rowRange
	for x in colRange
		if grid[y][x] > 3 and toTopple.indexOf([x,y]) == null then toTopple.push([x,y])
	end for
end for
tt = time
while toTopple.len > 0
	nextGen = []
	for cell in toTopple
		x = cell[0]; y = cell[1]
		grid[y][x] -= 4
		for delta in deltas
			x1 = (x + delta[0]) % cols; y1 = (y + delta[1]) % rows
			grid[y1][x1] += 1
		end for
	end for
	for y in rowRange
		for x in colRange
			if grid[y][x] > 3 and nextGen.indexOf([x,y]) == null then nextGen.push([x,y])
		end for
	end for
	toTopple = nextGen
	displayGrid(grid, td)
end while
key.get()
File:Miniscript abelian sandpille.png
Image for 3456 particles grid 48*72

Nim

Library: nimPNG

Our program uses Rust algorithm (and also its colors 🙂) and the formula to compute grid size from number of particles comes from Pascal algorithm. Number of particles is an input from user. The program displays the values on the terminal if there are not too many and produce a PNG image. Code to produce a PPM image is also provided but not used.

# Abelian sandpile.

from math import sqrt
from nimPNG import savePNG24
from sequtils import repeat
from strformat import fmt
from strutils import strip, addSep, parseInt

# The grid represented as a sequence of sequences of int32.
type Grid = seq[seq[int32]]

# Colors to use for PPM and PNG files.
const Colors = [[byte 100,  40,  15],
                [byte 117,  87,  30],
                [byte 181, 134,  47],
                [byte 245, 182,  66]]

#---------------------------------------------------------------------------------------------------

func sideLength(initVal: int32): int32 =
  # Return the grid side length needed for "initVal" particles.
  # We make sure that the returned value is odd.
  result = sqrt(initVal.toFloat / 1.75).int32 + 3
  result += result and 1 xor 1

#---------------------------------------------------------------------------------------------------

func doOneStep(grid: var Grid; boundary: var array[4, int]): bool =
  ## Compute one step.

  result = false

  for y in boundary[0]..boundary[2]:
    for x in boundary[1]..boundary[3]:
      if grid[y][x] >= 4:

        let rem = grid[y][x] div 4
        grid[y][x] = grid[y][x] mod 4

        if y - 1 >= 0:
          inc grid[y - 1][x], rem
          if y == boundary[0]:
            dec boundary[0]

        if x - 1 >= 0:
          inc grid[y][x - 1], rem
          if x == boundary[1]:
            dec boundary[1]

        if y + 1 < grid.len:
          inc grid[y + 1][x], rem
          if y == boundary[2]:
            inc boundary[2]

        if x + 1 < grid.len:
          inc grid[y][x + 1], rem
          if x == boundary[3]:
            inc boundary[3]

        result = true

#---------------------------------------------------------------------------------------------------

proc display(grid: Grid; initVal: int) =
  ## Display the grid as an array of values.

  echo fmt"Starting with {initVal} particles."
  echo ""

  var line = newStringOfCap(2 * grid.len - 1)
  for row in grid:
    for value in row:
      line.addSep(" ", 0)
      line.add($value)
    echo line
    line.setLen(0)
  echo ""

#---------------------------------------------------------------------------------------------------

proc writePpmFile(grid: Grid; name: string) =
  ## Write a grid representation in a PPM file.

  var file = open(name, fmWrite)
  file.write(fmt"P6 {grid.len} {grid.len} 255 ")

  for row in grid:
    for value in row:
      discard file.writeBytes(Colors[value], 0, 3)

  file.close()
  echo fmt"PPM image written in ""{name}""."

#---------------------------------------------------------------------------------------------------

proc writePngFile(grid: Grid; name: string) =
  ## Write a grid representation in a PNG file.

  var pixels = newSeq[byte](3 * grid.len * grid.len)

  # Build pixel list.
  var idx = 0
  for row in grid:
    for value in row:
      pixels[idx..idx+2] = Colors[value]
      inc idx, 3

  discard savePNG24(name, pixels, grid.len, grid.len)
  echo fmt"PNG image written in ""{name}""."

#---------------------------------------------------------------------------------------------------

proc askInitVal(): int32 =
  # Ask user for the number of particles.

  while true:
    stdout.write("Number of particles? ")
    try:
      let input = stdin.readLine().strip().parseInt()
      if input in 4..int32.high:
        return input.int32
      echo fmt"Value not in expected range: 4..{int32.high}"
    except ValueError:
      echo "Invalid input"
    except EOFError:
      quit(QuitSuccess)

#---------------------------------------------------------------------------------------------------

# Initialize the grid.
let initVal = askInitVal()
let sideLen = sideLength(initVal)
var grid = repeat(newSeq[int32](sideLen), sideLen)
let origin = grid.len div 2
var boundaries: array[4, int] = [origin, origin, origin, origin]
grid[origin][origin] = initVal

# Run the simulation.
while doOneStep(grid, boundaries):
  discard

# Display grid.
if grid.len <= 20:
  grid.display(initVal)
#grid.writePpmFile(fmt"grid_{initVal}.ppm")
grid.writePngFile(fmt"grid_{initVal}.png")
Output:
Number of particles? 100
Starting with 100 particles.

0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 2 1 0 0 0 0
0 0 0 3 2 0 2 3 0 0 0
0 0 3 0 3 2 3 0 3 0 0
0 1 2 3 0 3 0 3 2 1 0
0 2 0 2 3 0 3 2 0 2 0
0 1 2 3 0 3 0 3 2 1 0
0 0 3 0 3 2 3 0 3 0 0
0 0 0 3 2 0 2 3 0 0 0
0 0 0 0 1 2 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0

PNG image written in "grid_100.png".

Pascal

Works with: Free Pascal
The main optimization was to spread the sand immediatly.
mul := val DIV 4;//not only := val -4 
so that only (sand mod 4) stays in place.runtime for abelian(1e6) down to 1min 20 secs from 9 min

Memorizing the used colums of the rows has little effect when choosing the right size of the grid.Only 11 secs for abelian(1e6) -> 1min 9sec
Python shows 64 too.

program Abelian2;
{$IFDEF FPC}
   {$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$CODEALIGN proc=16}{$ALIGN 16}
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}
uses
  SysUtils;

type
  Tlimit = record
             lmtLow,LmtHigh : LongWord;
           end;
  TRowlimits = array of Tlimit;
  tOneRow  = pLongWord;
  tGrid = array of LongWord;

var
  Grid: tGrid;
  Rowlimits:TRowlimits;
  s : AnsiString;
  maxval,maxCoor : NativeUint;

function CalcMaxCoor(maxVal : NativeUint):NativeUint;
//  maxVal = 10000;maxCoor = 77-2;// maxCoor*maxCoor    *1,778;     0.009sec
//  maxVal = 100000;maxCoor = 236-2;// maxCoor*maxCoor  *1.826;     0.825sec
//  maxVal = 1000000;maxCoor = 732-2;// maxCoor*maxCoor *1.877;    74    sec
Begin
  result := trunc(sqrt(maxval/1.75))+3;
end;

procedure clear;
begin
  setlength(Grid,0);
  setlength(Rowlimits,0);
  s := '';
end;

procedure InitGrid(var G:tGrid;InitVal:NativeUint);
var
  row,middle: nativeINt;
begin
//  setlength(Rowlimits,0);   setlength(G,0);
  MaxCoor :=  CalcMaxCoor(InitVal);
  setlength(G,sqr(maxCoor));
  setlength(Rowlimits,maxCoor);
  fillchar(G[0],length(G)*SizeOf(G[0]),#0);

  middle := (maxCoor) div 2;
  Grid[middle*maxcoor+middle] := InitVal;
  For row := 1 to maxCoor do
    with Rowlimits[row] do
    Begin
      lmtLow := middle;
      lmtHigh := middle;
    end;

  with Rowlimits[middle] do
  Begin
    lmtLow := middle;
    lmtHigh := middle;
  end;
end;
procedure OutGridPPM(const G:tGrid;maxValue : NativeUint);
const
  color : array[0..3] of array[0..2] of Byte =
             //R,G,B)
            ((0,0,0),
             (255,0,0),
             (0,255,0),
             (0,0,255));
var
  f :text;
  pActRow: tOneRow;
  col,row,sIdx,value : NativeInt;
Begin
  Assignfile(f,'ppm/Grid_'+IntToStr(maxValue)+'.ppm');
  rewrite(f);
  write(f,Format('P6 %d %d %d ',[maxCoor-1,maxCoor-1,255]));
  setlength(s,(maxCoor-1)*3);
  pActRow :=@G[0];
  For row := maxCoor-2 downto 0 do
  Begin
    inc(pActRow,maxCoor);
    sIdx := 1;
    For col := 1 to maxCoor-1 do
    Begin
      value := pActRow[col];
      s[sIdx]   := CHR(color[value,0]);
      s[sIdx+1] := CHR(color[value,1]);
      s[sIdx+2] := CHR(color[value,2]);
      inc(sIdx,3);
    end;
    write(f,s);
  end;
  CloseFile(f);
end;

procedure OutGrid(const G:tGrid);
//output of grid and test, if no sand is lost
var
  pActRow: tOneRow;
  col,row,sum,value : NativeUint;
Begin
  setlength(s,maxcoor-1);
  pActRow := @G[0];
  sum := 0;
  For row := maxCoor-1 downto 1 do
  Begin
    inc(pActRow,maxcoor);
    For col := 1 to maxCoor-1 do
    Begin
      value := pActRow[col];
//      IF value>=4 then writeln(row:5,col:5,value:13);
      s[col] := chr(value+48);
      inc(sum,value);
    end;
    if maxCoor <80 then
      writeln(s);
  end;
  writeln('columns ',maxcoor-1,' checksum ',maxVal,' ?=? ',sum);
{
  For row := 1 to maxCoor do
    with Rowlimits[row] do
      writeln(lmtLow:10,lmtHigh:10);
      * }
end;

procedure Evolution(var G:tGrid);
var
  pActRow,pRowBefore,pRowAfter : tOneRow;
  col,row,mul,val,done : NativeUint;
begin
  repeat
    pRowBefore := @G[0];
    pActRow    := @G[maxcoor];
    pRowAfter  := @G[2*maxcoor];
    done := 0;
    For row := maxCoor-1 downto 1 do
    Begin
      with RowLimits[row] do
      Begin
      while (LmtLow >1) AND (pActRow[lmtLow]<> 0) do
        dec(lmtLow);
      while (lmtHigh < maxCoor) AND (pActRow[lmtHigh]<> 0) do
        inc(lmtHigh);
      For col := lmtLow to lmtHigh do
      Begin
        val := pActRow[col];
        IF val >=4 then
        Begin
          mul := val DIV 4;
          done := val;
          inc(pRowBefore[col],mul);
          inc(pActRow[col-1],mul);
          pActRow[col] := val-4*Mul;
          inc(pActRow[col+1],mul);
          inc(pRowAfter[col],mul);
        end;
      end;
      pRowBefore:= pActRow;
      pActRow := pRowAfter;
      inc(pRowAfter,maxcoor);
    end;
    end;
  until done=0;
end;

procedure OneTurn(count:NativeUint);
begin
  Writeln(' Test abelian sandpile( ',count,' )');
  MaxVal := count;
  InitGrid(Grid,count);
  Evolution(Grid);
  OutGrid(Grid);
  OutGridPPM(Grid,count);
  clear;
end;

BEGIN
  OneTurn(4);
  OneTurn(16);
  OneTurn(64);
  OneTurn(1000);
  OneTurn(10000);
  OneTurn(100000);
END.
Output:
 Test abelian sandpile( 4 )
010
101
010
columns 3 checksum 4 ?=? 4
 Test abelian sandpile( 16 )
00100
02120
11011
02120
00100
columns 5 checksum 16 ?=? 16
 Test abelian sandpile( 64 )
00121000
02222200
12222210
22202220
12222210
02222200
00121000
00000000
columns 8 checksum 64 ?=? 64
 Test abelian sandpile( 1000 )
0000000001111111000000000
0000000130233320310000000
0000013223313133223100000
0000213222130312223120000
0002220123332333210222000
0011223233123213323221100
0033032313221223132303300
0122123203311133023212210
0322231023333333201322230
1032333332231322333332301
1231312332232322332131321
1313322133322233312233131
1330231131220221311320331
1313322133322233312233131
1231312332232322332131321
1032333332231322333332301
0322231023333333201322230
0122123203311133023212210
0033032313221223132303300
0011223233123213323221100
0002220123332333210222000
0000213222130312223120000
0000013223313133223100000
0000000130233320310000000
0000000001111111000000000
columns 25 checksum 1000 ?=? 1000
 Test abelian sandpile( 10000 )
--shortened
columns 77 checksum 10000 ?=? 10000
 Test abelian sandpile( 100000 )
columns 241 checksum 100000 ?=? 100000

real    0m0,815s

OCaml

Works with: OCaml version 4.11

In Sandpile module (sandpile.ml)

module Make =
  functor (M : sig val m : int val n : int end)
  -> struct

    let grid = Array.init M.m (fun _ -> Array.make M.n 0)

    let print () =
      for i = 0 to M.m - 1
      do for j = 0 to M.n - 1
         do Printf.printf "%d " grid.(i).(j)
         done
       ; print_newline ()
       done

    let unstable = Hashtbl.create 10

    let add_grain x y
      = grid.(x).(y) <- grid.(x).(y) + 1
      ; if grid.(x).(y) >= 4 then
          Hashtbl.replace unstable (x,y) () (* Use Hashtbl.replace for uniqueness *)
      
    let topple x y
      = grid.(x).(y) <- grid.(x).(y) - 4
      ; if grid.(x).(y) < 4
        then Hashtbl.remove unstable (x,y)
      ; match (x,y) with
        (* corners *)
        | (0,0) -> add_grain 1 0
                 ; add_grain 0 1
        | (0,n) when n = M.n - 1
          -> add_grain 1 n
           ; add_grain 0 (n-1)
        | (m,0) when m = M.m - 1
          -> add_grain m 1
           ; add_grain (m-1) 0
        | (m,n) when m = M.m - 1 && n = M.n - 1
          -> add_grain ( m ) (n-1)
           ; add_grain (m-1) ( n )
        (* sides *)
        | (0,y) -> add_grain 1 y
                 ; add_grain 0 (y+1)
                 ; add_grain 0 (y-1)
        | (m,y) when m = M.m - 1
          -> add_grain ( m ) (y-1)
           ; add_grain ( m ) (y+1)
           ; add_grain (m-1) ( y )
        | (x,0) -> add_grain (x+1) 0
                 ; add_grain (x-1) 0
                 ; add_grain ( x ) 1
        | (x,n) when n = M.n - 1
          -> add_grain (x-1) ( n )
           ; add_grain (x+1) ( n )
           ; add_grain ( x ) (n-1)
        (* else *)
        | (x,y) -> add_grain ( x ) (y+1)
                 ; add_grain ( x ) (y-1)
                 ; add_grain (x+1) ( y )
                 ; add_grain (x-1) ( y )
                 
    let add_sand n x y
      = for i = 1 to n
        do add_grain x y
        done

    let avalanche ()
      = while Hashtbl.length unstable > 0 
        do
         let unstable' = Hashtbl.fold  (fun (x,y) () r -> (x,y) :: r) unstable []
         in List.iter (fun (x,y) -> topple x y ) unstable'
        done 
  end

(* testing *)
   
let ()
  = let module S = Make (struct let m = 11 let n = 11 end)
    in S.add_sand 500 5 5
     ; S.avalanche ()
     ; S.print ()

Perl

use strict;
use warnings;
use feature 'bitwise';

my ($high, $wide) = split ' ', qx(stty size);
my $mask = "\0" x $wide . ("\0" . "\177" x ($wide - 2) . "\0") x ($high - 5) .
  "\0" x $wide;
my $pile = $mask =~ s/\177/ rand() < 0.02 ? chr 64 + rand 20 : "\0" /ger;

for (1 .. 1e6)
  {
  print "\e[H", $pile =~ tr/\0-\177/ 1-~/r, "\n$_";
  my $add = $pile =~ tr/\1-\177/\0\0\0\200/r; # set high bit for >=4
  $add =~ /\200/ or last;
  $pile =~ tr/\4-\177/\0-\173/; # subtract 4 if >=4
  for ("\0$add", "\0" x $wide . $add, substr($add, 1), substr $add, $wide)
    {
    $pile |.= $_;
    $pile =~ tr/\200-\377/\1-\176/; # add one to each neighbor of >=4
    $pile &.= $mask;
    }
  select undef, undef, undef, 0.1; # comment out for full speed
  }

Phix

Library: Phix/pGUI

Generates moving images similar to the julia output. The distributed version also has variable speed, additional display modes, and a random dropping toggle.

-- demo\rosetta\Abelian_sandpile_model.exw
include pGUI.e
 
Ihandle dlg, canvas
cdCanvas cddbuffer
 
sequence board = {{0,0,0},
                  {0,0,0},
                  {0,0,0}}
 
procedure drop(integer y, x)
    sequence moves = {}
    while true do
        board[y,x] += 1
        if board[y,x]>=4 then
            board[y,x] -= 4
            moves &= {{y,x-1},{y,x+1},{y-1,x},{y+1,x}}
        end if
        -- extend board if rqd (maintain a border of zeroes)
        if x=1 then                             -- extend left
            for i=1 to length(board) do
                board[i] = prepend(board[i],0)
            end for
            for i=1 to length(moves) do
                moves[i][2] += 1
            end for
        elsif x=length(board[1]) then           -- extend right
            for i=1 to length(board) do
                board[i] = append(board[i],0)
            end for
        end if
        -- (copy the all-0 lines from the other end...)
        if y=1 then                             -- extend up
            board = prepend(board,board[$])
            for i=1 to length(moves) do
                moves[i][1] += 1
            end for
        elsif y=length(board) then              -- extend down
            board = append(board,board[1])
        end if
        if length(moves)=0 then exit end if
        {y,x} = moves[$]
        moves = moves[1..$-1]
    end while   
    IupUpdate(canvas)
end procedure
 
function timer_cb(Ihandle /*ih*/)
    integer y = floor(length(board)/2)+1,
            x = floor(length(board[1])/2)+1
    drop(y,x)
    return IUP_DEFAULT
end function
 
function redraw_cb(Ihandle ih, integer /*posx*/, integer /*posy*/)
    IupGLMakeCurrent(ih)
    cdCanvasActivate(cddbuffer)
    cdCanvasClear(cddbuffer)
    for y=1 to length(board) do
        for x=1 to length(board[1]) do 
            integer c = board[y][x]
            if c!=0 then
                integer colour = {CD_VIOLET,CD_RED,CD_BLUE}[c]
                cdCanvasPixel(cddbuffer, x, y, colour)
            end if
        end for
    end for
    cdCanvasFlush(cddbuffer)
    return IUP_DEFAULT
end function
 
function map_cb(Ihandle ih)
    IupGLMakeCurrent(ih)
    atom res = IupGetDouble(NULL, "SCREENDPI")/25.4
    cddbuffer = cdCreateCanvas(CD_GL, "300x100 %g", {res})
    cdCanvasSetBackground(cddbuffer, CD_PARCHMENT)
    return IUP_DEFAULT
end function
 
procedure main()
    IupOpen()
    canvas = IupGLCanvas("RASTERSIZE=300x100")
    IupSetCallbacks({canvas}, {"ACTION", Icallback("redraw_cb"),
                               "MAP_CB", Icallback("map_cb")})
    dlg = IupDialog(canvas,"TITLE=\"Abelian sandpile model\"")
    IupShow(dlg)
    Ihandle timer = IupTimer(Icallback("timer_cb"), 10)
    IupMainLoop()
    IupClose()
end procedure
 
main()

PicoLisp

(load "@lib/simul.l")
(symbols 'simul 'pico)
(de sandpile (A B)
   (let
      (Grid (grid A A)
         Size (/ (inc A) 2)
         Center (get Grid Size Size)
         Done T )
      (for G Grid
         (for This G
            (=: V 0) ) )
      (with Center
         (=: V B)
         (while Done
            (off Done)
            (for G Grid
               (for This G
                  (when (>= (: V) 4)
                     (=: V (- (: V) 4))
                     (on Done)
                     (mapc
                        '((Dir)
                           (with (Dir This) (=: V (inc (: V)))) )
                        '(north south west east) ) ) ) ) ) )
      (disp Grid 0
         '((This) (if (: V) (pack " " @ " ") "   ")) ) ) )
(sandpile 10 64)
Output:
   +---+---+---+---+---+---+---+---+---+---+
 10 | 0   0   0   0   0   0   0   0   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  9 | 0   0   0   0   0   0   0   0   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  8 | 0   0   0   1   2   1   0   0   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  7 | 0   0   2   2   2   2   2   0   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  6 | 0   1   2   2   2   2   2   1   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  5 | 0   2   2   2   0   2   2   2   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  4 | 0   1   2   2   2   2   2   1   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  3 | 0   0   2   2   2   2   2   0   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  2 | 0   0   0   1   2   1   0   0   0   0 |
    +   +   +   +   +   +   +   +   +   +   +
  1 | 0   0   0   0   0   0   0   0   0   0 |
    +---+---+---+---+---+---+---+---+---+---+
      a   b   c   d   e   f   g   h   i   j

Python

Python: Original, with output

import numpy as np
import matplotlib.pyplot as plt


def iterate(grid):
    changed = False
    for ii, arr in enumerate(grid):
        for jj, val in enumerate(arr):
            if val > 3:
                grid[ii, jj] -= 4
                if ii > 0:
                    grid[ii - 1, jj] += 1
                if ii < len(grid)-1:
                    grid[ii + 1, jj] += 1
                if jj > 0:
                    grid[ii, jj - 1] += 1
                if jj < len(grid)-1:
                    grid[ii, jj + 1] += 1
                changed = True
    return grid, changed


def simulate(grid):
    while True:
        grid, changed = iterate(grid)
        if not changed:
            return grid


if __name__ == '__main__':
    start_grid = np.zeros((10, 10))
    start_grid[4:5, 4:5] = 64
    final_grid = simulate(start_grid.copy())
    plt.figure()
    plt.gray()
    plt.imshow(start_grid)
    plt.figure()
    plt.gray()
    plt.imshow(final_grid)
Output:

Before:

[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0.64. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

After:

[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 2. 1. 0. 0. 0. 0.]
 [0. 0. 2. 2. 2. 2. 2. 0. 0. 0.]
 [0. 1. 2. 2. 2. 2. 2. 1. 0. 0.]
 [0. 2. 2. 2. 0. 2. 2. 2. 0. 0.]
 [0. 1. 2. 2. 2. 2. 2. 1. 0. 0.]
 [0. 0. 2. 2. 2. 2. 2. 0. 0. 0.]
 [0. 0. 0. 1. 2. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

An interactive variant to the above solution:

from os import system, name
from time import sleep

def clear():
	if name == 'nt':
		_ = system('cls')
	else: _ = system('clear')

def exit():
	import sys
	clear()
	sys.exit()

def make_area(x, y):
	area = [[0]*x for _ in range(y)]
	return area

def make_sandpile(area, loc, height):
	loc=list(n-1 for n in loc)
	x, y = loc

	try: area[y][x]+=height
	except IndexError: pass
	
def run(area, by_frame=False):
	def run_frame():
		for y_index, group in enumerate(area):
			y = y_index+1

			for x_index, height in enumerate(group):
				x = x_index+1

				if height < 4: continue

				else:
					make_sandpile(area, (x+1, y), 1)
					make_sandpile(area, (x, y+1), 1)

					if x_index-1 >= 0:
						make_sandpile(area, (x-1, y), 1)
					if y_index-1 >= 0:
						make_sandpile(area, (x, y-1), 1)

					make_sandpile(area, (x, y), -4)

	while any([any([pile>=4 for pile in group]) for group in area]):
		if by_frame:
			clear()
		run_frame()
		if by_frame:
			show_area(area); sleep(.05)

def show_area(area):
	display = [' '.join([str(item) if item else ' ' for item in group]) 
			   for group in area]
	[print(i) for i in display]

clear()
if __name__ == '__main__':
	area = make_area(10, 10)
	print('\nBefore:')
	show_area(area)
	make_sandpile(area, (5, 5), 64)
	run(area)
	print('\nAfter:')
	show_area(area)
Output:
Before:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0

After:
0 0 0 0 0 0 0 0 0 0
0 0 0 1 2 1 0 0 0 0
0 0 2 2 2 2 2 0 0 0
0 1 2 2 2 2 2 1 0 0
0 2 2 2 0 2 2 2 0 0
0 1 2 2 2 2 2 1 0 0
0 0 2 2 2 2 2 0 0 0
0 0 0 1 2 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0

Python: using tkinter

''' Python 3.6.5 code using Tkinter graphical user
    interface (Canvas widget) to display final results.'''
from tkinter import *

# given a grid, display it on a tkinter Canvas:
class Sandpile:
    def __init__(self, wn, grid):
        self.window = wn
        self.grid = grid
        self.canvas = Canvas(wn, bg='lemon chiffon')
        self.canvas.pack(fill=BOTH, expand=1)       

        colors = {0:'dodger blue',
                  1:'red',
                  2:'green',
                  3:'lemon chiffon'}
        
        x = 10
        y = 10
        d = 5
        
        for row in self.grid:
            for value in row:
                clr = colors[value]
                self.canvas.create_rectangle(
                    x, y, x+d, y+d,
                    outline=clr,
                    fill = clr)
                x += 5
            x = 10
            y += 5
            
class Grid:
    def __init__(self, size, center):
        self.size = size        # rows/cols in (best if odd)
        self.center = center    # start value at center of grid
        self.grid = [[0]*self.size for i in range(self.size)]
        self.grid[self.size // 2][self.size // 2] = self.center

    # print the grid:
    def show(self, msg):
        print('  ' + msg + ':')
        for row in self.grid:
            print(' '.join(str(x) for x in row))       
        print()
        return

    # dissipate piles of sand as required:
    def abelian(self):
        while True:
            found = False
            for r in range(self.size):
                for c in range(self.size):
                    if self.grid[r][c] > 3:
                        self.distribute(self.grid[r][c], r, c)
                        found = True
            if not found:
                return

    # distribute sand from a single pile to its neighbors:
    def distribute(self, nbr, row, col):
        qty, remain = divmod(nbr, 4)
        self.grid[row][col] = remain
        for r, c in [(row+1, col),
                     (row-1, col),
                     (row, col+1),
                     (row, col-1)]:
            self.grid[r][c] += qty
        return

    # display the grid using tkinter:
    def display(self):
        root = Tk()
        root.title('Sandpile')
        root.geometry('700x700+100+50')
        sp = Sandpile(root, self.grid)
        root.mainloop()

# execute program for size, center value pair:
# just print results for a small grid
g = Grid(9,17)
g.show('BEFORE')
g.abelian()          # scatter the sand
g.show('AFTER')

# just show results in tkinter for a large grid
# I wish there was a way to attach a screen shot
# of the tkinter result here
g = Grid(131,25000)
g.abelian()          # scatter the sand
g.display()          # display result using tkinter
  
##  OUTPUT:
##
##      BEFORE:
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 17 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##
##      AFTER:
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 1 0 0 0 0
##    0 0 0 2 1 2 0 0 0
##    0 0 1 1 1 1 1 0 0
##    0 0 0 2 1 2 0 0 0
##    0 0 0 0 1 0 0 0 0
##    0 0 0 0 0 0 0 0 0
##    0 0 0 0 0 0 0 0 0

R

# Return (x,y) index from a grid from an index in a list based on the grid size 
pos_to_index <- function(n) {
  f1 <- n/gridsize
  col <- ifelse(n%%gridsize == 0, f1,as.integer(f1)+1)
  row <- n - ((col-1)*gridsize)
  list(row=row,col=col)
}

# Return adjacent indexes (north, east, south, west)
adjacent_indexes <- function(r,c) {
  rup <- r - 1
  rdn <- ifelse(r == gridsize,0,r + 1)
  cleft <- c - 1
  cright <- ifelse(c==gridsize,0,c+1)
  list(up=c(rup,c),right=c(r,cright),left=c(r,cleft),down=c(rdn,c))
}

# Generate Abelian pattern
abelian <- function(gridsize,sand) {
  mat_ <- matrix(rep(0,gridsize^2),gridsize)
  midv <- as.integer(gridsize/2) + 1
  mat_[midv,midv] <- sand
  cat("Before\n")
  print(mat_)

  while(T) {
    cnt <- cnt + 1
    tgt <- which(mat_ >= 4)
    if (length(tgt) == 0) break
    pos   <- pos_to_index(tgt[1])
    idxes <- adjacent_indexes(pos$row,pos$col)
    mat_[pos$row,pos$col] <- mat_[pos$row,pos$col] - 4

    for (i in idxes) if (0 %in% i == F)  mat_[i[1],i[2]] <- mat_[i[1],i[2]] +1
  }
  cat("After\n")
  print(mat_) 
}

# Main

abelian(10,64)

Output:

Before
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    0    0    0    0    0    0     0
 [4,]    0    0    0    0    0    0    0    0    0     0
 [5,]    0    0    0    0    0    0    0    0    0     0
 [6,]    0    0    0    0    0   64    0    0    0     0
 [7,]    0    0    0    0    0    0    0    0    0     0
 [8,]    0    0    0    0    0    0    0    0    0     0
 [9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0
After
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    0    1    2    1    0    0     0
 [4,]    0    0    0    2    2    2    2    2    0     0
 [5,]    0    0    1    2    2    2    2    2    1     0
 [6,]    0    0    2    2    2    0    2    2    2     0
 [7,]    0    0    1    2    2    2    2    2    1     0
 [8,]    0    0    0    2    2    2    2    2    0     0
 [9,]    0    0    0    0    1    2    1    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

Raku

(formerly Perl 6)

Terminal based

Works with: Rakudo version 2019.07.1

Defaults to a stack of 1000 and showing progress. Pass in a custom stack size if desired and -hide-progress to run without displaying progress (much faster.)

sub cleanup { print "\e[0m\e[?25h\n"; exit(0) }

signal(SIGINT).tap: { cleanup(); exit(0) }

unit sub MAIN ($stack = 1000, :$hide-progress = False );

my @color = "\e[38;2;0;0;0m█",
            "\e[38;2;255;0;0m█",
            "\e[38;2;255;255;0m█",
            "\e[38;2;0;0;255m█",
            "\e[38;2;255;255;255m█"
            ;

my ($h, $w) = qx/stty size/.words».Int;
my $buf = $w * $h;
my @buffer = 0 xx $buf;
my $done;

@buffer[$w * ($h div 2) + ($w div 2) - 1] = $stack;

print "\e[?25l\e[48;5;232m";

repeat {
    $done = True;
    loop (my int $row; $row < $h; $row = $row + 1) {
        my int $rs = $row * $w; # row start
        my int $re = $rs  + $w; # row end
        loop (my int $idx = $rs; $idx < $re; $idx = $idx + 1) {
            if @buffer[$idx] >= 4 {
                my $grains = @buffer[$idx] div 4;
                @buffer[ $idx - $w ] += $grains if $row > 0;
                @buffer[ $idx - 1  ] += $grains if $idx - 1 >= $rs;
                @buffer[ $idx + $w ] += $grains if $row < $h - 1;
                @buffer[ $idx + 1  ] += $grains if $idx + 1 < $re;
                @buffer[ $idx ] %= 4;
                $done = False;
            }
        }
    }
    unless $hide-progress {
        print "\e[1;1H", @buffer.map( { @color[$_ min 4] }).join;
    }
} until $done;

print "\e[1;1H", @buffer.map( { @color[$_ min 4] }).join;

cleanup;

Passing in 2048 as a stack size results in: Abelian-sandpile-model-perl6.png (offsite .png image)

SDL2 Animation

use NativeCall;
use SDL2::Raw;

my ($width, $height) = 900, 900;

unit sub MAIN ($stack = 10000);

my int ($w, $h) = 160, 160;

my $buf = $w * $h;
my @buffer = 0 xx $buf;

@buffer[$w * ($h div 2) + ($w div 2) - 1] = $stack;


SDL_Init(VIDEO);

my SDL_Window $window = SDL_CreateWindow(
    "Abelian sandpile - Raku",
    SDL_WINDOWPOS_CENTERED_MASK, SDL_WINDOWPOS_CENTERED_MASK,
    $width, $height,
    RESIZABLE
);

my SDL_Renderer $renderer = SDL_CreateRenderer( $window, -1, ACCELERATED +| TARGETTEXTURE );

my $asp_texture = SDL_CreateTexture($renderer, %PIXELFORMAT<RGB332>, STREAMING, $w, $h);

my $pixdatabuf = CArray[int64].new(0, $w, $h, $w);

my @color = 0x00, 0xDE, 0x14, 0xAA, 0xFF;

sub render {
    my int $pitch;
    my int $cursor;

    # work-around to pass the pointer-pointer.
    my $pixdata = nativecast(Pointer[int64], $pixdatabuf);
    SDL_LockTexture($asp_texture, SDL_Rect, $pixdata, $pitch);

    $pixdata = nativecast(CArray[int8], Pointer.new($pixdatabuf[0]));

    loop (my int $row; $row < $h; $row = $row + 1) {
        my int $rs = $row * $w; # row start
        my int $re = $rs  + $w; # row end
        loop (my int $idx = $rs; $idx < $re; $idx = $idx + 1) {
            $pixdata[$idx] =  @buffer[$idx] < 4 ?? @color[@buffer[$idx]] !! @color[4];
            if @buffer[$idx] >= 4 {
                my $grains = floor @buffer[$idx] / 4;
                @buffer[ $idx - $w ] += $grains if $row > 0;
                @buffer[ $idx - 1  ] += $grains if $idx - 1 >= $rs;
                @buffer[ $idx + $w ] += $grains if $row < $h - 1;
                @buffer[ $idx + 1  ] += $grains if $idx + 1 < $re;
                @buffer[ $idx ] %= 4;
            }
        }
    }

    SDL_UnlockTexture($asp_texture);

    SDL_RenderCopy($renderer, $asp_texture, SDL_Rect, SDL_Rect.new(:x(0), :y(0), :w($width), :h($height)));
    SDL_RenderPresent($renderer);
}

my $event = SDL_Event.new;

main: loop {

    while SDL_PollEvent($event) {
        my $casted_event = SDL_CastEvent($event);

        given $casted_event {
            when *.type == QUIT {
                last main;
            }
        }
    }

    render();
    print fps;
}

say '';

sub fps {
    state $fps-frames = 0;
    state $fps-now    = now;
    state $fps        = '';
    $fps-frames++;
    if now - $fps-now >= 1 {
        $fps = [~] "\b" x 40, ' ' x 20, "\b" x 20 ,
            sprintf "FPS: %5.2f  ", ($fps-frames / (now - $fps-now)).round(.01);
        $fps-frames = 0;
        $fps-now = now;
    }
    $fps
}

Passing in a stack size of 20000 results in: Abelian-sandpile-sdl2.png (offsite .png image)

RPL

Using the built-in matrix data structure fulfils the requirements of the task:

[[1 2 0][2 1 1][0 1 3]] [[2 1 3][1 0 1][0 1 0]] +

Output:

1: [[3 3 3]
    [3 1 2]
    [0 2 3]]
RPL code Comment
≪ ROT OVER RE + { } + ROT ROT IM + +
≫ ‘→IDX’ STO

≪ DUP SIZE 1 GET → n
   ≪ DO 
        1 CF 1 n FOR h 1 n FOR j
           IF DUP h j 2 →LIST GET 3 > THEN
              1 SF DUP 0 CON
              h j 2 →LIST -4 PUT
              1 4 FOR a
                 h j (0,1) a ^ →IDX
                 IFERR 1 PUT THEN DROP2 END
              NEXT +
           END NEXT NEXT
      UNTIL 1 FC? END
≫ ≫ ‘SPILE’ STO
→IDX ( a b (c,d) → { a+c b+d } ) 


SPILE ( aa ) 
loop
  for h, j = 1 to n
    if a[h,j] > 3 then
       set flag, create empty matrix b
       b[h,j] = -4
       for a = 1 to 4
         (x,y) = (h,j) + i^a
          b[x,y] = 1 only if x > 0 and y > 0
       a += b
    end if, next h, j
until all elements <= 3
return a

It is sometimes necessary to run the program several times to reach stability: user's eye is much faster than a program to detect a remaining unstable sandpile. This is the way in RPL. It may nevertheless make sense when working on large matrices to have to run the program only once. In this case, the addtional line below shall be inserted after the END NEXT NEXT line:

1 n FOR h 1 n FOR j IF DUP h j 2 →LIST GET 3 > THEN 1 SF END NEXT NEXT
{3 3} 3 CON 'S3' STO
[[2 1 2][1 0 1][2 1 2]] 'S3ID' STO
S3 S3ID + SPILE SPILE
S3ID DUP + SPILE SPILE
Output:
2: [[ 3 3 3 ]
    [ 3 3 3 ]
    [ 3 3 3 ]]
1: [[ 2 1 2 ]
    [ 1 0 1 ]
    [ 2 1 2 ]]

Rust

// This is the main algorithm.
//
// It loops over the current state of the sandpile and updates it on-the-fly.
fn advance(field: &mut Vec<Vec<usize>>, boundary: &mut [usize; 4]) -> bool
{ 
    // This variable is used to check whether we changed anything in the array. If no, the loop terminates.
    let mut done = false;

    for y in boundary[0]..boundary[2]
    {
        for x in boundary[1]..boundary[3]
        {
            if field[y][x] >= 4
            {
                // This part was heavily inspired by the Pascal version. We subtract 4 as many times as we can
                // and distribute it to the neighbors. Also, in case we have outgrown the current boundary, we
                // update it to once again contain the entire sandpile.

                // The amount that gets added to the neighbors is the amount here divided by four and (implicitly) floored.
                // The remaining sand is just current modulo 4.
                let rem: usize = field[y][x] / 4;
                field[y][x] %= 4;

                // The isize casts are necessary because usize can not go below 0.
                // Also, the reason why x and y are compared to boundary[2]-1 and boundary[3]-1 is because for loops in
                // Rust are upper bound exclusive. This means a loop like 0..5 will only go over 0,1,2,3 and 4.
                if y as isize - 1 >= 0 {field[y-1][x] += rem; if y == boundary[0] {boundary[0]-=1;}}
                if x as isize - 1 >= 0 {field[y][x-1] += rem; if x == boundary[1] {boundary[1]-=1;}}
                if y+1 < field.len() {field[y+1][x] += rem; if x == boundary[2]-1 {boundary[2]+=1;}}
                if x+1 < field.len() {field[y][x+1] += rem; if y == boundary[3]-1 {boundary[3]+=1;}}

                done = true;
            }
        }
    }

    done
}

// This function can be used to display the sandpile in the console window.
//
// Each row is mapped onto chars and those characters are then collected into a string.
// These are then printed to the console.
//
// Eg.: [0,1,1,2,3,0] -> [' ','░','░','▒','▓',' ']-> " ░░▒▓ "
fn display(field: &Vec<Vec<usize>>)
{
    for row in field
    {
        let char_row = {
            row.iter().map(|c| {match c {
                0 => ' ',
                1 => '░',
                2 => '▒',
                3 => '▓',
                _ => '█'
            }}).collect::<String>()
        };
        println!("{}", char_row);
    }
}

// This function writes the end result to a file called "output.ppm".
//
// PPM is a very simple image format, however, it entirely uncompressed which leads to huge image sizes.
// Even so, for demonstrative purposes it's perfectly good to use. For something more robust, look into PNG libraries.
//
// Read more about the format here: http://netpbm.sourceforge.net/doc/ppm.html
fn write_pile(pile: &Vec<Vec<usize>>) {
    use std::fs::File;
    use std::io::Write;

    // We first create the file (or erase its contents if it already existed).
    let mut file = File::create("./output.ppm").unwrap();

    // Then we add the image signature, which is "P3 <newline>[width of image] [height of image]<newline>[maximum value of color]<newline>".
    write!(file, "P3\n{} {}\n255\n", pile.len(), pile.len()).unwrap();

    for row in pile {
        // For each row, we create a new string which has more or less enough capacity to hold the entire row.
        // This is for performance purposes, but shouldn't really matter much.
        let mut line = String::with_capacity(row.len() * 14);

        // We map each value in the field to a color.
        // These are just simple RGB values, 0 being the background, the rest being the "sand" itself.
        for elem in row {
            line.push_str(match elem {
                0 => "100 40 15 ",
                1 => "117 87 30 ",
                2 => "181 134 47 ",
                3 => "245 182 66 ",
                _ => unreachable!(),
            });
        }

        // Finally we write this string into the file.
        write!(file, "{}\n", line).unwrap();
    }
}

fn main() {
    // This is how big the final image will be. Currently the end result would be a 16x16 picture.
    let field_size = 16;
    let mut playfield = vec![vec![0; field_size]; field_size];

    // We put the initial sand in the exact middle of the field.
    // This isn't necessary per se, but it ensures that sand can fully topple.
    //
    // The boundary is initially just the single tile which has the sand in it, however, as the algorithm
    // progresses, this will grow larger too.
    let mut boundary = [field_size/2-1, field_size/2-1, field_size/2, field_size/2];
    playfield[field_size/2 - 1][field_size/2 - 1] = 16;

    // This is the main loop. We update the field until it returns false, signalling that the pile reached its
    // final state.
    while advance(&mut playfield, &mut boundary) {};

    // Once this happens, we simply display the result. Uncomment the line below to write it to a file.
    // Calling display with large field sizes is not recommended as it can easily become too large for the console.
    display(&playfield);
    //write_pile(&playfield);
}

Output:

                
                
                
                
                
       ░        
      ▒░▒       
     ░░ ░░      
      ▒░▒       
       ░        
                
                
                
                
                
                

Scheme

Works with: Chez Scheme
; A two-dimensional grid of values...

; Create an empty (all cells 0) grid of the specified size.
; Optionally, fill all cells with given value.
(define make-grid
  (lambda (size-x size-y . opt-value)
    (cons size-x (make-vector (* size-x size-y) (if (null? opt-value) 0 (car opt-value))))))

; Return the vector of all values of a grid.
(define grid-vector
  (lambda (grid)
    (cdr grid)))

; Return the X size of a grid.
(define grid-size-x
  (lambda (grid)
    (car grid)))

; Return the Y size of a grid.
(define grid-size-y
  (lambda (grid)
    (/ (vector-length (cdr grid)) (car grid))))

; Return #t if the specified x/y is within the range of the given grid.
(define grid-in-range
  (lambda (grid x y)
    (and (>= x 0) (>= y 0) (< x (grid-size-x grid)) (< y (grid-size-y grid)))))

; Return the value from the specified cell of the given grid.
; Note:  Returns 0 if x/y is out of range.
(define grid-ref
  (lambda (grid x y)
    (if (grid-in-range grid x y)
      (vector-ref (cdr grid) (+ x (* y (car grid))))
      0)))

; Store the given value into the specified cell of the given grid.
; Note:  Does nothing if x/y is out of range.
(define grid-set!
  (lambda (grid x y val)
    (when (grid-in-range grid x y)
      (vector-set! (cdr grid) (+ x (* y (car grid))) val))))

; Display the given grid, leaving correct spacing for maximum value.
; Optionally, uses a specified digit count for spacing.
; Returns the digit count of the largest grid value.
; Note:  Assumes the values in the grid are all non-negative integers.
(define grid-display
  (lambda (grid . opt-digcnt)
    ; Return count of digits in printed representation of integer.
    (define digit-count
      (lambda (int)
        (if (= int 0) 1 (1+ (exact (floor (log int 10)))))))
    ; Display the grid, leaving correct spacing for maximum value.
    (let* ((maxval (fold-left max 0 (vector->list (grid-vector grid))))
           (digcnt (if (null? opt-digcnt) (digit-count maxval) (car opt-digcnt))))
      (do ((y 0 (1+ y)))
          ((>= y (grid-size-y grid)))
        (do ((x 0 (1+ x)))
            ((>= x (grid-size-x grid)))
          (printf " ~vd" digcnt (grid-ref grid x y)))
        (printf "~%"))
      digcnt)))

; Implementation of the Abelian Sandpile Model using the above grid...

; Topple the specified cell of the given Abelian Sandpile Model grid.
; If number of grains in cell is less than 4, does nothing and returns #f.
; Otherwise, distributes 4 grains from the cell to its nearest neighbors and returns #t.
(define asm-topple
  (lambda (asm x y)
    (if (< (grid-ref asm x y) 4)
      #f
      (begin
        (grid-set! asm x y (- (grid-ref asm x y) 4))
        (grid-set! asm (1- x) y (1+ (grid-ref asm (1- x) y)))
        (grid-set! asm (1+ x) y (1+ (grid-ref asm (1+ x) y)))
        (grid-set! asm x (1- y) (1+ (grid-ref asm x (1- y))))
        (grid-set! asm x (1+ y) (1+ (grid-ref asm x (1+ y))))
        #t))))

; Repeatedly topple unstable cells in the given Abelian Sandpile Model grid
; until all cells are stable.
(define asm-stabilize
  (lambda (asm)
    (let loop ((any-toppled #f))
      (do ((y 0 (1+ y)))
          ((>= y (grid-size-y asm)))
        (do ((x 0 (1+ x)))
            ((>= x (grid-size-x asm)))
          (when (asm-topple asm x y)
            (set! any-toppled #t))))
      (when any-toppled
        (loop #f)))))

; Test the Abelian Sandpile Model on a simple grid...

(let ((asm (make-grid 9 9)))
  (grid-set! asm 4 4 64)
  (printf "Before:~%")
  (let ((digcnt (grid-display asm)))
    (asm-stabilize asm)
    (printf "After:~%")
    (grid-display asm digcnt)))
Output:
Before:
  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0
  0  0  0  0 64  0  0  0  0
  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0
  0  0  0  0  0  0  0  0  0
After:
  0  0  0  0  0  0  0  0  0
  0  0  0  1  2  1  0  0  0
  0  0  2  2  2  2  2  0  0
  0  1  2  2  2  2  2  1  0
  0  2  2  2  0  2  2  2  0
  0  1  2  2  2  2  2  1  0
  0  0  2  2  2  2  2  0  0
  0  0  0  1  2  1  0  0  0
  0  0  0  0  0  0  0  0  0

VBA

Sub SetupPile(a As Integer, b As Integer)
Application.ScreenUpdating = False
For i = 1 To a
For j = 1 To b
Cells(i, j).value = ""
Cells(i, j).Select

With Selection.Borders(xlEdgeLeft)
    .LineStyle = xlContinuous
    .Weight = xlMedium
End With
With Selection.Borders(xlEdgeTop)
    .LineStyle = xlContinuous
    .Weight = xlMedium
End With
With Selection.Borders(xlEdgeBottom)
    .LineStyle = xlContinuous
    .Weight = xlMedium
End With
With Selection.Borders(xlEdgeRight)
    .LineStyle = xlContinuous
    .Weight = xlMedium
End With

With Selection
    .HorizontalAlignment = xlCenter
    .VerticalAlignment = xlCenter
End With

Next j
Next i
Application.ScreenUpdating = True
End Sub


Sub Abelian_Sandpile()
Dim PileWidth As Integer
Dim PileHeight As Integer
Dim FieldArray() As Integer

Debug.Print "Start:" & Now()

'Set Size of Playing Field
PileWidth = 25
PileHeight = 25

ReDim FieldArray(PileWidth - 1, PileHeight - 1)

'Paint Basic Grid
SetupPile PileWidth, PileHeight

'Drop sand amount into middle of playing field
SandDropAmount = 1000
'Get around excel's incorrect rounding
SandDropColumn = Round((PileWidth / 2) + 0.001, 0)
SandDropRow = Round((PileHeight / 2) + 0.001, 0)

Cells(SandDropRow, SandDropColumn) = SandDropAmount
FieldArray(SandDropRow - 1, SandDropColumn - 1) = SandDropAmount

Continue = False

'Check if Pile is already stabilized at the start
For i = 1 To PileWidth 'Col
For j = 1 To PileHeight 'Row
If FieldArray(j - 1, i - 1) > 3 Then Continue = True
Next j
Next i

'While not stabilized
While Continue
For i = 1 To PileWidth
For j = 1 To PileHeight
    If FieldArray(j - 1, i - 1) > 3 Then
    'Reduce by 4
    FieldArray(j - 1, i - 1) = FieldArray(j - 1, i - 1) - 4
    'Increase Neighbours
    't
    If j >= 2 Then FieldArray(j - 2, i - 1) = FieldArray(j - 2, i - 1) + 1
    'r
    If i < PileWidth Then FieldArray(j - 1, i) = FieldArray(j - 1, i) + 1
    'b
    If j < PileHeight Then FieldArray(j, i - 1) = FieldArray(j, i - 1) + 1
    'l
    If i >= 2 Then FieldArray(j - 1, i - 2) = FieldArray(j - 1, i - 2) + 1
    'Next round
    GoTo Nextone
    End If
Next j
Next i

Nextone:

'Check if now stabilized
Continue = False
For i = 1 To PileWidth
For j = 1 To PileHeight
'Paint every step if needed
'Cells(j, i) = FieldArray(j - 1, i - 1)

If FieldArray(j - 1, i - 1) > 3 Then Continue = True
Next j
Next i

Wend

'Print out final step
For i = 1 To PileWidth
For j = 1 To PileHeight
Cells(j, i) = FieldArray(j - 1, i - 1)
Next j
Next i

'Make field square and remove 0
Cells.Select
Selection.ColumnWidth = 2
Selection.RowHeight = 13.5
Selection.Replace What:="0", Replacement:="", LookAt:=xlPart, SearchOrder:=xlByRows, MatchCase:=False, SearchFormat:=False, ReplaceFormat:=False
Range("A1").Select

Range(Cells(1, 1), Cells(PileHeight, PileWidth)).Select

'Conditional Format
Selection.FormatConditions.AddColorScale ColorScaleType:=3
Selection.FormatConditions(Selection.FormatConditions.Count).SetFirstPriority
Selection.FormatConditions(1).ColorScaleCriteria(1).Type = xlConditionValueLowestValue
With Selection.FormatConditions(1).ColorScaleCriteria(1).FormatColor
    .Color = 8109667
    .TintAndShade = 0
End With
Selection.FormatConditions(1).ColorScaleCriteria(2).Type = xlConditionValuePercentile
Selection.FormatConditions(1).ColorScaleCriteria(2).value = 50
With Selection.FormatConditions(1).ColorScaleCriteria(2).FormatColor
    .Color = 8711167
    .TintAndShade = 0
End With
Selection.FormatConditions(1).ColorScaleCriteria(3).Type = xlConditionValueHighestValue
With Selection.FormatConditions(1).ColorScaleCriteria(3).FormatColor
    .Color = 7039480
    .TintAndShade = 0
End With
Range("A1").Select

Debug.Print "W,H,A:" & PileWidth & "," & PileHeight & "," & SandDropAmount
Debug.Print "End:" & Now()

End Sub

Output:

On Excel Page

V (Vlang)

Translation of: Go
import os
import strings

const dim = 16

// Outputs the result to the terminal using UTF-8 block characters.
fn draw_pile(pile [][]int) {
    chars:= [` `,`░`,`▓`,`█`]
    for row in pile {
        mut line := []rune{len: row.len}
        for i, e in row {
            mut elem := e
            if elem > 3 { // only possible when algorithm not yet completed.
                elem = 3
            }
            line[i] = chars[elem]
        }
        println(line.string())
    }
}
 
// Creates a .ppm file in the current directory, which contains
// a colored image of the pile.
fn write_pile(pile [][]int) {
    mut file := os.create("output.ppm") or {panic('ERROR creating file')}
    defer {
file.close()
}
    // Write the signature, image dimensions and maximum color value to the file.
    file.writeln("P3\n$dim $dim\n255") or {panic('ERROR writing ln')}
    bcolors := ["125 0 25 ", "125 80 0 ", "186 118 0 ", "224 142 0 "]
    mut line := strings.new_builder(128)
    for row in pile {        
        for elem in row {
            line.write_string(bcolors[elem])
        }
        file.write_string('${line.str()}\n') or {panic('ERROR writing str')}
        line = strings.new_builder(128)
    }
}
 
// Main part of the algorithm, a simple, recursive implementation of the model.
fn handle_pile(x int, y int, mut pile [][]int) {
    if pile[y][x] >= 4 {
        pile[y][x] -= 4
        // Check each neighbor, whether they have enough "sand" to collapse and if they do,
        // recursively call handle_pile on them.
        if y > 0 {
            pile[y-1][x]++
            if pile[y-1][x] >= 4 {
                handle_pile(x, y-1, mut pile)
            }
        }
        if x > 0 {
            pile[y][x-1]++
            if pile[y][x-1] >= 4 {
                handle_pile(x-1, y, mut pile)
            }
        }
        if y < dim-1 {
            pile[y+1][x]++
            if pile[y+1][x] >= 4 {
                handle_pile(x, y+1, mut pile)
            }
        }
        if x < dim-1 {
            pile[y][x+1]++
            if pile[y][x+1] >= 4 {
                handle_pile(x+1, y, mut pile)
            }
        }
 
        // Uncomment this line to show every iteration of the program.
        // Not recommended with large input values.
        // draw_pile(pile)
 
        // Finally call the fntion on the current cell again,
        // in case it had more than 4 particles.
        handle_pile(x, y, mut pile)
    }
}
 
fn main() {
    // Create 2D grid and set size using the 'dim' constant.
    mut pile := [][]int{len: dim, init: []int{len: dim}}
 
    // Place some sand particles in the center of the grid and start the algorithm.
    hdim := int(dim/2 - 1)
    pile[hdim][hdim] = 16
    handle_pile(hdim, hdim, mut pile)
    draw_pile(pile)
 
    // Uncomment this to save the final image to a file
    // after the recursive algorithm has ended.
    // write_pile(pile)
}
Output:
                
                
                
                
                
       ░        
      ▓░▓       
     ░░ ░░      
      ▓░▓       
       ░        
                
                
                
                
                
                
       

Wren

Library: Wren-fmt
import "./fmt" for Fmt

class Sandpile {
    // 'a' is a list of integers in row order
    construct new(a) {
        var count = a.count
        _rows = count.sqrt.floor
        if (_rows * _rows != count) Fiber.abort("The matrix of values must be square.")
        _a = a
        _neighbors = List.filled(count, 0)
        for (i in 0...count) {
            _neighbors[i] = []
            if (i % _rows > 0)     _neighbors[i].add(i-1)
            if ((i + 1)%_rows > 0) _neighbors[i].add(i+1)
            if (i - _rows >= 0)    _neighbors[i].add(i-_rows)
            if (i + _rows < count) _neighbors[i].add(i+_rows)
        }
    }

    isStable { _a.all { |i| i <= 3 } }

    // topples until stable
    topple() {
        while (!isStable) {
            for (i in 0..._a.count) {
                if (_a[i] > 3) {
                    _a[i] = _a[i] - 4
                    for (j in _neighbors[i]) _a[j] = _a[j] + 1
                }
            }
        }
    }

    toString {
        var s = ""
        for (i in 0..._rows) {
            for (j in 0..._rows) s = s + Fmt.swrite("$2d ", _a[_rows*i + j])
            s = s + "\n"
        }
        return s
    }
}

var printAcross = Fn.new { |str1, str2|
    var r1 = str1.split("\n")
    var r2 = str2.split("\n")
    var rows = r1.count - 1
    var cr = (rows/2).floor
    for (i in 0...rows) {
        var symbol = (i == cr) ? "->" : "  "
        Fmt.print("$s $s $s", r1[i], symbol, r2[i])
    }
    System.print()
}

var a1 = List.filled(25, 0)
a1[12] = 4
var a2 = List.filled(25, 0)
a2[12] = 6
var a3 = List.filled(25, 0)
a3[12] = 16
var a4 = List.filled(100, 0)
a4[55] = 64
for (a in [a1, a2, a3, a4]) {
    var s = Sandpile.new(a)
    var str1 = s.toString
    s.topple()
    var str2 = s.toString
    printAcross.call(str1, str2)
}
Output:
 0  0  0  0  0      0  0  0  0  0 
 0  0  0  0  0      0  0  1  0  0 
 0  0  4  0  0  ->  0  1  0  1  0 
 0  0  0  0  0      0  0  1  0  0 
 0  0  0  0  0      0  0  0  0  0 

 0  0  0  0  0      0  0  0  0  0 
 0  0  0  0  0      0  0  1  0  0 
 0  0  6  0  0  ->  0  1  2  1  0 
 0  0  0  0  0      0  0  1  0  0 
 0  0  0  0  0      0  0  0  0  0 

 0  0  0  0  0      0  0  1  0  0 
 0  0  0  0  0      0  2  1  2  0 
 0  0 16  0  0  ->  1  1  0  1  1 
 0  0  0  0  0      0  2  1  2  0 
 0  0  0  0  0      0  0  1  0  0 

 0  0  0  0  0  0  0  0  0  0      0  0  0  0  0  0  0  0  0  0 
 0  0  0  0  0  0  0  0  0  0      0  0  0  0  0  0  0  0  0  0 
 0  0  0  0  0  0  0  0  0  0      0  0  0  0  1  2  1  0  0  0 
 0  0  0  0  0  0  0  0  0  0      0  0  0  2  2  2  2  2  0  0 
 0  0  0  0  0  0  0  0  0  0      0  0  1  2  2  2  2  2  1  0 
 0  0  0  0  0 64  0  0  0  0  ->  0  0  2  2  2  0  2  2  2  0 
 0  0  0  0  0  0  0  0  0  0      0  0  1  2  2  2  2  2  1  0 
 0  0  0  0  0  0  0  0  0  0      0  0  0  2  2  2  2  2  0  0 
 0  0  0  0  0  0  0  0  0  0      0  0  0  0  1  2  1  0  0  0 
 0  0  0  0  0  0  0  0  0  0      0  0  0  0  0  0  0  0  0  0 

XPL0

def  Size = 200;
char Pile1(Size*Size), Pile2(Size*Size);
int  Spigot, I, X, Y;
[SetVid($13);   \VGA 320x200x8
FillMem(Pile1, 0, Size*Size);
FillMem(Pile2, 0, Size*Size);
Spigot:= 400_000;
repeat  I:= 0;
        for Y:= 0 to Size-1 do
            for X:= 0 to Size-1 do
                [if X=Size/2 & Y=Size/2 then
                    [Spigot:= Spigot-4;
                     Pile1(I):= Pile1(I)+4;
                    ];
                if Pile1(I) >= 4 then
                    [Pile1(I):= Pile1(I)-4;
                     Pile2(I-1):= Pile2(I-1)+1;
                     Pile2(I+1):= Pile2(I+1)+1;
                     Pile2(I-Size):= Pile2(I-Size)+1;
                     Pile2(I+Size):= Pile2(I+Size)+1;
                    ];
                Point(X, Y, Pile1(I)*2);
                I:= I+1;
                ];
        I:= Pile1;  Pile1:= Pile2;  Pile2:= I;
until   Spigot < 4;
]
Output:
[http://www.xpl0.org/rpi/sand.gif]