Abelian sandpile model/Identity

Abelian sandpile model/Identity
You are encouraged to solve this task according to the task description, using any language you may know.

Our sandpiles are based on a 3 by 3 rectangular grid giving nine areas that contain a number from 0 to 3 inclusive. (The numbers are said to represent grains of sand in each area of the sandpile).

E.g. s1 =

    
    1 2 0
    2 1 1
    0 1 3

and s2 =

    2 1 3
    1 0 1
    0 1 0

Addition on sandpiles is done by adding numbers in corresponding grid areas, so for the above:

              1 2 0     2 1 3     3 3 3
    s1 + s2 = 2 1 1  +  1 0 1  =  3 1 2
              0 1 3     0 1 0     0 2 3

If the addition would result in more than 3 "grains of sand" in any area then those areas cause the whole sandpile to become "unstable" and the sandpile areas are "toppled" in an "avalanche" until the "stable" result is obtained.

Any unstable area (with a number >= 4), is "toppled" by loosing one grain of sand to each of its four horizontal or vertical neighbours. Grains are lost at the edge of the grid, but otherwise increase the number in neighbouring cells by one, whilst decreasing the count in the toppled cell by four in each toppling.

A toppling may give an adjacent area more than four grains of sand leading to a chain of topplings called an "avalanche". E.g.

    
    4 3 3     0 4 3     1 0 4     1 1 0     2 1 0
    3 1 2 ==> 4 1 2 ==> 4 2 2 ==> 4 2 3 ==> 0 3 3
    0 2 3     0 2 3     0 2 3     0 2 3     1 2 3

The final result is the stable sandpile on the right.

Note: The order in which cells are toppled does not affect the final result.

Task

• Create a class or datastructure and functions to represent and operate on sandpiles.
• Confirm the result of the avalanche of topplings shown above
• Confirm that s1 + s2 == s2 + s1 # Show the stable results
• If s3 is the sandpile with number 3 in every grid area, and s3_id is the following sandpile:

    2 1 2  
    1 0 1  
    2 1 2

• Show that s3 + s3_id == s3
• Show that s3_id + s3_id == s3_id

Show confirming output here, with your examples.

References

• https://www.youtube.com/watch?v=1MtEUErz7Gg
• https://en.wikipedia.org/wiki/Abelian_sandpile_model

11l

Translation of: Python

T Sandpile
   DefaultDict[(Int, Int), Int] grid

   F (gridtext)
      V array = gridtext.split_py().map(x -> Int(x))
      L(x) array
         .grid[(L.index I/ 3, L.index % 3)] = x

   Set[(Int, Int)] _border = Set(cart_product(-1 .< 4, -1 .< 4).filter((r, c) -> !(r C 0..2) | !(c C 0..2)))
   _cell_coords = cart_product(0.<3, 0.<3)

   F topple()
      V& g = .grid
      L(r, c) ._cell_coords
         I g[(r, c)] >= 4
            g[(r - 1, c)]++
            g[(r + 1, c)]++
            g[(r, c - 1)]++
            g[(r, c + 1)]++
            g[(r, c)] -= 4
            R 1B
      R 0B

   F stabilise()
      L .topple() {}

      L(row_col) ._border.intersection(Set(.grid.keys()))
         .grid.pop(row_col)

   F ==(other)
      R all(._cell_coords.map(row_col -> @.grid[row_col] == @other.grid[row_col]))

   F +(other)
      V ans = Sandpile(‘’)
      L(row_col) ._cell_coords
         ans.grid[row_col] = .grid[row_col] + other.grid[row_col]
      ans.stabilise()
      R ans

   F String()
      [String] txt
      L(row) 3
         txt.append((0.<3).map(col -> String(@.grid[(@row, col)])).join(‘ ’))
      R txt.join("\n")

V unstable = Sandpile(‘4 3 3
                       3 1 2
                       0 2 3’)
V s1 = Sandpile(‘1 2 0
                 2 1 1
                 0 1 3’)
V s2 = Sandpile(‘2 1 3
                 1 0 1
                 0 1 0’)
V s3 = Sandpile(‘3 3 3  3 3 3  3 3 3’)
V s3_id = Sandpile(‘2 1 2  1 0 1  2 1 2’)

print(unstable)
print()
unstable.stabilise()
print(unstable)
print()
print(s1 + s2)
print()
print(s2 + s1)
print()
print(s1 + s2 == s2 + s1)
print()
print(s3 + s3_id)
print()
print(s3 + s3_id == s3)
print()
print(s3_id + s3_id)
print()
print(s3_id + s3_id == s3_id)

Output:

4 3 3
3 1 2
0 2 3

2 1 0
0 3 3
1 2 3

3 3 3
3 1 2
0 2 3

3 3 3
3 1 2
0 2 3

1B

3 3 3
3 3 3
3 3 3

1B

2 1 2
1 0 1
2 1 2

1B

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 abelianSum64.s   */ 

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

/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessValue:        .asciz "@ "
szMessAdd1:         .asciz "Add sandpile 1 to sandpile 2  \n"
szMessAdd2:         .asciz "Add sandpile 2 to sandpile 1  \n"
szMessAdd2A:        .asciz "Add sandpile 2A to sandpile result  \n"
szMessAdd3:         .asciz "Add sandpile 3 to sandpile 3ID  \n"
szMessAdd3ID:       .asciz "Add sandpile 3ID to sandpile 3ID  \n"

szCarriageReturn:   .asciz "\n"

qSandPile1:    .quad 1,2,0
               .quad 2,1,1
               .quad 0,1,3
               
qSandPile2:    .quad 2,1,3
               .quad 1,0,1
               .quad 0,1,0

qSandPile2A:    .quad 1,0,0
               .quad 0,0,0
               .quad 0,0,0
               
qSandPile3:    .quad 3,3,3
               .quad 3,3,3
               .quad 3,3,3
               
qSandPile3ID:  .quad 2,1,2
               .quad 1,0,1
               .quad 2,1,2
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss
sZoneConv:        .skip 24
qSandPilex1:      .skip 8 * MAXI * MAXI
qSandPilex2:      .skip 8 * MAXI * MAXI
/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:                            // entry of program 

    ldr x0,qAdrqSandPile1        // sandpile1 address
    ldr x1,qAdrqSandPile2        // sandpile2 address
    ldr x2,qAdrqSandPilex1       // sandpile result address
    bl addSandPile
    
    ldr x0,qAdrszMessAdd1        // display message
    bl affichageMess
    ldr x0,qAdrqSandPilex1       // display sandpile
    bl displaySandPile
    
    ldr x0,qAdrqSandPile2        // sandpile2 address
    ldr x1,qAdrqSandPile1        // sandpile1 address
    ldr x2,qAdrqSandPilex1       // sandpile result address
    bl addSandPile
    
    ldr x0,qAdrszMessAdd2
    bl affichageMess
    ldr x0,qAdrqSandPilex1
    bl displaySandPile
    
    ldr x0,qAdrqSandPilex1        // sandpile1 address
    ldr x1,qAdrqSandPile2A        // sandpile2A address
    ldr x2,qAdrqSandPilex2        // sandpile result address
    bl addSandPile
    
    ldr x0,qAdrszMessAdd2A
    bl affichageMess
    ldr x0,qAdrqSandPilex2
    bl displaySandPile
    
    ldr x0,qAdrqSandPile3          // sandpile3 address
    ldr x1,qAdrqSandPile3ID        // sandpile3ID address
    ldr x2,qAdrqSandPilex2         // sandpile result address
    bl addSandPile
    
    ldr x0,qAdrszMessAdd3
    bl affichageMess
    ldr x0,qAdrqSandPilex2
    bl displaySandPile
    
    ldr x0,qAdrqSandPile3ID        // sandpile3 address
    ldr x1,qAdrqSandPile3ID        // sandpile3ID address
    ldr x2,qAdrqSandPilex2         // sandpile result address
    bl addSandPile
    
    ldr x0,qAdrszMessAdd3ID
    bl affichageMess
    ldr x0,qAdrqSandPilex2
    bl displaySandPile
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
qAdrszMessAdd1:           .quad szMessAdd1
qAdrszMessAdd2:           .quad szMessAdd2
qAdrszMessAdd2A:          .quad szMessAdd2A
qAdrszMessAdd3:           .quad szMessAdd3
qAdrszMessAdd3ID:         .quad szMessAdd3ID
qAdrqSandPile1:           .quad qSandPile1
qAdrqSandPilex1:          .quad qSandPilex1
qAdrqSandPilex2:          .quad qSandPilex2
qAdrqSandPile2:           .quad qSandPile2
qAdrqSandPile2A:          .quad qSandPile2A
qAdrqSandPile3:           .quad qSandPile3
qAdrqSandPile3ID:         .quad qSandPile3ID
/***************************************************/
/*     add two  sandpile               */
/***************************************************/
// x0 contains address to sandpile 1
// x1 contains address to sandpile 2
// x2 contains address to sandpile result 
addSandPile:
    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,x1                 // save addresse sandpile2
    mov x1,x2                 // and copy sandpile 1 to sandpile result
    bl copySandPile
    mov x0,x2                 // sanspile result
    mov x2,#0                 // indice y
    mov x4,#MAXI
1:
    mov x1,#0                  // indice x
2:
    madd x5,x2,x4,x1           // compute offset
    ldr x7,[x0,x5,lsl #3]      // load value at pos x,y sanspile result
    ldr x3,[x6,x5,lsl #3]      // load value at pos x,y sandpile 2
    add x7,x7,x3
    str x7,[x0,x5,lsl #3]      // store sum on sandpile result
    bl avalancheRisk
    add x1,x1,#1
    cmp x1,#MAXI
    blt 2b
    add x2,x2,#1
    cmp x2,#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
/***************************************************/
/*     copy sandpile                               */
/***************************************************/
// x0 contains address to sandpile 
// x1 contains address to sandpile result 
copySandPile:
    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 x2,#0                   // indice y
    mov x3,#MAXI
1:
    mov x4,#0                   // indice x
2:
    madd x5,x2,x3,x4            // compute offset
    ldr x6,[x0,x5,lsl #3]       // load value at pos x,y sanspile
    str x6,[x1,x5,lsl #3]       // store value at pos x,y sandpile result
    add x4,x4,#1
    cmp x4,#MAXI
    blt 2b
    add x2,x2,#1
    cmp x2,#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
/***************************************************/
/*     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
/***************************************************/
/*     avalanche risk              */
/***************************************************/
// x0 contains address to sanspile
// x1 contains position x
// x2 contains position y
avalancheRisk:
    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
    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 avalancheRisk                  // and compute new pile
    sub x1,x1,#1
2:
    cmp x1,#0                   // left position ok ?
    beq 3f
    sub x1,x1,#1
    bl add1Sand
    bl avalancheRisk
    add x1,x1,#1
3:
    cmp x2,#0                   // higt position ok ?
    beq 4f
    sub x2,x2,#1
    bl add1Sand
    bl avalancheRisk
    add x2,x2,#1
4:
    cmp x2,#MAXI-1               // low position ok ?
    beq 5f
    add x2,x2,#1
    bl add1Sand
    bl avalancheRisk
    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
    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 $ abelianSum64
Add sandpile 1 to sandpile 2
3 3 3
3 1 2
0 2 3
Add sandpile 2 to sandpile 1
3 3 3
3 1 2
0 2 3
Add sandpile 2A to sandpile result
2 1 0
0 3 3
1 2 3
Add sandpile 3 to sandpile 3ID
3 3 3
3 3 3
3 3 3
Add sandpile 3ID to sandpile 3ID
2 1 2
1 0 1
2 1 2

ARM Assembly

Works with: as version Raspberry Pi
or android 32 bits with application Termux

/* ARM assembly Raspberry PI  or android 32 bits */
/*  program abelianSum.s   */ 

/* 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, 3

/*********************************/
/* Initialized data              */
/*********************************/
.data
szMessValue:        .asciz "@ "
szMessAdd1:         .asciz "Add sandpile 1 to sandpile 2  \n"
szMessAdd2:         .asciz "Add sandpile 2 to sandpile 1  \n"
szMessAdd2A:        .asciz "Add sandpile 2A to sandpile result  \n"
szMessAdd3:         .asciz "Add sandpile 3 to sandpile 3ID  \n"
szMessAdd3ID:       .asciz "Add sandpile 3ID to sandpile 3ID  \n"

szMessFin:          .asciz "End display :\n"
szCarriageReturn:   .asciz "\n"

iSandPile1:    .int 1,2,0
               .int 2,1,1
               .int 0,1,3
               
iSandPile2:    .int 2,1,3
               .int 1,0,1
               .int 0,1,0

iSandPile2A:    .int 1,0,0
               .int 0,0,0
               .int 0,0,0
               
iSandPile3:    .int 3,3,3
               .int 3,3,3
               .int 3,3,3
               
iSandPile3ID:  .int 2,1,2
               .int 1,0,1
               .int 2,1,2
/*********************************/
/* UnInitialized data            */
/*********************************/
.bss
sZoneConv:        .skip 24
iSandPileR1:      .skip 4 * MAXI * MAXI
iSandPileR2:      .skip 4 * MAXI * MAXI
/*********************************/
/*  code section                 */
/*********************************/
.text
.global main 
main:                            @ entry of program 

    ldr r0,iAdriSandPile1        @ sandpile1 address
    ldr r1,iAdriSandPile2        @ sandpile2 address
    ldr r2,iAdriSandPileR1       @ sandpile result address
    bl addSandPile
    
    ldr r0,iAdrszMessAdd1        @ display message
    bl affichageMess
    ldr r0,iAdriSandPileR1       @ display sandpile
    bl displaySandPile
    
    ldr r0,iAdriSandPile2        @ sandpile2 address
    ldr r1,iAdriSandPile1        @ sandpile1 address
    ldr r2,iAdriSandPileR1       @ sandpile result address
    bl addSandPile
    
    ldr r0,iAdrszMessAdd2
    bl affichageMess
    ldr r0,iAdriSandPileR1
    bl displaySandPile
    
    ldr r0,iAdriSandPileR1        @ sandpile1 address
    ldr r1,iAdriSandPile2A        @ sandpile2A address
    ldr r2,iAdriSandPileR2        @ sandpile result address
    bl addSandPile
    
    ldr r0,iAdrszMessAdd2A
    bl affichageMess
    ldr r0,iAdriSandPileR2
    bl displaySandPile
    
    ldr r0,iAdriSandPile3          @ sandpile3 address
    ldr r1,iAdriSandPile3ID        @ sandpile3ID address
    ldr r2,iAdriSandPileR2         @ sandpile result address
    bl addSandPile
    
    ldr r0,iAdrszMessAdd3
    bl affichageMess
    ldr r0,iAdriSandPileR2
    bl displaySandPile
    
    ldr r0,iAdriSandPile3ID        @ sandpile3 address
    ldr r1,iAdriSandPile3ID        @ sandpile3ID address
    ldr r2,iAdriSandPileR2         @ sandpile result address
    bl addSandPile
    
    ldr r0,iAdrszMessAdd3ID
    bl affichageMess
    ldr r0,iAdriSandPileR2
    bl displaySandPile
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
iAdrszMessFin:            .int szMessFin
iAdrszMessAdd1:           .int szMessAdd1
iAdrszMessAdd2:           .int szMessAdd2
iAdrszMessAdd2A:          .int szMessAdd2A
iAdrszMessAdd3:           .int szMessAdd3
iAdrszMessAdd3ID:         .int szMessAdd3ID
iAdriSandPile1:           .int iSandPile1
iAdriSandPileR1:          .int iSandPileR1
iAdriSandPileR2:          .int iSandPileR2
iAdriSandPile2:           .int iSandPile2
iAdriSandPile2A:          .int iSandPile2A
iAdriSandPile3:           .int iSandPile3
iAdriSandPile3ID:         .int iSandPile3ID
/***************************************************/
/*     add two  sandpile               */
/***************************************************/
// r0 contains address to sandpile 1
// r1 contains address to sandpile 2
// r2 contains address to sandpile result 
addSandPile:
    push {r1-r7,lr}           @ save  registers 
    mov r6,r1                 @ save addresse sandpile2
    mov r1,r2                 @ and copy sandpile 1 to sandpile result
    bl copySandPile
    mov r0,r2                 @ sanspile result
    mov r2,#0                 @ indice y
    mov r4,#MAXI
1:
    mov r1,#0                  @ indice x
2:
    mla r5,r2,r4,r1            @ compute offset
    ldr r7,[r0,r5,lsl #2]      @ load value at pos x,y sanspile result
    ldr r3,[r6,r5,lsl #2]      @ load value at pos x,y sandpile 2
    add r7,r3
    str r7,[r0,r5,lsl #2]      @ store sum on sandpile result
    bl avalancheRisk
    add r1,r1,#1
    cmp r1,#MAXI
    blt 2b
    add r2,r2,#1
    cmp r2,#MAXI
    blt 1b
100:
    pop {r1-r7,lr}             @ restaur registers
    bx lr                      @ return
/***************************************************/
/*     copy sandpile                               */
/***************************************************/
// r0 contains address to sandpile 
// r1 contains address to sandpile result 
copySandPile:
    push {r1-r6,lr}           @ save  registers 
    mov r2,#0                 @ indice y
    mov r3,#MAXI
1:
    mov r4,#0                   @ indice x
2:
    mla r5,r2,r3,r4            @ compute offset
    ldr r6,[r0,r5,lsl #2]      @ load value at pos x,y sanspile
    str r6,[r1,r5,lsl #2]      @ store value at pos x,y sandpile result
    add r4,r4,#1
    cmp r4,#MAXI
    blt 2b
    add r2,r2,#1
    cmp r2,#MAXI
    blt 1b
100:
    pop {r1-r6,lr}             @ restaur registers
    bx lr                      @ return
/***************************************************/
/*     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
/***************************************************/
/*     avalanche risk              */
/***************************************************/
// r0 contains address to sanspile
// r1 contains position x
// r2 contains position y
avalancheRisk:
    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 avalancheRisk                  @ and compute new pile
    sub r1,#1
2:
    cmp r1,#0                   @ left position ok ?
    beq 3f
    sub r1,#1
    bl add1Sand
    bl avalancheRisk
    add r1,#1
3:
    cmp r2,#0                   @ higt position ok ?
    beq 4f
    sub r2,#1
    bl add1Sand
    bl avalancheRisk
    add r2,#1
4:
    cmp r2,#MAXI-1               @ low position ok ?
    beq 5f
    add r2,#1
    bl add1Sand
    bl avalancheRisk
    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:

Add sandpile 1 to sandpile 2
3  3  3
3  1  2
0  2  3
Add sandpile 2 to sandpile 1
3  3  3
3  1  2
0  2  3
Add sandpile 2A to sandpile result
2  1  0
0  3  3
1  2  3
Add sandpile 3 to sandpile 3ID
3  3  3
3  3  3
3  3  3
Add sandpile 3ID to sandpile 3ID
2  1  2
1  0  1
2  1  2

Ada

Translation of: C++

This Ada example works with Ada 2012 because of the use the aspect with Default_Component_Value.

The package specification for Abelian_Sandpile is:

-- Works with Ada 2012

package Abelian_Sandpile is
   Limit : constant Integer := 4;

   type Sandpile is array (0 .. 2, 0 .. 2) of Natural with
      Default_Component_Value => 0;

   procedure Stabalize (Pile : in out Sandpile);
   function Is_Stable (Pile : in Sandpile) return Boolean;
   procedure Topple (Pile : in out Sandpile);
   function "+" (Left, Right : Sandpile) return Sandpile;
   procedure Print(PIle : in Sandpile);

end Abelian_Sandpile;

The package body for Abelian_Sandpile is

with Ada.Text_Io; use Ada.Text_IO;

package body Abelian_Sandpile is

   ---------------
   -- Stabalize --
   ---------------

   procedure Stabalize (Pile : in out Sandpile) is
   begin
      while not Is_Stable(Pile) loop
         Topple(Pile);
      end loop;
   end Stabalize;

   ---------------
   -- Is_Stable --
   ---------------

   function Is_Stable (Pile : in Sandpile) return Boolean is
   begin
     return (for all E of Pile => E < Limit);
   end Is_Stable;

   ------------
   -- Topple --
   ------------

   procedure Topple (Pile : in out Sandpile) is
   begin
      outer:
      for Row in Pile'Range(1) loop
         for Col in Pile'Range(2) loop
            if Pile(Row, Col) >= Limit then
               Pile(Row, Col) := Pile(Row, Col) - Limit;
               if Row > 0 then
                  Pile(Row - 1, Col) := Pile(Row -1, Col) + 1;
               end if;
               if Row < Pile'Last(1) then
                  Pile(Row + 1, Col) := Pile(Row + 1, Col) + 1;
               end if;
               if Col > 0 then
                  Pile(Row, Col - 1) := Pile(Row, Col - 1) + 1;
               end if;
               if Col < Pile'Last(2) then
                  Pile(Row, Col + 1) := Pile(Row, Col + 1) + 1;
               end if;

               exit outer;
            end if;
         end loop;
      end loop outer;
   end Topple;

   ---------
   -- "+" --
   ---------

   function "+" (Left, Right : Sandpile) return Sandpile is
      Result : Sandpile;
   begin
      for I in Sandpile'Range(1) loop
         for J in Sandpile'Range(2) loop
            Result(I, J) := Left(I, J) + Right(I, J);
         end loop;
      end loop;
      Stabalize(Result);
      return Result;
   end "+";

   -----------
   -- Print --
   -----------

   procedure Print(Pile : in Sandpile) is
   begin
      for row in Pile'Range(1) loop
         for col in Pile'Range(2) loop
            Put(Integer'Image(Pile(row, col)));
         end loop;
         New_Line;
      end loop;
      New_Line;
   end Print;
end Abelian_Sandpile;

The main procedure performing the same tests as the C++ example is

with Ada.Text_IO;      use Ada.Text_IO;
with Abelian_Sandpile; use Abelian_Sandpile;

procedure Main is
   sp    : Sandpile := ((4, 3, 3), (3, 1, 2), (0, 2, 3));
   s1    : Sandpile := ((1, 2, 0), (2, 1, 1), (0, 1, 3));
   s2    : Sandpile := ((2, 1, 3), (1, 0, 1), (0, 1, 0));
   s3    : Sandpile := ((3, 3, 3), (3, 3, 3), (3, 3, 3));
   s3_id : Sandpile := ((2, 1, 2), (1, 0, 1), (2, 1, 2));
   sum1  : Sandpile := s1 + s2;
   sum2  : Sandpile := s2 + s1;
   sum3  : Sandpile := s3 + s3_id;
   sum4  : Sandpile := s3_id + s3_id;

begin
   Put_Line ("Avalanche:");
   while not Is_Stable (sp) loop
      Print (sp);
      Put_Line ("stable? " & Boolean'Image (Is_Stable (sp)));
      New_Line;
      Topple (sp);
   end loop;
   Print (sp);
   Put_Line ("stable? " & Boolean'Image (Is_Stable (sp)));
   New_Line;

   Put_Line ("Commutivity:");
   Put_Line ("s1 + s2 equals s2 + s2? " & Boolean'Image (sum1 = sum2));
   New_Line;
   Put_Line ("S1 : s2 =");
   Print (sum1);
   New_Line;
   Put_Line ("s2 + s1 =");
   Print (sum2);
   New_Line;
   
   Put_Line ("Identity:");
   Put_Line ("s3 + s3_id equals s3? " & Boolean'Image (sum3 = s3));
   Put_Line ("s3_id + s3_id equals s3_id? " & Boolean'Image (sum4 = s3_id));
   New_Line;
   Put_Line ("s3 + s3_id =");
   Print (sum3);
   Put_Line ("s3_id + s3_id =");
   Print (sum4);

end Main;

Output:

Avalanche:
 4 3 3
 3 1 2
 0 2 3

stable? FALSE

 0 4 3
 4 1 2
 0 2 3

stable? FALSE

 1 0 4
 4 2 2
 0 2 3

stable? FALSE

 1 1 0
 4 2 3
 0 2 3

stable? FALSE

 2 1 0
 0 3 3
 1 2 3

stable? TRUE

Commutivity:
s1 + s2 equals s2 + s2? TRUE

S1 : s2 =
 3 3 3
 3 1 2
 0 2 3


s2 + s1 =
 3 3 3
 3 1 2
 0 2 3


Identity:
s3 + s3_id equals s3? TRUE
s3_id + s3_id equals s3_id? TRUE

s3 + s3_id =
 3 3 3
 3 3 3
 3 3 3

s3_id + s3_id =
 2 1 2
 1 0 1
 2 1 2

ALGOL 68

BEGIN # model Abelian sandpiles #
    # represents a sandpile #
    INT elements = 3;
    MODE SANDPILE = [ 1 : elements, 1 : elements ]INT;
    # returns TRUE if the sandpiles a and b have the same values, FALSE otherwise #
    OP = = ( SANDPILE a, b )BOOL:
    BEGIN
        BOOL result := TRUE;
        FOR i TO elements WHILE result DO
            FOR j TO elements WHILE ( result := a[ i, j ] = b[ i, j ] ) DO SKIP OD
        OD;
        result
    END # = # ;
    # returns TRUE if the sandpile s is stable, FALSE otherwise #
    OP STABLE = ( SANDPILE s )BOOL:
       BEGIN
           BOOL result := TRUE;
           FOR i TO elements WHILE result DO
               FOR j TO elements WHILE result := s[ i, j ] < 4 DO SKIP OD
           OD;
           result
       END # STABLE # ;
    # returns the sandpile s after avalanches #
    OP AVALANCHE = ( SANDPILE s )SANDPILE:
       BEGIN
           SANDPILE result := s;
           WHILE BOOL had avalanche := FALSE;
                 FOR i TO elements DO
                     FOR j TO elements 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 < elements THEN result[ i + 1, j     ] +:= 1 FI;
                             IF j > 1        THEN result[ i,     j - 1 ] +:= 1 FI;
                             IF j < elements THEN result[ i,     j + 1 ] +:= 1 FI
                         FI
                     OD
                 OD;
                 had avalanche
           DO SKIP OD;
           result
       END # AVALANCHE # ;
    # returns the result of adding the sandpile b to a, handling avalanches #
    OP + = ( SANDPILE a, b )SANDPILE:
    BEGIN
        SANDPILE result;
        FOR i TO elements DO
            FOR j TO elements DO result[ i, j ] := a[ i, j ] + b[ i, j ] OD
        OD;
        # handle avalanches #
        AVALANCHE result
    END # + # ;
    # prints the sandpile s #
    PROC show sandpile = ( STRING title, SANDPILE s )VOID:
    BEGIN
        print( ( title, newline ) );
        FOR i TO elements DO
            FOR j TO elements DO
                print( ( " ", whole( s[ i, j ], 0 ) ) )
            OD;
            print( ( newline ) )
        OD
    END # show sandpile # ;
    # task test cases #
    SANDPILE us    = ( ( 4, 3, 3 )
                     , ( 3, 1, 2 )
                     , ( 0, 2, 3 )
                     );
    SANDPILE s1    = ( ( 1, 2, 0 )
                     , ( 2, 1, 1 )
                     , ( 0, 1, 3 )
                     );
    SANDPILE s2    = ( ( 2, 1, 3 )
                     , ( 1, 0, 1 )
                     , ( 0, 1, 0 )
                     );
    SANDPILE s3    = ( ( 3, 3, 3 )
                     , ( 3, 3, 3 )
                     , ( 3, 3, 3 )
                     );
    SANDPILE s3_id = ( ( 2, 1, 2 )
                     , ( 1, 0, 1 )
                     , ( 2, 1, 2 )
                     );
    SANDPILE t := us;
    WHILE NOT STABLE t DO
        show sandpile( "unstable:", t );
        t := AVALANCHE t
    OD;
    show sandpile( "stable: ", t );
    print( ( newline ) );
    show sandpile( "s1:", s1 );
    show sandpile( "s2:", s2 );
    show sandpile( "s1 + s2:", s1 + s2 );
    show sandpile( "s2 + s1:", s2 + s1 );
    print( ( newline ) );
    show sandpile( "s3:", s3 );
    show sandpile( "s3_id:", s3_id );
    print( ( newline ) );
    print( ( "s3 + s3_id = s3 is ", IF s3 + s3_id = s3 THEN "TRUE" ELSE "FALSE" FI, newline ) );
    show sandpile( "s3 + s3_id", s3 + s3_id );
    print( ( "s3_id + s3 = s3 is ", IF s3 + s3_id = s3 THEN "TRUE" ELSE "FALSE" FI, newline ) );
    show sandpile( "s3_id + s3", s3_id + s3 );
    show sandpile( "s3_id + s3_id:", s3_id + s3_id )

END

Output:

unstable:
 4 3 3
 3 1 2
 0 2 3
stable:
 2 1 0
 0 3 3
 1 2 3

s1:
 1 2 0
 2 1 1
 0 1 3
s2:
 2 1 3
 1 0 1
 0 1 0
s1 + s2:
 3 3 3
 3 1 2
 0 2 3
s2 + s1:
 3 3 3
 3 1 2
 0 2 3

s3:
 3 3 3
 3 3 3
 3 3 3
s3_id:
 2 1 2
 1 0 1
 2 1 2

s3 + s3_id = s3 is TRUE
s3 + s3_id
 3 3 3
 3 3 3
 3 3 3
s3_id + s3 = s3 is TRUE
s3_id + s3
 3 3 3
 3 3 3
 3 3 3
s3_id + s3_id:
 2 1 2
 1 0 1
 2 1 2

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <initializer_list>
#include <iostream>

constexpr size_t sp_rows = 3;
constexpr size_t sp_columns = 3;
constexpr size_t sp_cells = sp_rows * sp_columns;
constexpr int sp_limit = 4;

class abelian_sandpile {
    friend std::ostream& operator<<(std::ostream&, const abelian_sandpile&);

public:
    abelian_sandpile();
    explicit abelian_sandpile(std::initializer_list<int> init);
    void stabilize();
    bool is_stable() const;
    void topple();
    abelian_sandpile& operator+=(const abelian_sandpile&);
    bool operator==(const abelian_sandpile&);

private:
    int& cell_value(size_t row, size_t column) {
        return cells_[cell_index(row, column)];
    }
    static size_t cell_index(size_t row, size_t column) {
        return row * sp_columns + column;
    }
    static size_t row_index(size_t cell_index) {
        return cell_index/sp_columns;
    }
    static size_t column_index(size_t cell_index) {
        return cell_index % sp_columns;
    }

    std::array<int, sp_cells> cells_;
};

abelian_sandpile::abelian_sandpile() {
    cells_.fill(0);
}

abelian_sandpile::abelian_sandpile(std::initializer_list<int> init) {
    assert(init.size() == sp_cells);
    std::copy(init.begin(), init.end(), cells_.begin());
}

abelian_sandpile& abelian_sandpile::operator+=(const abelian_sandpile& other) {
    for (size_t i = 0; i < sp_cells; ++i)
        cells_[i] += other.cells_[i];
    stabilize();
    return *this;
}

bool abelian_sandpile::operator==(const abelian_sandpile& other) {
    return cells_ == other.cells_;
}

bool abelian_sandpile::is_stable() const {
    return std::none_of(cells_.begin(), cells_.end(),
                        [](int a) { return a >= sp_limit; });
}

void abelian_sandpile::topple() {
    for (size_t i = 0; i < sp_cells; ++i) {
        if (cells_[i] >= sp_limit) {
            cells_[i] -= sp_limit;
            size_t row = row_index(i);
            size_t column = column_index(i);
            if (row > 0)
                ++cell_value(row - 1, column);
            if (row + 1 < sp_rows)
                ++cell_value(row + 1, column);
            if (column > 0)
                ++cell_value(row, column - 1);
            if (column + 1 < sp_columns)
                ++cell_value(row, column + 1);
            break;
        }
    }
}

void abelian_sandpile::stabilize() {
    while (!is_stable())
        topple();
}

abelian_sandpile operator+(const abelian_sandpile& a, const abelian_sandpile& b) {
    abelian_sandpile c(a);
    c += b;
    return c;
}

std::ostream& operator<<(std::ostream& out, const abelian_sandpile& as) {
    for (size_t i = 0; i < sp_cells; ++i) {
        if (i > 0)
            out << (as.column_index(i) == 0 ? '\n' : ' ');
        out << as.cells_[i];
    }
    return out << '\n';
}

int main() {
    std::cout << std::boolalpha;

    std::cout << "Avalanche:\n";
    abelian_sandpile sp{4,3,3, 3,1,2, 0,2,3};
    while (!sp.is_stable()) {
        std::cout << sp << "stable? " << sp.is_stable() << "\n\n";
        sp.topple();
    }
    std::cout << sp << "stable? " << sp.is_stable() << "\n\n";

    std::cout << "Commutativity:\n";
    abelian_sandpile s1{1,2,0, 2,1,1, 0,1,3};
    abelian_sandpile s2{2,1,3, 1,0,1, 0,1,0};
    abelian_sandpile sum1(s1 + s2);
    abelian_sandpile sum2(s2 + s1);
    std::cout << "s1 + s2 equals s2 + s1? " << (sum1 == sum2) << "\n\n";
    std::cout << "s1 + s2 = \n" << sum1;
    std::cout << "\ns2 + s1 = \n" << sum2;
    std::cout << '\n';

    std::cout << "Identity:\n";
    abelian_sandpile s3{3,3,3, 3,3,3, 3,3,3};
    abelian_sandpile s3_id{2,1,2, 1,0,1, 2,1,2};
    abelian_sandpile sum3(s3 + s3_id);
    abelian_sandpile sum4(s3_id + s3_id);
    std::cout << "s3 + s3_id equals s3? " << (sum3 == s3) << '\n';
    std::cout << "s3_id + s3_id equals s3_id? " << (sum4 == s3_id) << "\n\n";
    std::cout << "s3 + s3_id = \n" << sum3;
    std::cout << "\ns3_id + s3_id = \n" << sum4;

    return 0;
}

Output:

Avalanche:
4 3 3
3 1 2
0 2 3
stable? false

0 4 3
4 1 2
0 2 3
stable? false

1 0 4
4 2 2
0 2 3
stable? false

1 1 0
4 2 3
0 2 3
stable? false

2 1 0
0 3 3
1 2 3
stable? true

Commutativity:
s1 + s2 equals s2 + s1? true

s1 + s2 = 
3 3 3
3 1 2
0 2 3

s2 + s1 = 
3 3 3
3 1 2
0 2 3

Identity:
s3 + s3_id equals s3? true
s3_id + s3_id equals s3_id? true

s3 + s3_id = 
3 3 3
3 3 3
3 3 3

s3_id + s3_id = 
2 1 2
1 0 1
2 1 2

F#

This task uses Abelian Sandpile Model (F#)

let s1=Sandpile(3,3,[|1;2;0;2;1;1;0;1;3|])
let s2=Sandpile(3,3,[|2;1;3;1;0;1;0;1;0|])
printfn "%s\n" ((s1+s2).toS)
printfn "%s\n" ((s2+s1).toS);;
printfn "%s\n" ((s1+s1).toS)
printfn "%s\n" ((s2+s2).toS);;
printfn "%s\n" (Sandpile(3,3,[|4;3;3;3;1;2;0;2;3|])).toS;;
let s3=Sandpile(3,3,(Array.create 9 3))
let s3_id=Sandpile(3,3,[|2;1;2;1;0;1;2;1;2|])
printfn "%s\n" (s3+s3_id).toS
printfn "%s\n" (s3_id+s3_id).toS
//Add together 2 5x5 Sandpiles
let e1=Array.zeroCreate<int> 25 in e1.[12]<-6
let e2=Array.zeroCreate<int> 25 in e2.[12]<-16
printfn "%s\n" ((Sandpile(5,5,e1)+Sandpile(5,5,e2)).toS)

Output:

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

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

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

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

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

[[3; 3; 3]
 [3; 3; 3]
 [3; 3; 3]]

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

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

Factor

I wouldn't call it a translation, but the idea of storing sandpiles as flat arrays came from the Wren entry.

Works with: Factor version 0.99 2020-07-03

USING: arrays grouping io kernel math math.vectors prettyprint
qw sequences ;

CONSTANT: neighbors {
    { 1 3 } { 0 2 4 } { 1 5 } { 0 4 6 } { 1 3 5 7 }
    { 2 4 8 } { 3 7 } { 4 6 8 } { 5 7 }
}

! Sandpile words
: find-tall ( seq -- n ) [ 3 > ] find drop ;
: tall? ( seq -- ? ) find-tall >boolean ;
: distribute ( ind seq -- ) [ [ 1 + ] change-nth ] curry each ;
: adjacent ( n seq -- ) [ neighbors nth ] dip distribute ;
: shrink ( n seq -- ) [ 4 - ] change-nth ;
: (topple) ( n seq -- ) [ shrink ] [ adjacent ] 2bi ;
: topple ( seq -- seq' ) [ find-tall ] [ (topple) ] [ ] tri ;
: avalanche ( seq -- ) [ dup tall? ] [ topple ] while drop ;
: s+ ( seq1 seq2 -- seq3 ) v+ dup avalanche ;

! Output words
: mappend ( seq1 seq2 -- seq3 ) [ flip ] bi@ append flip ;
: sym ( seq str -- seq ) 1array " " 1array tuck 3array mappend ;
: arrow ( seq -- new-seq ) ">" sym ;
: plus  ( seq -- new-seq ) "+" sym ;
: eq    ( seq -- new-seq ) "=" sym ;
: topple> ( seq seq -- seq seq ) arrow over topple 3 group mappend ;
: (.s+) ( seq seq seq -- seq ) [ plus ] [ mappend eq ] [ mappend ] tri* ;
: .s+ ( seq1 seq2 -- ) 2dup s+ [ 3 group ] tri@ (.s+) simple-table. ;

! Task
CONSTANT: s1 { 1 2 0 2 1 1 0 1 3 }
CONSTANT: s2 { 2 1 3 1 0 1 0 1 0 }
CONSTANT: s3 { 3 3 3 3 3 3 3 3 3 }
CONSTANT: id { 2 1 2 1 0 1 2 1 2 }

"Avalanche:" print nl
{ 4 3 3 3 1 2 0 2 3 }
dup 3 group topple> topple> topple> topple> nip simple-table. nl

"s1 + s2 = s2 + s1" print nl
s1 s2 .s+ nl s2 s1 .s+ nl

"s3 + s3_id = s3" print nl
s3 id .s+ nl

"s3_id + s3_id = s3_id" print nl
id id .s+

Output:

Avalanche:

4 3 3   0 4 3   1 0 4   1 1 0   2 1 0
3 1 2 > 4 1 2 > 4 2 2 > 4 2 3 > 0 3 3
0 2 3   0 2 3   0 2 3   0 2 3   1 2 3

s1 + s2 = s2 + s1

1 2 0   2 1 3   3 3 3
2 1 1 + 1 0 1 = 3 1 2
0 1 3   0 1 0   0 2 3

2 1 3   1 2 0   3 3 3
1 0 1 + 2 1 1 = 3 1 2
0 1 0   0 1 3   0 2 3

s3 + s3_id = s3

3 3 3   2 1 2   3 3 3
3 3 3 + 1 0 1 = 3 3 3
3 3 3   2 1 2   3 3 3

s3_id + s3_id = s3_id

2 1 2   2 1 2   2 1 2
1 0 1 + 1 0 1 = 1 0 1
2 1 2   2 1 2   2 1 2

FutureBasic

void local fn SandpilePrint( s(2,2) as ^long )
  long r, c
  for r = 0 to 2
    for c = 0 to 2
      printf @"%ld\t",s(r,c)
    next
    print
  next
  print
end fn

void local fn SandpileTopple( s(2,2) as ^long )
  BOOL stable = NO
  long r, c, value
  
  while ( stable == NO )
    stable = YES
    for r = 0 to 2
      for c = 0 to 2
        value = s(r,c)
        if ( value > 3 )
          s(r,c) -= 4
          if ( r > 0 ) then s(r-1,c)++
          if ( r < 2 ) then s(r+1,c)++
          if ( c > 0 ) then s(r,c-1)++
          if ( c < 2 ) then s(r,c+1)++
          print @"⇣ ⇣ ⇣ ⇣ ⇣"
          print
          fn SandpilePrint( s(0,0) )
          stable = NO : break
        end if
      next
      if ( stable == NO ) then break
    next
  wend
end fn

void local fn SandpileLoad( s(2,2) as ^long, values as CFStringRef )
  long r, c, i = 0
  for r = 0 to 2
    for c = 0 to 2
      s(r,c) = intval(mid(values,i,1))
      i++
    next
  next
end fn

void local fn DoIt
  long r, c, s(2,2), s1(2,2), s2(2,2), s3(2,2), s3_id(2,2)
  
  // s
  text @"Menlo-Bold" : print @"avalanche"
  text @"Menlo" : print @"----------"
  fn SandpileLoad( s(0,0), @"433312023" )
  fn SandpilePrint( s(0,0) )
  fn SandpileTopple( s(0,0) )
  
  // s1
  fn SandpileLoad( s1(0,0), @"120211013" )
  
  // s2
  fn SandpileLoad( s2(0,0), @"213101010" )
  
  // s1 + s2
  for r = 0 to 2
    for c = 0 to 2
      s(r,c) = s1(r,c) + s2(r,c)
    next
  next
  text @"Menlo-Bold" : print @"s1 + s2"
  text @"Menlo" : print @"----------"
  fn SandpileTopple( s(0,0) )
  fn SandpilePrint( s(0,0) )
  
  // s2 + s1
  for r = 0 to 2
    for c = 0 to 2
      s(r,c) = s2(r,c) + s1(r,c)
    next
  next
  text @"Menlo-Bold" : print @"s2 + s1"
  text @"Menlo" : print @"----------"
  fn SandpileTopple( s(0,0) )
  fn SandpilePrint( s(0,0) )
  
  // s3
  fn SandpileLoad( s3(0,0), @"333333333" )
  text @"Menlo-Bold" : print @"s3"
  text @"Menlo" : print @"----------"
  fn SandpilePrint( s3(0,0) )
  
  // s3_id
  fn SandpileLoad( s3_id(0,0), @"212101212" )
  text @"Menlo-Bold" : print @"s3_id"
  text @"Menlo" : print @"----------"
  fn SandpilePrint( s3_id(0,0) )
  
  // s3 + s3_id
  for r = 0 to 2
    for c = 0 to 2
      s(r,c) = s3(r,c) + s3_id(r,c)
    next
  next
  text @"Menlo-Bold" : print @"s3+s3_id"
  text @"Menlo" : print @"----------"
  fn SandpilePrint( s(0,0) )
  fn SandpileTopple( s(0,0) )
  
  // s3_id + s3_id
  for r = 0 to 2
    for c = 0 to 2
      s(r,c) = s3_id(r,c) + s3_id(r,c)
    next
  next
  text @"Menlo-Bold" : print @"s3_id+s3_id"
  text @"Menlo" : print @"-----------"
  fn SandpilePrint( s(0,0) )
  fn SandpileTopple( s(0,0) )
  
end fn

fn DoIt

HandleEvents

Output:

avalanche
----------
4   3   3   
3   1   2   
0   2   3   

⇣ ⇣ ⇣ ⇣ ⇣

0   4   3   
4   1   2   
0   2   3   

⇣ ⇣ ⇣ ⇣ ⇣

1   0   4   
4   2   2   
0   2   3   

⇣ ⇣ ⇣ ⇣ ⇣

1   1   0   
4   2   3   
0   2   3   

⇣ ⇣ ⇣ ⇣ ⇣

2   1   0   
0   3   3   
1   2   3   

s1 + s2
----------
3   3   3   
3   1   2   
0   2   3   

s2 + s1
----------
3   3   3   
3   1   2   
0   2   3   

s3
----------
3   3   3   
3   3   3   
3   3   3   

s3_id
----------
2   1   2   
1   0   1   
2   1   2   

s3+s3_id
----------
5   4   5   
4   3   4   
5   4   5   

⇣ ⇣ ⇣ ⇣ ⇣

1   5   5   
5   3   4   
5   4   5   

⇣ ⇣ ⇣ ⇣ ⇣

2   1   6   
5   4   4   
5   4   5   

⇣ ⇣ ⇣ ⇣ ⇣

2   2   2   
5   4   5   
5   4   5   

⇣ ⇣ ⇣ ⇣ ⇣

3   2   2   
1   5   5   
6   4   5   

⇣ ⇣ ⇣ ⇣ ⇣

3   3   2   
2   1   6   
6   5   5   

⇣ ⇣ ⇣ ⇣ ⇣

3   3   3   
2   2   2   
6   5   6   

⇣ ⇣ ⇣ ⇣ ⇣

3   3   3   
3   2   2   
2   6   6   

⇣ ⇣ ⇣ ⇣ ⇣

3   3   3   
3   3   2   
3   2   7   

⇣ ⇣ ⇣ ⇣ ⇣

3   3   3   
3   3   3   
3   3   3   

s3_id+s3_id
----------
4   2   4   
2   0   2   
4   2   4   

⇣ ⇣ ⇣ ⇣ ⇣

0   3   4   
3   0   2   
4   2   4   

⇣ ⇣ ⇣ ⇣ ⇣

0   4   0   
3   0   3   
4   2   4   

⇣ ⇣ ⇣ ⇣ ⇣

1   0   1   
3   1   3   
4   2   4   

⇣ ⇣ ⇣ ⇣ ⇣

1   0   1   
4   1   3   
0   3   4   

⇣ ⇣ ⇣ ⇣ ⇣

2   0   1   
0   2   3   
1   3   4   

⇣ ⇣ ⇣ ⇣ ⇣

2   0   1   
0   2   4   
1   4   0   

⇣ ⇣ ⇣ ⇣ ⇣

2   0   2   
0   3   0   
1   4   1   

⇣ ⇣ ⇣ ⇣ ⇣

2   0   2   
0   4   0   
2   0   2   

⇣ ⇣ ⇣ ⇣ ⇣

2   1   2   
1   0   1   
2   1   2

Go

Translation of: Wren

package main

import (
    "fmt"
    "strconv"
    "strings"
)

type sandpile struct{ a [9]int }

var neighbors = [][]int{
    {1, 3}, {0, 2, 4}, {1, 5}, {0, 4, 6}, {1, 3, 5, 7}, {2, 4, 8}, {3, 7}, {4, 6, 8}, {5, 7},
}

// 'a' is in row order
func newSandpile(a [9]int) *sandpile { return &sandpile{a} }

func (s *sandpile) plus(other *sandpile) *sandpile {
    b := [9]int{}
    for i := 0; i < 9; i++ {
        b[i] = s.a[i] + other.a[i]
    }
    return &sandpile{b}
}

func (s *sandpile) isStable() bool {
    for _, e := range s.a {
        if e > 3 {
            return false
        }
    }
    return true
}

// just topples once so we can observe intermediate results
func (s *sandpile) topple() {
    for i := 0; i < 9; i++ {
        if s.a[i] > 3 {
            s.a[i] -= 4
            for _, j := range neighbors[i] {
                s.a[j]++
            }
            return
        }
    }
}

func (s *sandpile) String() string {
    var sb strings.Builder
    for i := 0; i < 3; i++ {
        for j := 0; j < 3; j++ {
            sb.WriteString(strconv.Itoa(s.a[3*i+j]) + " ")
        }
        sb.WriteString("\n")
    }
    return sb.String()
}

func main() {
    fmt.Println("Avalanche of topplings:\n")
    s4 := newSandpile([9]int{4, 3, 3, 3, 1, 2, 0, 2, 3})
    fmt.Println(s4)
    for !s4.isStable() {
        s4.topple()
        fmt.Println(s4)
    }

    fmt.Println("Commutative additions:\n")
    s1 := newSandpile([9]int{1, 2, 0, 2, 1, 1, 0, 1, 3})
    s2 := newSandpile([9]int{2, 1, 3, 1, 0, 1, 0, 1, 0})
    s3_a := s1.plus(s2)
    for !s3_a.isStable() {
        s3_a.topple()
    }
    s3_b := s2.plus(s1)
    for !s3_b.isStable() {
        s3_b.topple()
    }
    fmt.Printf("%s\nplus\n\n%s\nequals\n\n%s\n", s1, s2, s3_a)
    fmt.Printf("and\n\n%s\nplus\n\n%s\nalso equals\n\n%s\n", s2, s1, s3_b)

    fmt.Println("Addition of identity sandpile:\n")
    s3 := newSandpile([9]int{3, 3, 3, 3, 3, 3, 3, 3, 3})
    s3_id := newSandpile([9]int{2, 1, 2, 1, 0, 1, 2, 1, 2})
    s4 = s3.plus(s3_id)
    for !s4.isStable() {
        s4.topple()
    }
    fmt.Printf("%s\nplus\n\n%s\nequals\n\n%s\n", s3, s3_id, s4)

    fmt.Println("Addition of identities:\n")
    s5 := s3_id.plus(s3_id)
    for !s5.isStable() {
        s5.topple()
    }
    fmt.Printf("%s\nplus\n\n%s\nequals\n\n%s", s3_id, s3_id, s5)
}

Output:

Avalanche of topplings:

4 3 3 
3 1 2 
0 2 3 

0 4 3 
4 1 2 
0 2 3 

1 0 4 
4 2 2 
0 2 3 

1 1 0 
4 2 3 
0 2 3 

2 1 0 
0 3 3 
1 2 3 

Commutative additions:

1 2 0 
2 1 1 
0 1 3 

plus

2 1 3 
1 0 1 
0 1 0 

equals

3 3 3 
3 1 2 
0 2 3 

and

2 1 3 
1 0 1 
0 1 0 

plus

1 2 0 
2 1 1 
0 1 3 

also equals

3 3 3 
3 1 2 
0 2 3 

Addition of identity sandpile:

3 3 3 
3 3 3 
3 3 3 

plus

2 1 2 
1 0 1 
2 1 2 

equals

3 3 3 
3 3 3 
3 3 3 

Addition of identities:

2 1 2 
1 0 1 
2 1 2 

plus

2 1 2 
1 0 1 
2 1 2 

equals

2 1 2 
1 0 1 
2 1 2 

Haskell

{-# LANGUAGE TupleSections #-}

import Data.List (findIndex, transpose)
import Data.List.Split (chunksOf)

--------------------------- TEST ---------------------------
main :: IO ()
main = do
  let s0 = [[4, 3, 3], [3, 1, 2], [0, 2, 3]]
      s1 = [[1, 2, 0], [2, 1, 1], [0, 1, 3]]
      s2 = [[2, 1, 3], [1, 0, 1], [0, 1, 0]]
      s3_id = [[2, 1, 2], [1, 0, 1], [2, 1, 2]]
      s3 = replicate 3 (replicate 3 3)
      x:xs = reverse $ cascade s0
  mapM_
    putStrLn
    [ "Cascade:"
    , showCascade $ ([], x) : fmap ("->", ) xs
    
    , "s1 + s2 == s2 + s1 -> " <> show (addSand s1 s2 == addSand s2 s1)
    , showCascade [([], s1), (" +", s2), (" =", addSand s1 s2)]
    , showCascade [([], s2), (" +", s1), (" =", addSand s2 s1)]
    
    , "s3 + s3_id == s3 -> " <> show (addSand s3 s3_id == s3)
    , showCascade [([], s3), (" +", s3_id), (" =", addSand s3 s3_id)]
    
    , "s3_id + s3_id == s3_id -> " <> show (addSand s3_id s3_id == s3_id)
    , showCascade [([], s3_id), (" +", s3_id), (" =", addSand s3_id s3_id)]
    ]

------------------------ SAND PILES ------------------------
addSand :: [[Int]] -> [[Int]] -> [[Int]]
addSand xs ys =
  (head . cascade . chunksOf (length xs)) $ zipWith (+) (concat xs) (concat ys)

cascade :: [[Int]] -> [[[Int]]]
cascade xs = chunksOf w <$> convergence (==) (iterate (tumble w) (concat xs))
  where
    w = length xs

convergence :: (a -> a -> Bool) -> [a] -> [a]
convergence p = go
  where
    go (x:ys@(y:_))
      | p x y = [x]
      | otherwise = go ys <> [x]

tumble :: Int -> [Int] -> [Int]
tumble w xs = maybe xs go $ findIndex (w <) xs
  where
    go i = zipWith f [0 ..] xs
      where
        neighbours = indexNeighbours w i
        f j x
          | j `elem` neighbours = succ x
          | i == j = x - succ w
          | otherwise = x

indexNeighbours :: Int -> Int -> [Int]
indexNeighbours w = go
  where
    go i =
      concat
        [ [ j
          | j <- [i - w, i + w] 
          , -1 < j 
          , wSqr > j ]
        , [ pred i
          | 0 /= col ]
        , [ succ i
          | pred w /= col ]
        ]
      where
        wSqr = w * w
        col = rem i w

------------------------- DISPLAY --------------------------
showCascade :: [(String, [[Int]])] -> String
showCascade pairs =
  unlines $
  fmap unwords $
  transpose $
  fmap
    (\(pfx, xs) ->
        unwords <$> transpose (centered pfx : transpose (fmap (fmap show) xs)))
    pairs

centered :: String -> [String]
centered s = [pad, s, pad <> replicate r ' ']
  where
    lng = length s
    pad = replicate lng ' '
    (q, r) = quotRem (2 + lng) 2

Cascade:
 4 3 3    0 4 3    1 0 4    1 1 0    2 1 0
 3 1 2 -> 4 1 2 -> 4 2 2 -> 4 2 3 -> 0 3 3
 0 2 3    0 2 3    0 2 3    0 2 3    1 2 3

s1 + s2 == s2 + s1 -> True
 1 2 0    2 1 3    3 3 3
 2 1 1  + 1 0 1  = 3 1 2
 0 1 3    0 1 0    0 2 3

 2 1 3    1 2 0    3 3 3
 1 0 1  + 2 1 1  = 3 1 2
 0 1 0    0 1 3    0 2 3

s3 + s3_id == s3 -> True
 3 3 3    2 1 2    3 3 3
 3 3 3  + 1 0 1  = 3 3 3
 3 3 3    2 1 2    3 3 3

s3_id + s3_id == s3_id -> True
 2 1 2    2 1 2    2 1 2
 1 0 1  + 1 0 1  = 1 0 1
 2 1 2    2 1 2    2 1 2

J

While=:2 :'u^:(0-.@:-:v)^:_'
index_of_maximum=: $ #: (i. >./)@:,

frame=: ({.~ -@:>:@:$)@:({.~ >:@:$) :. ([;.0~ (1,:_2+$))
NEIGHBORS=: _2]\_1  0  0 _1  0  0  0  1  1  0
AVALANCHE =: 1 1 _4 1 1

avalanche=: (AVALANCHE + {)`[`]}~ ([: <"1 NEIGHBORS +"1 index_of_maximum)
erode=: avalanche&.:frame While(3 < [: >./ ,)

   NB. common ways to construct a matrix in j from directly entered vectors
   s3_id=: >2 1 2;1 0 1;2 1 2 NB. 3 3$2 1 2 1 0 1 2 1 2 NB. _3]\2 1 2 1 0 1 2 1 2 NB. 2 1 2,1 0 1,:2 1 2
   s3=: 3 3 $ 3 NB. ($~,~)3 NB. 3"0 i.3 3

   matches =: -:
   Commutes=: adverb def '(u matches u~)~'

   NB. demonstrate Commutes adbverb
   4 - Commutes 3
0
   4 + Commutes 3
1

   NB. confirmation
   <"2 A , ] avalanche&.:frame@:([ 3 :'A=:A,y') While(3 < [: >./ ,) 10#.inv 433 312 023 [ A=:0 3 3$0
┌─────┬─────┬─────┬─────┬─────┐
│4 3 3│0 4 3│1 0 4│1 1 0│2 1 0│
│3 1 2│4 1 2│4 2 2│4 2 3│0 3 3│
│0 2 3│0 2 3│0 2 3│0 2 3│1 2 3│
└─────┴─────┴─────┴─────┴─────┘

   NB. matrix addition commutes
   's1 s2'=: 120 211 013 ;&:(10&#.inv) 213 101 010
   s1 + Commutes s2
1
   erode s1 + s2
3 3 3
3 1 2
0 2 3

   NB. use: IDENTITY verify_identity MATRIX
   verify_identity=: (erode@:+ matches ]) erode

   raku_id verify_identity raku
1

   (; erode) raku
┌─────────┬─────────┐
│4 1 0 5 1│1 3 2 1 0│
│9 3 6 1 0│2 2 3 3 1│
│8 1 2 5 3│1 1 2 0 3│
│3 0 1 7 5│2 0 3 2 0│
│4 2 2 4 0│3 2 3 2 1│
└─────────┴─────────┘

Java

import java.util.ArrayList;
import java.util.List;

public final class AbelianSandpileModel {

	public static void main(String[] aArgs) {
		Sandpile avalanche = new Sandpile(List.of( 4, 3, 3, 3, 1, 2, 0, 2, 3 ));
		System.out.println("Avalanche reduction to stable state:");
		avalanche.display();
		System.out.println(" ==> ");
		avalanche.stabilise();
		avalanche.display();
		
		Sandpile s1 = new Sandpile(List.of( 1, 2, 0, 2, 1, 1, 0, 1, 3 ));
		Sandpile s2 = new Sandpile(List.of( 2, 1, 3, 1, 0, 1, 0, 1, 0 ));
		Sandpile sum1 = s1.add(s2);
		Sandpile sum2 = s2.add(s1);
		System.out.println(System.lineSeparator() + "Commutativity of addition" + System.lineSeparator());
		System.out.println("Sandpile1 + Sandpile2:");
		sum1.display();
		System.out.println("Sandpile2 + Sandpile1:");
		sum2.display();
		System.out.println("Sandpile1 + Sandpile2 = Sandpile2 + Sandpile1: " + sum1.equals(sum2));
		
		Sandpile s3 = new Sandpile(List.of( 3, 3, 3, 3, 3, 3, 3, 3, 3 ));
		Sandpile s3_id = new Sandpile(List.of( 2, 1, 2, 1, 0, 1, 2, 1, 2 ));
		Sandpile sum3 = s3.add(s3_id);
		Sandpile sum4 = s3_id.add(s3_id);
		System.out.println(System.lineSeparator() + "Identity Sandpile" + System.lineSeparator());
		System.out.println("Sandpile3 + Sandpile3_id:");
		sum3.display();
		System.out.println("Sandpile3_id + Sandpile3_id:");
		sum4.display();		
	}

}

final class Sandpile {
	
    public Sandpile(List<Integer> aList) {
    	if ( aList.size() != CELL_COUNT ) {
    		throw new IllegalArgumentException("Initialiser list must contain " + CELL_COUNT + " elements");
    	}    	
    	cells = new ArrayList<Integer>(aList);    	
    }
    
    public void stabilise() {
    	while ( ! isStable() ) {
            topple();
    	}
    }
    
    public boolean isStable() {
    	return cells.stream().noneMatch( i -> i >= CELL_LIMIT );
    }
     
    public void topple() {
    	for ( int i = 0; i < CELL_COUNT; i++ ) {
    		if ( cells.get(i) >= CELL_LIMIT ) {
	            cells.set(i, cells.get(i) - CELL_LIMIT);
	            final int row = rowIndex(i);
	            final int col = colIndex(i);
	            if ( row > 0 ) {
	            	increment(row - 1, col);
	            }
	            if ( row + 1 < ROW_COUNT ) {
	                increment(row + 1, col);
	            }
	            if ( col > 0 ) {
	                increment(row, col - 1);
	            }
	            if ( col + 1 < COL_COUNT ) {
	                increment(row, col + 1);
	            }
    		}
	    }
    }
    
    public Sandpile add(Sandpile aOther) {
    	List<Integer> list = new ArrayList<Integer>();
    	for ( int i = 0; i < CELL_COUNT; i++ ) {
            list.add(cells.get(i) + aOther.cells.get(i));
    	}
    	Sandpile result = new Sandpile(list);
        result.stabilise();
        return result;
    }
    
    public boolean equals(Sandpile aOther) {
    	return cells.equals(aOther.cells);
    } 
    
    public void display() {
    	for ( int i = 0; i < CELL_COUNT; i++ ) {
    		System.out.print(cells.get(i));
    		System.out.print( ( colIndex(i + 1) == 0 ) ? System.lineSeparator() : " ");
    	}
    }
  
	private void increment(int aRow, int aCol) {
		final int index = cellIndex(aRow, aCol);
		cells.set(index, cells.get(index) + 1);
	}
    
    private static int cellIndex(int aRow, int aCol) {
        return aRow * COL_COUNT + aCol;
    }
    
    private static int rowIndex(int aCellIndex) {
        return aCellIndex / COL_COUNT;
    }
    
    private static int colIndex(int aCellIndex) {
        return aCellIndex % COL_COUNT;
    }
	
    private List<Integer> cells;

	private static final int ROW_COUNT = 3;
	private static final int COL_COUNT = 3;
	private static final int CELL_COUNT = ROW_COUNT * COL_COUNT;
	private static final int CELL_LIMIT = 4;
	
}

Output:

Avalanche reduction to stable state:
4 3 3
3 1 2
0 2 3
 ==> 
2 1 0
0 3 3
1 2 3

Commutativity of addition

Sandpile1 + Sandpile2:
3 3 3
3 1 2
0 2 3
Sandpile2 + Sandpile1:
3 3 3
3 1 2
0 2 3
Sandpile1 + Sandpile2 = Sandpile2 + Sandpile1: true

Identity Sandpile

Sandpile3 + Sandpile3_id:
3 3 3
3 3 3
3 3 3
Sandpile3_id + Sandpile3_id:
2 1 2
1 0 1
2 1 2

jq

Adapted from Wren

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

# `whilst/2` is like `while/2` but emits the final term rather than the first one
def whilst(cond; update):
     def _whilst:
         if cond then update | (., _whilst) else empty end;
     _whilst;

# module Sandpile

def new($a): {$a};

def neighbors: [
   [1, 3], [0, 2, 4], [1, 5],
   [0, 4, 6], [1, 3, 5, 7], [2, 4, 8],
   [3, 7], [4, 6, 8], [5, 7]
];

def add($other):
  . as $in
  | reduce range(0; .a|length) as $i ($in; .a[$i] += $other.a[$i] );

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

# 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 tos:
  . as $in
  | reduce range(0;3) as $i ("";
      reduce range(0;3) as $j (.;
        . + " \($in.a[3*$i + $j])" )
      | . +"\n" );

# Some sandpiles:
def s1: new([1, 2, 0, 2, 1, 1, 0, 1, 3]);
def s2: new([2, 1, 3, 1, 0, 1, 0, 1, 0]);
def s3: new([range(0;9)|3]);
def s4: new([4, 3, 3, 3, 1, 2, 0, 2, 3]);

def s3_id: new([2, 1, 2, 1, 0, 1, 2, 1, 2]);

# For brevity
def report_add($s1; $s2):
  "\($s1|tos)\nplus\n\n\($s2|tos)\nequals\n\n\($s1 | add($s2) | until(isStable; topple) | tos)";

def task1:
  "Avalanche of topplings:\n",
  (s4
   | (., whilst(isStable|not; topple))
   | tos ) ;

def task2:
  def s3_a: s1 | add(s2);
  def s3_b: s2 | add(s1);

  "Commutative additions:\n",
  ( (s3_b | until(isStable; topple)) as $s3_b
    | report_add(s1; s2),
      "and\n\n\(s2|tos)\nplus\n\n\(s1|tos)\nalso equals\n\n\($s3_b|tos)" ) ;

def task3:
  "Addition of identity sandpile:\n",
  report_add(s3; s3_id);

def task4:
  "Addition of identities:\n",
  report_add(s3_id; s3_id);

task1, task2, task3, task4

Output:

As for Wren.

Julia

import Base.+, Base.print

struct Sandpile
    pile::Matrix{UInt8}
end

function Sandpile(s::String)
    arr = [parse(UInt8, x.match) for x in eachmatch(r"\d+", s)]
    siz = isqrt(length(arr))
    return Sandpile(reshape(UInt8.(arr), siz, siz)')
end

const HMAX = 3

function avalanche!(s::Sandpile, lim=HMAX)
    nrows, ncols = size(s.pile)
    while any(x -> x > lim, s.pile)
        for j in 1:ncols, i in 1:nrows
            if s.pile[i, j] > lim
                i > 1 && (s.pile[i - 1, j] += 1)
                i < nrows && (s.pile[i + 1, j] += 1)
                j > 1 && (s.pile[i, j - 1] += 1)
                j < ncols && (s.pile[i, j + 1] += 1)
                s.pile[i, j] -= 4
            end
        end
    end
    s
end

+(s1::Sandpile, s2::Sandpile) = avalanche!(Sandpile((s1.pile + s2.pile)))

function print(io::IO, s::Sandpile)
    for row in 1:size(s.pile)[1]
        for col in 1:size(s.pile)[2]
            print(io, lpad(s.pile[row, col], 4))
        end
        println()
    end
end

const s1 = Sandpile("""
    1 2 0
    2 1 1
    0 1 3""")

const s2 = Sandpile("""
    2 1 3
    1 0 1
    0 1 0""")

const s3 = Sandpile("""
    3 3 3
    3 3 3
    3 3 3""")

const s3_id = Sandpile("""
    2 1 2
    1 0 1
    2 1 2""")

const s3a = Sandpile("""
    4 3 3
    3 1 2
    0 2 3""")

println("Avalanche reduction to group:\n", s3a, "   =>")
println(avalanche!(s3a), "\n")

println("Commutative Property:\ns1 + s2 =\n", s1 + s2, "\ns2 + s1 =\n", s2 + s1, "\n")

println("Addition:\n", s3, "   +\n", s3_id, "   =\n", s3 + s3_id, "\n")
println(s3_id, "   +\n", s3_id, "   =\n", s3_id + s3_id, "\n")

Output:

Avalanche reduction to group:
   4   3   3
   3   1   2
   0   2   3
   =>
   2   1   0
   0   3   3
   1   2   3


Commutative Property:
s1 + s2 =
   3   3   3
   3   1   2
   0   2   3

s2 + s1 =
   3   3   3
   3   1   2
   0   2   3


Addition:
   3   3   3
   3   3   3
   3   3   3
   +
   2   1   2
   1   0   1
   2   1   2
   =
   3   3   3
   3   3   3
   3   3   3


   2   1   2
   1   0   1
   2   1   2
   +
   2   1   2
   1   0   1
   2   1   2
   =
   2   1   2
   1   0   1
   2   1   2

Lua

Uses Abelian sandpile model here, then extends..

sandpile.__index = sandpile
sandpile.new = function(self, vals)
  local inst = setmetatable({},sandpile)
  inst.cell, inst.dim = {}, #vals
  for r = 1, inst.dim do
    inst.cell[r] = {}
    for c = 1, inst.dim do
      inst.cell[r][c] = vals[r][c]
    end
  end
  return inst
end
sandpile.add = function(self, other)
  local vals = {}
  for r = 1, self.dim do
    vals[r] = {}
    for c = 1, self.dim do
      vals[r][c] = self.cell[r][c] + other.cell[r][c]
    end
  end
  local inst = sandpile:new(vals)
  inst:iter()
  return inst
end

local s1 = sandpile:new{{1,2,0},{2,1,1},{0,1,3}}
local s2 = sandpile:new{{2,1,3},{1,0,1},{0,1,0}}
print("s1 =")  s1:draw()
print("\ns2 =")  s2:draw()
local s1ps2 = s1:add(s2)
print("\ns1 + s2 =")  s1ps2:draw()
local s2ps1 = s2:add(s1)
print("\ns2 + s1 =")  s2ps1:draw()
local topple = sandpile:new{{4,3,3},{3,1,2},{0,2,3}}
print("\ntopple, before =")  topple:draw()
topple:iter()
print("\ntopple, after =")  topple:draw()
local s3 = sandpile:new{{3,3,3},{3,3,3},{3,3,3}}
print("\ns3 =")  s3:draw()
local s3_id = sandpile:new{{2,1,2},{1,0,1},{2,1,2}}
print("\ns3_id =") s3_id:draw()
local s3ps3_id = s3:add(s3_id)
print("\ns3 + s3_id =")  s3ps3_id:draw()
local s3_idps3_id = s3_id:add(s3_id)
print("\ns3_id + s3_id =")  s3_idps3_id:draw()

Output:

s1 =
1 2 0
2 1 1
0 1 3

s2 =
2 1 3
1 0 1
0 1 0

s1 + s2 =
3 3 3
3 1 2
0 2 3

s2 + s1 =
3 3 3
3 1 2
0 2 3

topple, before =
4 3 3
3 1 2
0 2 3

topple, after =
2 1 0
0 3 3
1 2 3

s3 =
3 3 3
3 3 3
3 3 3

s3_id =
2 1 2
1 0 1
2 1 2

s3 + s3_id =
3 3 3
3 3 3
3 3 3

s3_id + s3_id =
2 1 2
1 0 1
2 1 2

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]
  ]
Format[x_sp] := Grid[x[[1]]]

s1 = sp[{{1, 2, 0}, {2, 1, 1}, {0, 1, 3}}];
s2 = sp[{{2, 1, 3}, {1, 0, 1}, {0, 1, 0}}];
s3 = sp[ConstantArray[3, {3, 3}]];
s3id = sp[{{2, 1, 2}, {1, 0, 1}, {2, 1, 2}}];

s1 + s2
s2 + s1
sp[{{4, 3, 3}, {3, 1, 2}, {0, 2, 3}}] + sp[0]
s3 + s3id === s3
s3id + s3id === s3id

Output:

3	3	3
3	1	2
0	2	3

3	3	3
3	1	2
0	2	3

2	1	0
0	3	3
1	2	3

True

True

Nim

import sequtils
import strutils

type SandPile = array[3, array[3, int]]

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

iterator neighbors(i, j: int): tuple[a, b: int] =
  ## Yield the indexes of the neighbours of cell at indexes (i, j).
  if i > 0:
    yield (i - 1, j)
  if i < 2:
    yield (i + 1, j)
  if j > 0:
    yield (i, j - 1)
  if j < 2:
    yield (i, j + 1)

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

proc print(s: openArray[SandPile]) =
  ## Print a list of sandpiles.
  for i in 0..2:
    for n, sp in s:
      if n != 0:
        stdout.write(if i == 1: " ⇨ " else: "   ")
      stdout.write(sp[i].join(" "))
    stdout.write('\n')

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

proc printSum(s1, s2, s3: SandPile) =
  ## Print "s1 + s2 = s3".
  for i in 0..2:
    stdout.write(s1[i].join(" "))
    stdout.write(if i == 1: " + " else: "   ", s2[i].join(" "))
    stdout.write(if i == 1: " = " else: "   ", s3[i].join(" "))
    stdout.write('\n')

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

func isStable(sandPile: SandPile): bool =
  ## Return true if the sandpile is stable, else false.
  result = true
  for row in sandPile:
    if row.anyit(it > 3):
      return false

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

proc topple(sandPile: var SandPile) =
  ## Eliminate one value > 3, propagating a grain to each neighbor.
  for i, row in sandPile:
    for j, val in row:
      if val > 3:
        dec sandPile[i][j], 4
        for (i, j) in neighbors(i, j):
          inc sandPile[i][j]
        return

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

proc stabilize(sandPile: var SandPile) =
  ## Stabilize a sandpile.
  while not sandPile.isStable():
    sandPile.topple()

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

proc `+`(s1, s2: SandPile): SandPile =
  ## Add two sandpiles, stabilizing the result.
  for row in 0..2:
    for col in 0..2:
      result[row][col] = s1[row][col] + s2[row][col]
  result.stabilize()

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

const Separator = "\n-----\n"

echo "Avalanche\n"
var s: SandPile = [[4, 3, 3], [3, 1, 2], [0, 2, 3]]
var list = @[s]
while not s.isStable():
  s.topple()
  list.add(s)
list.print()
echo Separator

echo "s1 + s2 == s2 + s1\n"
let s1 = [[1, 2, 0], [2, 1, 1], [0, 1, 3]]
let s2 = [[2, 1, 3], [1, 0, 1], [0, 1, 0]]
printSum(s1, s2, s1 + s2)
echo ""
printSum(s2, s1, s2 + s1)
echo Separator

echo "s3 + s3_id == s3\n"
let s3 = [[3, 3, 3], [3, 3, 3], [3, 3, 3]]
let s3_id = [[2, 1, 2], [1, 0, 1], [2, 1, 2]]
printSum(s3, s3_id, s3 + s3_id)
echo Separator

echo "s3_id + s3_id = s3_id\n"
printSum(s3_id, s3_id, s3_id + s3_id)

Output:

Avalanche

4 3 3   0 4 3   1 0 4   1 1 0   2 1 0
3 1 2 ⇨ 4 1 2 ⇨ 4 2 2 ⇨ 4 2 3 ⇨ 0 3 3
0 2 3   0 2 3   0 2 3   0 2 3   1 2 3

-----

s1 + s2 == s2 + s1

1 2 0   2 1 3   3 3 3
2 1 1 + 1 0 1 = 3 1 2
0 1 3   0 1 0   0 2 3

2 1 3   1 2 0   3 3 3
1 0 1 + 2 1 1 = 3 1 2
0 1 0   0 1 3   0 2 3

-----

s3 + s3_id == s3

3 3 3   2 1 2   3 3 3
3 3 3 + 1 0 1 = 3 3 3
3 3 3   2 1 2   3 3 3

-----

s3_id + s3_id = s3_id

2 1 2   2 1 2   2 1 2
1 0 1 + 1 0 1 = 1 0 1
2 1 2   2 1 2   2 1 2

OCaml

(* https://en.wikipedia.org/wiki/Abelian_sandpile_model *)

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

    type t = { grid : int array array ; unstable : ((int*int),unit) Hashtbl.t }
           
    let make () = { grid = Array.init M.m (fun _ -> Array.make M.n 0); unstable = Hashtbl.create 10 }     

    let print {grid=grid} =
      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 add_grain {grid=grid;unstable=unstable} 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 ({grid=grid;unstable=unstable} as s) x y
      = grid.(x).(y) <- grid.(x).(y) - 4
      ; if grid.(x).(y) < 4
        then Hashtbl.remove unstable (x,y)
      ; let add_grain = add_grain s in 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 s n x y
      = for i = 1 to n
        do add_grain s x y
        done

    let avalanche ?(avalanche_print=fun _ -> ()) ({grid=grid;unstable=unstable} as s)
      = 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 s x y; avalanche_print s ) unstable'
        done

    let init ?(avalanche_print=fun _ -> ()) f
      = let s = { grid = Array.init M.m (fun x -> Array.init M.n (fun y -> f x y)) ; unstable = Hashtbl.create 10 }
        in Array.iteri (fun x -> Array.iteri (fun y e -> if e >= 4 then Hashtbl.replace s.unstable (x,y) ())) s.grid
         ; avalanche_print s
         ; avalanche ~avalanche_print s
         ; s
    
    let sandpile n
      = let s = make ()
        in add_sand s n (M.m/2) (M.n/2)
         ; avalanche s
         ; s
                                        
    let (+.) {grid=a} {grid=b}
      = let c = init (fun x y -> a.(x).(y) + b.(x).(y))
        in avalanche c
         ; c
  end

(* testing *)
   
let ()
  = let module S = Make (struct let m = 3 let n = 3 end)
    in let open S
       in print_endline "Avalanche example"
        ; begin
           let s0 = init ~avalanche_print:(fun s -> print s
                                                  ; print_endline "  ↓")
                      (fun x y -> [| [| 4 ; 3 ; 3 |] 
                                   ; [| 3 ; 1 ; 2 |]
                                   ; [| 0 ; 2 ; 3 |]
                                  |].(x).(y))
           in print s0
            ; print_endline "---------------"
         end
        ; print_endline "Addition example"
        ; begin 
            let  s1 = init (fun x y -> [| [| 1 ; 2 ; 0 |]
                                        ; [| 2 ; 1 ; 1 |]
                                        ; [| 0 ; 1 ; 3 |]
                                       |].(x).(y))
             and s2 = init (fun x y -> [| [| 2 ; 1 ; 3 |]
                                        ; [| 1 ; 0 ; 1 |]
                                        ; [| 0 ; 1 ; 0|]
                                       |].(x).(y))
             and s3 = init (fun _ _ -> 3)
             and s3_id = init (fun x y -> match (x,y) with
                                          | ((0,0)|(2,0)|(0,2)|(2,2)) -> 2
                                          | ((1,0)|(1,2)|(0,1)|(2,1)) -> 1
                                          | _ -> 0)
            in print s1
             ; print_endline "  +"
             ; print s2
             ; print_endline "  ="
             ; print (s1 +. s2)
             ; print_endline "------ Identity examples -----"
             ; print s3
             ; print_endline "  +"
             ; print s3_id
             ; print_endline "  ="
             ; print (s3 +. s3_id)
             ; print_endline "-----"
             ; print s3_id
             ; print_endline "  +"
             ; print s3_id
             ; print_endline "  ="
             ; print (s3_id +. s3_id)
          end

Output:

Avalanche example
4 3 3 
3 1 2 
0 2 3 
  ↓
0 4 3 
4 1 2 
0 2 3 
  ↓
1 4 3 
0 2 2 
1 2 3 
  ↓
2 0 4 
0 3 2 
1 2 3 
  ↓
2 1 0 
0 3 3 
1 2 3 
  ↓
2 1 0 
0 3 3 
1 2 3 
---------------
Addition example
1 2 0 
2 1 1 
0 1 3 
  +
2 1 3 
1 0 1 
0 1 0 
  =
3 3 3 
3 1 2 
0 2 3 
------ Identity examples -----
3 3 3 
3 3 3 
3 3 3 
  +
2 1 2 
1 0 1 
2 1 2 
  =
3 3 3 
3 3 3 
3 3 3 
-----
2 1 2 
1 0 1 
2 1 2 
  +
2 1 2 
1 0 1 
2 1 2 
  =
2 1 2 
1 0 1 
2 1 2 

Phix

constant s1 = {"1 2 0",
               "2 1 1",
               "0 1 3"},
 
         s2 = {"2 1 3",
               "1 0 1",
               "0 1 0"},
 
         s3 = {"3 3 3",
               "3 3 3",
               "3 3 3"},
 
      s3_id = {"2 1 2",
               "1 0 1",
               "2 1 2"},
 
         s4 = {"4 3 3",
               "3 1 2",
               "0 2 3"}
 
function add(sequence s, t)
    for i=1 to 3 do
        for j=1 to 5 by 2 do
            s[i][j] += t[i][j]-'0'
        end for
    end for
    return s
end function
 
function topple(sequence s, integer one=0)
    for i=1 to 3 do
        for j=1 to 5 by 2 do
            if s[i][j]>'3' then
                s[i][j] -= 4
                if i>1 then s[i-1][j] += 1 end if
                if i<3 then s[i+1][j] += 1 end if
                if j>1 then s[i][j-2] += 1 end if
                if j<5 then s[i][j+2] += 1 end if
                if one=1 then return s end if
                one = -1
            end if
        end for
    end for
    return iff(one=1?{}:iff(one=-1?topple(s):s))
end function
 
procedure shout(sequence s)
    sequence r = repeat("",5)
    for i=1 to length(s) do
        sequence si = s[i]
        if string(si) then
            string ti = repeat(' ',length(si))
            r[1] &= ti
            r[2] &= si
            r[3] &= ti
        else
            for j=1 to 3 do
                r[j] &= si[j]
            end for
        end if
    end for
    puts(1,join(r,"\n"))
end procedure
 
puts(1,"1. Show avalanche\n\n")
sequence s = s4,
         res = {"    ",s}
while true do
    s = topple(s,1)
    if s={} then exit end if
    res &= {" ==> ",s}
end while
shout(res)
 
puts(1,"2. Prove s1 + s2 = s2 + s1\n\n")
shout({"    ",s1,"  +  ",s2,"  =  ",topple(add(s1,s2))})
shout({"    ",s2,"  +  ",s1,"  =  ",topple(add(s2,s1))})
 
puts(1,"3. Show that s3 + s3_id == s3\n\n")
shout({"    ",s3,"  +  ",s3_id,"  =  ",topple(add(s3,s3_id))})
 
puts(1,"4. Show that s3_id + s3_id == s3_id\n\n")
shout({"    ",s3_id,"  +  ",s3_id,"  =  ",topple(add(s3_id,s3_id))})

Output:

1. Show avalanche

    4 3 3     0 4 3     1 0 4     1 1 0     2 1 0
    3 1 2 ==> 4 1 2 ==> 4 2 2 ==> 4 2 3 ==> 0 3 3
    0 2 3     0 2 3     0 2 3     0 2 3     1 2 3

2. Prove s1 + s2 = s2 + s1

    1 2 0     2 1 3     3 3 3
    2 1 1  +  1 0 1  =  3 1 2
    0 1 3     0 1 0     0 2 3

    2 1 3     1 2 0     3 3 3
    1 0 1  +  2 1 1  =  3 1 2
    0 1 0     0 1 3     0 2 3

3. Show that s3 + s3_id == s3

    3 3 3     2 1 2     3 3 3
    3 3 3  +  1 0 1  =  3 3 3
    3 3 3     2 1 2     3 3 3

4. Show that s3_id + s3_id == s3_id

    2 1 2     2 1 2     2 1 2
    1 0 1  +  1 0 1  =  1 0 1
    2 1 2     2 1 2     2 1 2

Python

Object Oriented

from itertools import product
from collections import defaultdict


class Sandpile():
    def __init__(self, gridtext):
        array = [int(x) for x in gridtext.strip().split()]
        self.grid = defaultdict(int,
                                {(i //3, i % 3): x 
                                 for i, x in enumerate(array)})

    _border = set((r, c) 
                  for r, c in product(range(-1, 4), repeat=2) 
                  if not 0 <= r <= 2 or not 0 <= c <= 2
                  )
    _cell_coords = list(product(range(3), repeat=2))
    
    def topple(self):
        g = self.grid
        for r, c in self._cell_coords:
            if g[(r, c)] >= 4:
                g[(r - 1, c)] += 1
                g[(r + 1, c)] += 1
                g[(r, c - 1)] += 1
                g[(r, c + 1)] += 1
                g[(r, c)] -= 4
                return True
        return False
    
    def stabilise(self):
        while self.topple():
            pass
        # Remove extraneous grid border
        g = self.grid
        for row_col in self._border.intersection(g.keys()):
            del g[row_col]
        return self
    
    __pos__ = stabilise     # +s == s.stabilise()
    
    def __eq__(self, other):
        g = self.grid
        return all(g[row_col] == other.grid[row_col]
                   for row_col in self._cell_coords)

    def __add__(self, other):
        g = self.grid
        ans = Sandpile("")
        for row_col in self._cell_coords:
            ans.grid[row_col] = g[row_col] + other.grid[row_col]
        return ans.stabilise()
       
    def __str__(self):
        g, txt = self.grid, []
        for row in range(3):
            txt.append(' '.join(str(g[(row, col)]) 
                                for col in range(3)))
        return '\n'.join(txt)
    
    def __repr__(self):
        return f'{self.__class__.__name__}(""""\n{self.__str__()}""")'
        

unstable = Sandpile("""
4 3 3
3 1 2
0 2 3""")
s1 = Sandpile("""
    1 2 0
    2 1 1
    0 1 3
""")
s2 = Sandpile("""
    2 1 3
    1 0 1
    0 1 0 
""")
s3 = Sandpile("3 3 3  3 3 3  3 3 3")
s3_id = Sandpile("2 1 2  1 0 1  2 1 2")

Command line session to complete task.

In [2]: unstable
Out[2]: 
Sandpile(""""
4 3 3
3 1 2
0 2 3""")

In [3]: unstable.stabilise()
Out[3]: 
Sandpile(""""
2 1 0
0 3 3
1 2 3""")

In [4]: s1 + s2
Out[4]: 
Sandpile(""""
3 3 3
3 1 2
0 2 3""")

In [5]: s2 + s1
Out[5]: 
Sandpile(""""
3 3 3
3 1 2
0 2 3""")

In [6]: s1 + s2 == s2 + s1
Out[6]: True

In [7]: s3
Out[7]: 
Sandpile(""""
3 3 3
3 3 3
3 3 3""")

In [8]: s3_id
Out[8]: 
Sandpile(""""
2 1 2
1 0 1
2 1 2""")

In [9]: s3 + s3_id
Out[9]: 
Sandpile(""""
3 3 3
3 3 3
3 3 3""")

In [10]: s3 + s3_id == s3
Out[10]: True

In [11]: s3_id + s3_id
Out[11]: 
Sandpile(""""
2 1 2
1 0 1
2 1 2""")

In [12]: s3_id + s3_id == s3_id
Out[12]: True

In [13]: 

Functional

'''Abelian Sandpile – Identity'''

from operator import add, eq


# -------------------------- TEST --------------------------
# main :: IO ()
def main():
    '''Tests of cascades and additions'''
    s0 = [[4, 3, 3], [3, 1, 2], [0, 2, 3]]
    s1 = [[1, 2, 0], [2, 1, 1], [0, 1, 3]]
    s2 = [[2, 1, 3], [1, 0, 1], [0, 1, 0]]
    s3 = [[3, 3, 3], [3, 3, 3], [3, 3, 3]]
    s3_id = [[2, 1, 2], [1, 0, 1], [2, 1, 2]]

    series = list(cascadeSeries(s0))
    for expr in [
            'Cascade:',
            showSandPiles(
                [(' ', series[0])] + [
                    (':', xs) for xs in series[1:]
                ]
            ),
            '',
            f's1 + s2 == s2 + s1 -> {addSand(s1)(s2) == addSand(s2)(s1)}',
            showSandPiles([
                (' ', s1),
                ('+', s2),
                ('=', addSand(s1)(s2))
            ]),
            '',
            showSandPiles([
                (' ', s2),
                ('+', s1),
                ('=', addSand(s2)(s1))
            ]),
            '',
            f's3 + s3_id == s3 -> {addSand(s3)(s3_id) == s3}',
            showSandPiles([
                (' ', s3),
                ('+', s3_id),
                ('=', addSand(s3)(s3_id))
            ]),
            '',
            f's3_id + s3_id == s3_id -> {addSand(s3_id)(s3_id) == s3_id}',
            showSandPiles([
                (' ', s3_id),
                ('+', s3_id),
                ('=', addSand(s3_id)(s3_id))
            ]),

    ]:
        print(expr)


# ----------------------- SANDPILES ------------------------

# addSand :: [[Int]] -> [[Int]] -> [[Int]]
def addSand(xs):
    '''The stabilised sum of two sandpiles.
    '''
    def go(ys):
        return cascadeSeries(
            chunksOf(len(xs))(
                map(
                    add,
                    concat(xs),
                    concat(ys)
                )
            )
        )[-1]
    return go


# cascadeSeries :: [[Int]] -> [[[Int]]]
def cascadeSeries(rows):
    '''The sequence of states from a given
       sand pile to a stable condition.
    '''
    xs = list(rows)
    w = len(xs)
    return [
        list(chunksOf(w)(x)) for x
        in convergence(eq)(
            iterate(nextState(w))(
                concat(xs)
            )
        )
    ]


# convergence :: (a -> a -> Bool) -> [a] -> [a]
def convergence(p):
    '''All items of xs to the point where the binary
       p returns True over two successive values.
    '''
    def go(xs):
        def conv(prev, ys):
            y = next(ys)
            return [prev] + (
                [] if p(prev, y) else conv(y, ys)
            )
        return conv(next(xs), xs)
    return go


# nextState Int -> Int -> [Int] -> [Int]
def nextState(w):
    '''The next state of a (potentially unstable)
       flattened sand-pile matrix of row length w.
    '''
    def go(xs):
        def tumble(i):
            neighbours = indexNeighbours(w)(i)
            return [
                1 + k if j in neighbours else (
                    k - (1 + w) if j == i else k
                ) for (j, k) in enumerate(xs)
            ]
        return maybe(xs)(tumble)(
            findIndex(lambda x: w < x)(xs)
        )
    return go


# indexNeighbours :: Int -> Int -> [Int]
def indexNeighbours(w):
    '''Indices vertically and horizontally adjoining the
       given index in a flattened matrix of dimension w.
    '''
    def go(i):
        lastCol = w - 1
        iSqr = (w * w)
        col = i % w
        return [
            j for j in [i - w, i + w]
            if -1 < j < iSqr
        ] + ([i - 1] if 0 != col else []) + (
            [1 + i] if lastCol != col else []
        )
    return go


# ------------------------ DISPLAY -------------------------

# showSandPiles :: [(String, [[Int]])] -> String
def showSandPiles(pairs):
    '''Indented multi-line representation
       of a sequence of matrices, delimited
       by preceding operators or indents.
    '''
    return '\n'.join([
        ' '.join([' '.join(map(str, seq)) for seq in tpl])
        for tpl in zip(*[
            zip(
                *[list(str(pfx).center(len(rows)))]
                + list(zip(*rows))
            )
            for (pfx, rows) in pairs
        ])
    ])


# ------------------------ GENERIC -------------------------

# chunksOf :: Int -> [a] -> [[a]]
def chunksOf(n):
    '''A series of lists of length n, subdividing the
       contents of xs. Where the length of xs is not evenly
       divible, the final list will be shorter than n.
    '''
    def go(xs):
        ys = list(xs)
        return (
            ys[i:n + i] for i in range(0, len(ys), n)
        ) if 0 < n else None
    return go


# concat :: [[a]] -> [a]
def concat(xs):
    '''The concatenation of all
       elements in a list.
    '''
    return [x for lst in xs for x in lst]


# findIndex :: (a -> Bool) -> [a] -> Maybe Int
def findIndex(p):
    '''Just the first index at which an
       element in xs matches p,
       or Nothing if no elements match.
    '''
    def go(xs):
        return next(
            (i for (i, x) in enumerate(xs) if p(x)),
            None
        )
    return go


# iterate :: (a -> a) -> a -> Gen [a]
def iterate(f):
    '''An infinite list of repeated
       applications of f to x.
    '''
    def go(x):
        v = x
        while True:
            yield v
            v = f(v)
    return go


# maybe :: b -> (a -> b) -> Maybe a -> b
def maybe(v):
    '''Either the default value v, if x is None,
       or the application of f to x.
    '''
    def go(f):
        def g(x):
            return v if None is x else f(x)
        return g
    return go


# MAIN ---
if __name__ == '__main__':
    main()

Output:

Cascade:
  4 3 3   0 4 3   1 0 4   1 1 0   2 1 0
  3 1 2 : 4 1 2 : 4 2 2 : 4 2 3 : 0 3 3
  0 2 3   0 2 3   0 2 3   0 2 3   1 2 3

s1 + s2 == s2 + s1 -> True
  1 2 0   2 1 3   3 3 3
  2 1 1 + 1 0 1 = 3 1 2
  0 1 3   0 1 0   0 2 3

  2 1 3   1 2 0   3 3 3
  1 0 1 + 2 1 1 = 3 1 2
  0 1 0   0 1 3   0 2 3

s3 + s3_id == s3 -> True
  3 3 3   2 1 2   3 3 3
  3 3 3 + 1 0 1 = 3 3 3
  3 3 3   2 1 2   3 3 3

s3_id + s3_id == s3_id -> True
  2 1 2   2 1 2   2 1 2
  1 0 1 + 1 0 1 = 1 0 1
  2 1 2   2 1 2   2 1 2

Raku

Works with: Rakudo version 2020.05

Most of the logic is lifted straight from the Abelian sandpile model task.

class ASP {
    has $.h = 3;
    has $.w = 3;
    has @.pile = 0 xx $!w * $!h;

    method topple {
        my $buf = $!w * $!h;
        my $done;
        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 self.pile[$idx] >= 4 {
                        my $grain = self.pile[$idx] div 4;
                        self.pile[ $idx - $!w ] += $grain if $row > 0;
                        self.pile[ $idx - 1  ]  += $grain if $idx - 1 >= $rs;
                        self.pile[ $idx + $!w ] += $grain if $row < $!h - 1;
                        self.pile[ $idx + 1  ]  += $grain if $idx + 1 < $re;
                        self.pile[ $idx ] %= 4;
                        $done = False;
                    }
                }
            }
        } until $done;
        self.pile;
    }
}

# some handy display layout modules
use Terminal::Boxer:ver<0.2+>;
use Text::Center;

for 3, (4,3,3,3,1,2,0,2,3), (2,1,2,1,0,1,2,1,2), # 3 square task example
    3, (2,1,2,1,0,1,2,1,2), (2,1,2,1,0,1,2,1,2), # 3 square identity
    5, (4,1,0,5,1,9,3,6,1,0,8,1,2,5,3,3,0,1,7,5,4,2,2,4,0), (2,3,2,3,2,3,2,1,2,3,2,1,0,1,2,3,2,1,2,3,2,3,2,3,2) # 5 square test
  -> $size, $pile, $identity {

    my $asp = ASP.new(:h($size), :w($size));

    $asp.pile = |$pile;

    my @display;

    my %p = :col($size), :3cw, :indent("\t");

    @display.push: rs-box |%p, |$identity;

    @display.push: rs-box |%p, $asp.pile;

    @display.push: rs-box |%p, $asp.topple;

    $asp.pile Z+= $identity.list;

    @display.push: rs-box |%p, $asp.pile;

    @display.push: rs-box |%p, $asp.topple;

    put %p<indent> ~ qww<identity 'test pile' toppled 'plus identity' toppled>».&center($size * 4 + 1).join: %p<indent>;

    .put for [Z~] @display».lines;

    put '';
}

Output:

	   identity  	  test pile  	   toppled   	plus identity	   toppled   
	╭───┬───┬───╮	╭───┬───┬───╮	╭───┬───┬───╮	╭───┬───┬───╮	╭───┬───┬───╮
	│ 2 │ 1 │ 2 │	│ 4 │ 3 │ 3 │	│ 2 │ 1 │ 0 │	│ 4 │ 2 │ 2 │	│ 2 │ 1 │ 0 │
	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤
	│ 1 │ 0 │ 1 │	│ 3 │ 1 │ 2 │	│ 0 │ 3 │ 3 │	│ 1 │ 3 │ 4 │	│ 0 │ 3 │ 3 │
	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤
	│ 2 │ 1 │ 2 │	│ 0 │ 2 │ 3 │	│ 1 │ 2 │ 3 │	│ 3 │ 3 │ 5 │	│ 1 │ 2 │ 3 │
	╰───┴───┴───╯	╰───┴───┴───╯	╰───┴───┴───╯	╰───┴───┴───╯	╰───┴───┴───╯

	   identity  	  test pile  	   toppled   	plus identity	   toppled   
	╭───┬───┬───╮	╭───┬───┬───╮	╭───┬───┬───╮	╭───┬───┬───╮	╭───┬───┬───╮
	│ 2 │ 1 │ 2 │	│ 2 │ 1 │ 2 │	│ 2 │ 1 │ 2 │	│ 4 │ 2 │ 4 │	│ 2 │ 1 │ 2 │
	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤
	│ 1 │ 0 │ 1 │	│ 1 │ 0 │ 1 │	│ 1 │ 0 │ 1 │	│ 2 │ 0 │ 2 │	│ 1 │ 0 │ 1 │
	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤	├───┼───┼───┤
	│ 2 │ 1 │ 2 │	│ 2 │ 1 │ 2 │	│ 2 │ 1 │ 2 │	│ 4 │ 2 │ 4 │	│ 2 │ 1 │ 2 │
	╰───┴───┴───╯	╰───┴───┴───╯	╰───┴───┴───╯	╰───┴───┴───╯	╰───┴───┴───╯

	       identity      	      test pile      	       toppled       	    plus identity    	       toppled       
	╭───┬───┬───┬───┬───╮	╭───┬───┬───┬───┬───╮	╭───┬───┬───┬───┬───╮	╭───┬───┬───┬───┬───╮	╭───┬───┬───┬───┬───╮
	│ 2 │ 3 │ 2 │ 3 │ 2 │	│ 4 │ 1 │ 0 │ 5 │ 1 │	│ 1 │ 3 │ 2 │ 1 │ 0 │	│ 3 │ 6 │ 4 │ 4 │ 2 │	│ 1 │ 3 │ 2 │ 1 │ 0 │
	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤
	│ 3 │ 2 │ 1 │ 2 │ 3 │	│ 9 │ 3 │ 6 │ 1 │ 0 │	│ 2 │ 2 │ 3 │ 3 │ 1 │	│ 5 │ 4 │ 4 │ 5 │ 4 │	│ 2 │ 2 │ 3 │ 3 │ 1 │
	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤
	│ 2 │ 1 │ 0 │ 1 │ 2 │	│ 8 │ 1 │ 2 │ 5 │ 3 │	│ 1 │ 1 │ 2 │ 0 │ 3 │	│ 3 │ 2 │ 2 │ 1 │ 5 │	│ 1 │ 1 │ 2 │ 0 │ 3 │
	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤
	│ 3 │ 2 │ 1 │ 2 │ 3 │	│ 3 │ 0 │ 1 │ 7 │ 5 │	│ 2 │ 0 │ 3 │ 2 │ 0 │	│ 5 │ 2 │ 4 │ 4 │ 3 │	│ 2 │ 0 │ 3 │ 2 │ 0 │
	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤	├───┼───┼───┼───┼───┤
	│ 2 │ 3 │ 2 │ 3 │ 2 │	│ 4 │ 2 │ 2 │ 4 │ 0 │	│ 3 │ 2 │ 3 │ 2 │ 1 │	│ 5 │ 5 │ 5 │ 5 │ 3 │	│ 3 │ 2 │ 3 │ 2 │ 1 │
	╰───┴───┴───┴───┴───╯	╰───┴───┴───┴───┴───╯	╰───┴───┴───┴───┴───╯	╰───┴───┴───┴───┴───╯	╰───┴───┴───┴───┴───╯

Red

Red [Purpose: "implement Abelian sandpile model"]

sadd: make object! [
	comb: function [pile1 [series!] pile2 [series!]] [
		repeat r 3 [
			repeat c 3 [
				pile2/:r/:c: pile2/:r/:c + pile1/:r/:c
			]
		]
		check pile2
	]
	check: func [pile [series!]] [
		stable: true row: col: none
		repeat r 3[
			repeat c 3[
				if pile/:r/:c >= 4 [
					stable: false 
					pile/:r/:c: pile/:r/:c - 4 
					row: r col: c 
					break]
			]
			if stable = false [break]
		]
		unless stable = false [print trim/with mold/only pile "[]" exit]
		spill pile row col
	]
	spill: func [pile [series!] r [integer!] c [integer!]] [
		neigh: reduce [
			right: reduce [r c - 1] up: reduce [r + 1 c]
			left: reduce [r c + 1] down: reduce [r - 1 c]	
		]
		foreach n neigh [
			unless any [(pile/(n/1) = none) (pile/(n/1)/(n/2) = none)] [
				pile/(n/1)/(n/2): pile/(n/1)/(n/2) + 1
			] 
		]
		check pile
	]
]

s1: [
	[1 2 0]
	[2 1 1]
	[0 1 3]
]

s2: [
	[2 1 3]
	[1 0 1]
	[0 1 0]
]

s3: [
	[3 3 3]
	[3 3 3]
	[3 3 3]
]

s3_id: [
	[2 1 2]
	[1 0 1]
	[2 1 2]
]

ex: [
	[4 3 3]
	[3 1 2]
	[0 2 3]
]

sadd/check copy/deep ex
sadd/comb copy/deep s1 copy/deep s2
sadd/comb copy/deep s2 copy/deep s1
sadd/comb copy/deep s3 copy/deep s3_id
sadd/comb copy/deep s3_id copy/deep s3_id

Output:

    2 1 0
    0 3 3
    1 2 3


    3 3 3
    3 1 2
    0 2 3


    3 3 3
    3 1 2
    0 2 3


    3 3 3
    3 3 3
    3 3 3


    2 1 2
    1 0 1
    2 1 2

REXX

/*REXX program demonstrates a 3x3 sandpile model by addition with toppling & avalanches.*/
@.= 0;   size= 3                                 /*assign 0 to all grid cells; grid size*/
call init   1,      1 2 0    2 1 1    0 1 3      /*   "   grains of sand──► sandpile 1. */
call init   2,      2 1 3    1 0 1    0 1 0      /*   "      "    "   "  "     "     2  */
call init   3,      3 3 3    3 3 3    3 3 3      /*   "      "    "   "  "     "     3  */
call init 's3_id',  2 1 2    1 0 1    2 1 2      /*   "      "    "   "  "     "    3_id*/
call show   1                                    /*display sandpile  1  to the terminal.*/
call show   2                                    /*   "        "     2   "  "      "    */
call  add   1, 2, 'sum1',           'adding sandpile  s1  and  s2  yields:'
call show 'sum1'
call  add   2, 1, 'sum2',           'adding sandpile  s2  and  s1  yields:'
call show 'sum2'
call  eq? 'sum1', 'sum2'                         /*is  sum1  the same as  sum2 ?        */
call show   3
call show 's3_id'
call  add   3,     's3_id', 'sum3', 'adding sandpile  s3  and  s3_id  yields:'
call show 'sum3'
call  add 's3_id', 's3_id', 'sum4', 'adding sandpile  s3_id  and  s3_id  yields:'
call show 'sum4'
call  eq? 'sum4',  's3_id'                       /*is  sum4  the same as  s3_id ?       */
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
@get: procedure expose @.; parse arg grid,r,c    ;                      return @.grid.r.c
@set: procedure expose @.; parse arg grid,r,c,val;  @.grid.r.c= val;    return
tran: procedure; parse arg a;  if datatype(a,'W')  then a='s'a;         return a
/*──────────────────────────────────────────────────────────────────────────────────────*/
add:  parse arg x, y, t;    if t==''  then t= 'sum';  xx= tran(x);  yy= tran(y)
        do r=1  for size;   do c=1  for size;   @.t.r.c= @.xx.r.c  +  @.yy.r.c
                            end   /*c*/
        end   /*r*/;    say arg(4);             call norm t;            return
/*──────────────────────────────────────────────────────────────────────────────────────*/
eq?:  parse arg x, y;   xx= tran(x);      yy= tran(y);      ?= 1
        do r=1  for size;         do c=1  for size;         ?= ?  &  (@.xx.r.c==@.yy.r.c)
                                  end   /*c*/
        end   /*r*/
      if ?  then say  'comparison of '      xx       " and "       yy':  same.'
            else say  'comparison of '      xx       " and "       yy':  not the same.'
      return
/*──────────────────────────────────────────────────────────────────────────────────────*/
init: parse arg x, $;   xx= tran(x);   #= 0;      pad= left('', 8);   ind= left('', 45)
        do r=1  for size;   do c=1  for size;     #= # + 1;           @.xx.r.c= word($, #)
                            end   /*c*/
        end   /*r*/;                              shows= 0;             return
/*──────────────────────────────────────────────────────────────────────────────────────*/
norm: procedure expose @. size; parse arg x;        xx= tran(x);            recurse= 0
        do r=1  for size;     do c=1  for size;     if @.xx.r.c<=size  then iterate
                              recurse= 1;              @.xx.r.c= @.xx.r.c - 4
                              call @set  xx, r-1, c  , @get(xx, r-1, c  )   +   1
                              call @set  xx, r+1, c  , @get(xx, r+1, c  )   +   1
                              call @set  xx, r  , c-1, @get(xx, r  , c-1)   +   1
                              call @set  xx, r  , c+1, @get(xx, r  , c+1)   +   1
                              end   /*c*/
        end   /*r*/;    if recurse  then call norm xx;                  return
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: parse arg x;  xx= tran(x);  say ind center("sandpile" xx,25,'─')  /*show the title*/
        do r=1  for size;  $=;      do c=1  for size;   $= $  @.xx.r.c  /*build a row.  */
                                    end   /*c*/
        say ind pad $                                                   /*display a row.*/
        end   /*r*/;     shows= shows + 1;     if shows==1  then say;   return

output   when using the default inputs:

(Shown at three-quarter size.)

                                              ──────sandpile s1  ──────
                                                        0 0 0
                                                        0 0 0
                                                        0 0 0

                                              ───────sandpile s2───────
                                                        2 1 3
                                                        1 0 1
                                                        0 1 0
adding sandpile  s1  and  s2  yields:
                                              ──────sandpile sum1──────
                                                        3 3 3
                                                        3 1 2
                                                        0 2 3
adding sandpile  s2  and  s1  yields:
                                              ──────sandpile sum2──────
                                                        3 3 3
                                                        3 1 2
                                                        0 2 3
comparison of  sum1  and  sum2:  same.
                                              ───────sandpile s3───────
                                                        3 3 3
                                                        3 3 3
                                                        3 3 3
                                              ─────sandpile s3_id──────
                                                        2 1 2
                                                        1 0 1
                                                        2 1 2
adding sandpile  s3  and  s3_id  yields:
                                              ──────sandpile sum3──────
                                                        3 3 3
                                                        3 3 3
                                                        3 3 3
adding sandpile  s3_id  and  s3_id  yields:
                                              ──────sandpile sum4──────
                                                        2 1 2
                                                        1 0 1
                                                        2 1 2
comparison of  sum4  and  s3_id:  same.

Ruby

class Sandpile
  
  def initialize(ar) = @grid = ar
 
  def to_a = @grid.dup
    
  def + (other)
    res = self.to_a.zip(other.to_a).map{|row1, row2| row1.zip(row2).map(&:sum) }
    Sandpile.new(res)
  end
  
  def stable? = @grid.flatten.none?{|v| v > 3}
  
  def avalanche 
    topple until stable?
    self
  end
  
  def == (other) = self.avalanche.to_a == other.avalanche.to_a
   
  def topple
    a = @grid
    a.each_index do |row|
       a[row].each_index do |col|
        next if a[row][col] < 4 
        a[row+1][col] += 1 unless row == a.size-1
        a[row-1][col] += 1 if row > 0
        a[row][col+1] += 1 unless col == a.size-1
        a[row][col-1] += 1 if col > 0
        a[row][col]   -= 4
      end
    end
    self
  end
    
  def to_s = "\n" + @grid.map {|row| row.join(" ") }.join("\n") 
      
end
    
puts "Sandpile:"
puts demo = Sandpile.new( [[4,3,3], [3,1,2],[0,2,3]] )
puts "\nAfter the avalanche:"
puts demo.avalanche
puts "_" * 30,""
  
s1 = Sandpile.new([[1, 2, 0], [2, 1, 1], [0, 1, 3]] )
puts "s1: #{s1}"
s2 = Sandpile.new([[2, 1, 3], [1, 0, 1], [0, 1, 0]] )
puts "\ns2: #{s2}"
puts "\ns1 + s2 == s2 + s1: #{s1 + s2 == s2 + s1}"
puts "_" * 30,""
  
s3    = Sandpile.new([[3, 3, 3], [3, 3, 3], [3, 3, 3]] )
s3_id = Sandpile.new([[2, 1, 2], [1, 0, 1], [2, 1, 2]] )
puts "s3 + s3_id == s3: #{s3 + s3_id == s3}"
puts "s3_id + s3_id == s3_id: #{s3_id + s3_id == s3_id}"

Output:

4 3 3
3 1 2
0 2 3

After the avalanche:

2 1 0
0 3 3
1 2 3
______________________________

s1: 
1 2 0
2 1 1
0 1 3

s2: 
2 1 3
1 0 1
0 1 0

s1 + s2 == s2 + s1: true
______________________________

s3 + s3_id == s3: true
s3_id + s3_id == s3_id: true

Rust

#[derive(Clone)]
struct Box {
    piles: [[u8; 3]; 3],
}

impl Box {
    fn init(piles: [[u8; 3]; 3]) -> Box {
        let a = Box { piles };

        if a.piles.iter().any(|&row| row.iter().any(|&pile| pile >= 4)) {
            return a.avalanche();
        } else {
            return a;
        }
    }

    fn avalanche(&self) -> Box {
        let mut a = self.clone();
        for (i, row) in self.piles.iter().enumerate() {
            for (j, pile) in row.iter().enumerate() {
                if *pile >= 4u8 {
                    if i > 0 {
                        a.piles[i - 1][j] += 1u8
                    }
                    if i < 2 {
                        a.piles[i + 1][j] += 1u8
                    }
                    if j > 0 {
                        a.piles[i][j - 1] += 1u8
                    }
                    if j < 2 {
                        a.piles[i][j + 1] += 1u8
                    }
                    a.piles[i][j] -= 4;
                }
            }
        }
        Box::init(a.piles)
    }

    fn add(&self, a: &Box) -> Box {
        let mut b = Box {
            piles: [[0u8; 3]; 3],
        };
        for (row, columns) in b.piles.iter_mut().enumerate() {
            for (col, pile) in columns.iter_mut().enumerate() {
                *pile = self.piles[row][col] + a.piles[row][col]
            }
        }
        Box::init(b.piles)
    }
}

fn main() {
    println!(
        "The piles demonstration avalanche starts as:\n{:?}\n{:?}\n{:?}",
        [4, 3, 3],
        [3, 1, 2],
        [0, 2, 3]
    );
    let s0 = Box::init([[4u8, 3u8, 3u8], [3u8, 1u8, 2u8], [0u8, 2u8, 3u8]]);
    println!(
        "And ends as:\n{:?}\n{:?}\n{:?}",
        s0.piles[0], s0.piles[1], s0.piles[2]
    );
    let s1 = Box::init([[1u8, 2u8, 0u8], [2u8, 1u8, 1u8], [0u8, 1u8, 3u8]]);
    let s2 = Box::init([[2u8, 1u8, 3u8], [1u8, 0u8, 1u8], [0u8, 1u8, 0u8]]);
    let s1_2 = s1.add(&s2);
    let s2_1 = s2.add(&s1);
    println!(
        "The piles in s1 + s2 are:\n{:?}\n{:?}\n{:?}",
        s1_2.piles[0], s1_2.piles[1], s1_2.piles[2]
    );
    println!(
        "The piles in s2 + s1 are:\n{:?}\n{:?}\n{:?}",
        s2_1.piles[0], s2_1.piles[1], s2_1.piles[2]
    );
    let s3 = Box::init([[3u8; 3]; 3]);
    let s3_id = Box::init([[2u8, 1u8, 2u8], [1u8, 0u8, 1u8], [2u8, 1u8, 2u8]]);
    let s4 = s3.add(&s3_id);
    println!(
        "The piles in s3 + s3_id are:\n{:?}\n{:?}\n{:?}",
        s4.piles[0], s4.piles[1], s4.piles[2]
    );
    let s5 = s3_id.add(&s3_id);
    println!(
        "The piles in s3_id + s3_id are:\n{:?}\n{:?}\n{:?}",
        s5.piles[0], s5.piles[1], s5.piles[2]
    );
}

Output:

The piles demonstration avalanche starts as:
[4, 3, 3]
[3, 1, 2]
[0, 2, 3]
And ends as:
[2, 1, 0]
[0, 3, 3]
[1, 2, 3]
The piles in s1 + s2 are:
[3, 3, 3]
[3, 1, 2]
[0, 2, 3]
The piles in s2 + s1 are:
[3, 3, 3]
[3, 1, 2]
[0, 2, 3]
The piles in s3 + s3_id are:
[3, 3, 3]
[3, 3, 3]
[3, 3, 3]
The piles in s3_id + s3_id are:
[2, 1, 2]
[1, 0, 1]
[2, 1, 2]

V (Vlang)

Translation of: Go

import strings

struct Sandpile {
mut:
    a [9]int
}
 
const (
neighbors = [
    [1, 3], [0, 2, 4], [1, 5], [0, 4, 6], [1, 3, 5, 7], [2, 4, 8], [3, 7], [4, 6, 8], [5, 7]
]
)
 
// 'a' is in row order
fn new_sandpile(a [9]int) Sandpile { return Sandpile{a} }
 
fn (s &Sandpile) plus(other &Sandpile) Sandpile {
    mut b := [9]int{}
    for i in 0..9 {
        b[i] = s.a[i] + other.a[i]
    }
    return Sandpile{b}
}
 
fn (s &Sandpile) is_stable() bool {
    for e in s.a {
        if e > 3 {
            return false
        }
    }
    return true
}
 
// just topples once so we can observe intermediate results
fn (mut s Sandpile) topple() {
    for i in 0..9 {
        if s.a[i] > 3 {
            s.a[i] -= 4
            for j in neighbors[i] {
                s.a[j]++
            }
            return
        }
    }
}
 
fn (s Sandpile) str() string {
    mut sb := strings.new_builder(64)
    for i in 0..3 {
        for j in 0..3 {
            sb.write_string("${u8(s.a[3*i+j])} ")
        }
        sb.write_string("\n")
    }
    return sb.str()
}
 
fn main() {
    println("Avalanche of topplings:\n")
    mut s4 := new_sandpile([4, 3, 3, 3, 1, 2, 0, 2, 3]!)
    println(s4)
    for !s4.is_stable() {
        s4.topple()
        println(s4)
    }
 
    println("Commutative additions:\n")
    s1 := new_sandpile([1, 2, 0, 2, 1, 1, 0, 1, 3]!)
    s2 := new_sandpile([2, 1, 3, 1, 0, 1, 0, 1, 0]!)
    mut s3_a := s1.plus(s2)
    for !s3_a.is_stable() {
        s3_a.topple()
    }
    mut s3_b := s2.plus(s1)
    for !s3_b.is_stable() {
        s3_b.topple()
    }
    println("$s1\nplus\n\n$s2\nequals\n\n$s3_a")
    println("and\n\n$s2\nplus\n\n$s1\nalso equals\n\n$s3_b")
 
    println("Addition of identity sandpile:\n")
    s3 := new_sandpile([3, 3, 3, 3, 3, 3, 3, 3, 3]!)
    s3_id := new_sandpile([2, 1, 2, 1, 0, 1, 2, 1, 2]!)
    s4 = s3.plus(s3_id)
    for !s4.is_stable() {
        s4.topple()
    }
    println("$s3\nplus\n\n$s3_id\nequals\n\n$s4")
 
    println("Addition of identities:\n")
    mut s5 := s3_id.plus(s3_id)
    for !s5.is_stable() {
        s5.topple()
    }
    print("$s3_id\nplus\n\n$s3_id\nequals\n\n$s5")
}

Output:

Same as Go entry

Wren

Library: Wren-fmt

import "./fmt" for Fmt

class Sandpile {
    static init() {
        __neighbors = [
            [1, 3], [0, 2, 4], [1, 5], [0, 4, 6], [1, 3, 5, 7], [2, 4, 8], [3, 7], [4, 6, 8], [5, 7]
        ]
    }

    // 'a' is a list of 9 integers in row order
    construct new(a) {
        _a = a
    }

    a { _a }

    +(other) {
        var b = List.filled(9, 0)
        for (i in 0..8) b[i] = _a[i] + other.a[i]
        return Sandpile.new(b)
    }

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

    // just topples once so we can observe intermediate results
    topple() {
        for (i in 0..8) {
            if (_a[i] > 3) {
                _a[i] = _a[i] - 4
                for (j in __neighbors[i]) _a[j] = _a[j] + 1
                return
            }
        }
    }

    toString {
        var s = ""
        for (i in 0..2) {
            for (j in 0..2) s = s + "%(a[3*i + j]) "
            s = s + "\n"
        }
        return s
    }
}

Sandpile.init()
System.print("Avalanche of topplings:\n")
var s4 = Sandpile.new([4, 3, 3, 3, 1, 2, 0, 2, 3])
System.print(s4)
while (!s4.isStable) {
    s4.topple()
    System.print(s4)
}

System.print("Commutative additions:\n")
var s1 = Sandpile.new([1, 2, 0, 2, 1, 1, 0, 1, 3])
var s2 = Sandpile.new([2, 1, 3, 1, 0, 1, 0, 1, 0])
var s3_a = s1 + s2
while (!s3_a.isStable) s3_a.topple()
var s3_b = s2 + s1
while (!s3_b.isStable) s3_b.topple()
Fmt.print("$s\nplus\n\n$s\nequals\n\n$s", s1, s2, s3_a)
Fmt.print("and\n\n$s\nplus\n\n$s\nalso equals\n\n$s", s2, s1, s3_b)

System.print("Addition of identity sandpile:\n")
var s3 = Sandpile.new(List.filled(9, 3))
var s3_id = Sandpile.new([2, 1, 2, 1, 0, 1, 2, 1, 2])
s4 = s3 + s3_id
while (!s4.isStable) s4.topple()
Fmt.print("$s\nplus\n\n$s\nequals\n\n$s", s3, s3_id, s4)

System.print("Addition of identities:\n")
var s5 = s3_id + s3_id
while (!s5.isStable) s5.topple()
Fmt.write("$s\nplus\n\n$s\nequals\n\n$s", s3_id, s3_id, s5)

Output:

Avalanche of topplings:

4 3 3 
3 1 2 
0 2 3 

0 4 3 
4 1 2 
0 2 3 

1 0 4 
4 2 2 
0 2 3 

1 1 0 
4 2 3 
0 2 3 

2 1 0 
0 3 3 
1 2 3 

Commutative additions:

1 2 0 
2 1 1 
0 1 3 

plus

2 1 3 
1 0 1 
0 1 0 

equals

3 3 3 
3 1 2 
0 2 3 

and

2 1 3 
1 0 1 
0 1 0 

plus

1 2 0 
2 1 1 
0 1 3 

also equals

3 3 3 
3 1 2 
0 2 3 

Addition of identity sandpile:

3 3 3 
3 3 3 
3 3 3 

plus

2 1 2 
1 0 1 
2 1 2 

equals

3 3 3 
3 3 3 
3 3 3 

Addition of identities:

2 1 2 
1 0 1 
2 1 2 

plus

2 1 2 
1 0 1 
2 1 2 

equals

2 1 2 
1 0 1 
2 1 2