I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Abelian sandpile model

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

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

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

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

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

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

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

11l

`V grid = [[0] * 10] * 10grid[5][5] = 64 print(‘Before:’)L(row) grid   print(row.map(c -> ‘#3’.format(c)).join(‘’)) F simulate(&grid)   L      V changed = 0B      L(arr) grid         V ii = L.index         L(val) arr            V jj = L.index            I val > 3               grid[ii][jj] -= 4               I ii > 0                  grid[ii - 1][jj]++               I ii < grid.len - 1                  grid[ii + 1][jj]++               I jj > 0                  grid[ii][jj - 1]++               I jj < grid.len - 1                  grid[ii][jj + 1]++               changed = 1B      I !changed         L.break simulate(&grid) print("\nAfter:")L(row) grid   print(row.map(c -> ‘#3’.format(c)).join(‘’)) grid = [[0] * 65] * 65grid[32][32] = 64 * 64 simulate(&grid) V ppm = File(‘sand_pile.ppm’, ‘w’)ppm.write_bytes(("P6\n#. #.\n255\n".format(grid.len, grid.len)).encode())V colors = [[Byte(0),   0, 0],            [Byte(255), 0, 0],            [Byte(0), 255, 0],            [Byte(0), 0, 255]]L(row) grid   L(c) row      ppm.write_bytes(colors[c])`
Output:
```Before:
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0 64  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0

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

AArch64 Assembly

Works with: as version Raspberry Pi 3B version Buster 64 bits or android 64 bits with application Termux
` /* ARM assembly AARCH64 Raspberry PI 3B or android 64 bits *//*  program abelian64.s   */  /* run : abelian 256 12 12  */ /*******************************************//* Constantes file                         *//*******************************************//* for this file see task include a file in language AArch64 assembly*/.include "../includeConstantesARM64.inc".equ MAXI, 25 /*********************************//* Initialized data              *//*********************************/.dataszMessValue:        .asciz "@ "szMessErrParam:     .asciz "error : command line = abelian size posx posy  \n"szMessFin:          .asciz "End display :\n"szCarriageReturn:   .asciz "\n" /*********************************//* UnInitialized data            *//*********************************/.bsssZoneConv:        .skip 24iSandPile:        .skip 8 * MAXI * MAXI/*********************************//*  code section                 *//*********************************/.text.global main main:                            // entry of program     mov fp,sp    ldr x4,[fp]                  // load number of parameters command line    cmp x4,#3                    // < 4 -> error    ble 99f    add x0,fp,#32                // load address param 4 = pos y    ldr x0,[x0]    bl conversionAtoD            // conversion ascii -> numeric    mov x3,x0    add x0,fp,#24                // load address param 3 = pos x    ldr x0,[x0]    bl conversionAtoD    mov x2,x0    add x0,fp,#16                 // load address param 2 = size begin pile    ldr x0,[x0]    bl conversionAtoD    ldr x4,qAdriSandPile    mov x5,#MAXI    madd x5,x3,x5,x2              // compute offset = maxi * y + x    str x0,[x4,x5,lsl #3]         // and store size in pos x,y    //mov x0,x4                   // display start position    //bl displaySandPile     mov x0,x4                     // sandpile address    mov x1,x2                     // pos x to start    mov x2,x3                     // pos y to start    bl addSand     ldr x0,qAdrszMessFin    bl affichageMess    mov x0,x4    bl displaySandPile    b 100f99:                               // line command error   ldr x0,qAdrszMessErrParam   bl affichageMess100:                              // 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 szCarriageReturnqAdrsZoneConv:            .quad sZoneConvqAdrszMessErrParam:       .quad szMessErrParamqAdrszMessFin:            .quad szMessFinqAdriSandPile:            .quad iSandPile/***************************************************//*     display  sandpile               *//***************************************************/// x0 contains address to sandpiledisplaySandPile:    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,#MAXI1:    mov x2,#0                   // indice x2:    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    retqAdrszMessValue:       .quad szMessValue/***************************************************//*     display  sandpile               *//***************************************************/// x0 contains address to sanspile// x1 contains position x// x2 contains position yaddSand:    stp x1,lr,[sp,-16]!         // save  registres    stp x2,x3,[sp,-16]!         // save  registres    stp x4,x5,[sp,-16]!         // save  registres    mov x3,#MAXI    madd x4,x3,x2,x1            // compute offset    ldr x5,[x0,x4,lsl #3]1:    cmp x5,#4                   // 4 grains ?    blt 100f    sub x5,x5,4                 // yes sustract    str x5,[x0,x4,lsl #3]    cmp x1,MAXI-1               // right position ok ?    beq 2f    add x1,x1,1                 // yes    bl add1Sand                 // add 1 grain    bl addSand                  // and compute new pile    sub x1,x1,12:    cmp x1,0                    // left position ok ?    beq 3f    sub x1,x1,1    bl add1Sand    bl addSand    add x1,x1,13:    cmp x2,0                    // higt position ok ?    beq 4f    sub x2,x2,1    bl add1Sand    bl addSand    add x2,x2,14:    cmp x2,MAXI-1               // low position ok ?    beq 5f    add x2,x2,1    bl add1Sand    bl addSand    sub x2,x2,15:   ldr x5,[x0,x4,lsl #3]       // reload value   b 1b                        // and loop100:    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 yadd1Sand:    stp x3,lr,[sp,-16]!       // save  registres    stp x4,x5,[sp,-16]!       // save  registres    mov x3,#MAXI    madd x4,x3,x2,x1          // compute offset    ldr x5,[x0,x4,lsl #3]     // load value at pos x,y    add x5,x5,1    str x5,[x0,x4,lsl #3]     // and store 100:    ldp x4,x5,[sp],16         // restaur des  2 registres    ldp x3,lr,[sp],16         // restaur des  2 registres    ret/********************************************************//*        File Include fonctions                        *//********************************************************//* for this file see task include a file in language AArch64 assembly */.include "../includeARM64.inc" `
Output:
```~/.../rosetta/asm1 \$ abelian64 64 12 12
End display :
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 2 2 2 0 2 2 2 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 2 2 2 2 2 1 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
```

ARM Assembly

Works with: as version Raspberry Pi
` /* ARM assembly Raspberry PI  or android 32 bits *//*  program abelian.s   */  /* run : abelian 256 12 12  */ /* REMARK 1 : this program use routines in a include file    see task Include a file language arm assembly    for the routine affichageMess conversion10    see at end of this program the instruction include *//* for constantes see task include a file in arm assembly *//************************************//* Constantes                       *//************************************/.include "../constantes.inc".equ MAXI, 25 /*********************************//* Initialized data              *//*********************************/.dataszMessValue:        .asciz "@ "szMessErrParam:     .asciz "error : command line = abelian size posx posy  \n"szMessFin:          .asciz "End display :\n"szCarriageReturn:   .asciz "\n" /*********************************//* UnInitialized data            *//*********************************/.bsssZoneConv:        .skip 24iSandPile:        .skip 4 * MAXI * MAXI/*********************************//*  code section                 *//*********************************/.text.global main main:                            @ entry of program     mov fp,sp    ldr r4,[fp]                  @ load number of parameters commend line    cmp r4,#3                    @ < 4 -> error    ble 99f    add r0,fp,#16                @ load address param 4 = pos y    ldr r0,[r0]    bl conversionAtoD            @ conversion ascii -> numeric    mov r3,r0    add r0,fp,#12                @ load address param 3 = pos x    ldr r0,[r0]    bl conversionAtoD    mov r2,r0    add r0,fp,#8                 @ load address param 2 = size begin pile    ldr r0,[r0]    bl conversionAtoD    ldr r4,iAdriSandPile    mov r5,#MAXI    mul r5,r3,r5                 @ compute offset = maxi * y    add r5,r2                    @ + x    str r0,[r4,r5,lsl #2]        @ and store size in pos x,y    //mov r0,r4                    @ display start position    //bl displaySandPile     mov r0,r4                   @ sandpile address    mov r1,r2                   @ pos x to start    mov r2,r3                   @ pos y to start    bl addSand     ldr r0,iAdrszMessFin    bl affichageMess    mov r0,r4    bl displaySandPile    b 100f99:                                  @ line command error   ldr r0,iAdrszMessErrParam   bl affichageMess100:                                  @ 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 szCarriageReturniAdrsZoneConv:            .int sZoneConviAdrszMessErrParam:       .int szMessErrParamiAdrszMessFin:            .int szMessFiniAdriSandPile:            .int iSandPile/***************************************************//*     display  sandpile               *//***************************************************/// r0 contains address to sandpiledisplaySandPile:    push {r1-r6,lr}             @ save  registers     mov r6,r0    mov r3,#0                   @ indice y    mov r4,#MAXI1:    mov r2,#0                   @ indice x2:    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                      @ returniAdrszMessValue:       .int szMessValue/***************************************************//*     display  sandpile               *//***************************************************/// r0 contains address to sanspile// r1 contains position x// r2 contains position yaddSand:    push {r1-r5,lr}             @ save  registers     mov r3,#MAXI    mul r4,r3,r2    add r4,r1    ldr r5,[r0,r4,lsl #2]1:    cmp r5,#4                   @ 4 grains ?    blt 100f    sub r5,#4                   @ yes sustract    str r5,[r0,r4,lsl #2]    cmp r1,#MAXI-1              @ right position ok ?    beq 2f    add r1,#1                   @ yes    bl add1Sand                 @ add 1 grain    bl addSand                  @ and compute new pile    sub r1,#12:    cmp r1,#0                   @ left position ok ?    beq 3f    sub r1,#1    bl add1Sand    bl addSand    add r1,#13:    cmp r2,#0                   @ higt position ok ?    beq 4f    sub r2,#1    bl add1Sand    bl addSand    add r2,#14:    cmp r2,#MAXI-1               @ low position ok ?    beq 5f    add r2,#1    bl add1Sand    bl addSand    sub r2,#15:   ldr r5,[r0,r4,lsl #2]       @ reload value   b 1b                        @ and loop100:    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 yadd1Sand:    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:
```[email protected]:~/asm32/rosetta32/ass10 \$ abelian 512 12 12
End display :
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  1  3  0  2  2  2  0  3  1  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  2  1  2  3  2  3  2  3  2  1  2  0  0  0  0  0  0  0
0  0  0  0  0  0  2  3  3  0  2  0  3  0  2  0  3  3  2  0  0  0  0  0  0
0  0  0  0  0  2  3  0  3  3  3  3  3  3  3  3  3  0  3  2  0  0  0  0  0
0  0  0  0  1  1  3  3  0  3  1  3  3  3  1  3  0  3  3  1  1  0  0  0  0
0  0  0  0  3  2  0  3  3  2  2  0  3  0  2  2  3  3  0  2  3  0  0  0  0
0  0  0  1  0  3  2  3  1  2  2  2  3  2  2  2  1  3  2  3  0  1  0  0  0
0  0  0  1  2  2  0  3  3  0  2  0  3  0  2  0  3  3  0  2  2  1  0  0  0
0  0  0  1  2  3  3  3  3  3  3  3  0  3  3  3  3  3  3  3  2  1  0  0  0
0  0  0  1  2  2  0  3  3  0  2  0  3  0  2  0  3  3  0  2  2  1  0  0  0
0  0  0  1  0  3  2  3  1  2  2  2  3  2  2  2  1  3  2  3  0  1  0  0  0
0  0  0  0  3  2  0  3  3  2  2  0  3  0  2  2  3  3  0  2  3  0  0  0  0
0  0  0  0  1  1  3  3  0  3  1  3  3  3  1  3  0  3  3  1  1  0  0  0  0
0  0  0  0  0  2  3  0  3  3  3  3  3  3  3  3  3  0  3  2  0  0  0  0  0
0  0  0  0  0  0  2  3  3  0  2  0  3  0  2  0  3  3  2  0  0  0  0  0  0
0  0  0  0  0  0  0  2  1  2  3  2  3  2  3  2  1  2  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  1  3  0  2  2  2  0  3  1  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
```

C

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

` #include<stdlib.h>#include<string.h>#include<stdio.h> int main(int argc, char** argv){	int i,j,sandPileEdge, centerPileHeight, processAgain = 1,top,down,left,right;		int** sandPile;	char* fileName;	static unsigned char colour[3]; 	if(argc!=3){		printf("Usage: %s <Sand pile side> <Center pile height>",argv[0]);		return 0;	} 	sandPileEdge = atoi(argv[1]);	centerPileHeight = atoi(argv[2]); 	if(sandPileEdge<=0 || centerPileHeight<=0){		printf("Sand pile and center pile dimensions must be positive integers.");		return 0;	} 	sandPile = (int**)malloc(sandPileEdge * sizeof(int*)); 	for(i=0;i<sandPileEdge;i++){		sandPile[i] = (int*)calloc(sandPileEdge,sizeof(int));	} 	sandPile[sandPileEdge/2][sandPileEdge/2] = centerPileHeight; 	printf("Initial sand pile :\n\n"); 	for(i=0;i<sandPileEdge;i++){		for(j=0;j<sandPileEdge;j++){			printf("%3d",sandPile[i][j]);		}		printf("\n");	} 	while(processAgain == 1){ 		processAgain = 0;		top = 0;		down = 0;		left = 0;		right = 0; 		for(i=0;i<sandPileEdge;i++){			for(j=0;j<sandPileEdge;j++){				if(sandPile[i][j]>=4){									if(i-1>=0){						top = 1;						sandPile[i-1][j]+=1;						if(sandPile[i-1][j]>=4)							processAgain = 1;					}					if(i+1<sandPileEdge){						down = 1;						sandPile[i+1][j]+=1;						if(sandPile[i+1][j]>=4)							processAgain = 1;					}					if(j-1>=0){						left = 1;						sandPile[i][j-1]+=1;						if(sandPile[i][j-1]>=4)							processAgain = 1;					}					if(j+1<sandPileEdge){						right = 1;						sandPile[i][j+1]+=1;						if(sandPile[i][j+1]>=4)							processAgain = 1;					}				sandPile[i][j] -= (top + down + left + right);				if(sandPile[i][j]>=4)					processAgain = 1;				}			}		}	} 	printf("Final sand pile : \n\n"); 	for(i=0;i<sandPileEdge;i++){		for(j=0;j<sandPileEdge;j++){			printf("%3d",sandPile[i][j]);		}		printf("\n");	} 	fileName = (char*)malloc((strlen(argv[1]) + strlen(argv[2]) + 23)*sizeof(char)); 	strcpy(fileName,"Final_Sand_Pile_");	strcat(fileName,argv[1]);	strcat(fileName,"_");	strcat(fileName,argv[2]);	strcat(fileName,".ppm"); 	FILE *fp = fopen(fileName,"wb"); 	fprintf(fp,"P6\n%d %d\n255\n",sandPileEdge,sandPileEdge); 	for(i=0;i<sandPileEdge;i++){		for(j=0;j<sandPileEdge;j++){			colour[0] = (sandPile[i][j] + i)%256;			colour[1] = (sandPile[i][j] + j)%256;			colour[2] = (sandPile[i][j] + i*j)%256;			fwrite(colour,1,3,fp);		}	} 	fclose(fp); 	printf("\nImage file written to %s\n",fileName); 	return 0;} `

Console output :

```[email protected]:~/doodles\$ ./a.out 10 64
Initial sand pile :

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

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

Image file written to Final_Sand_Pile_10_64.ppm
```

C++

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

`#include <iostream>#include "xtensor/xarray.hpp"#include "xtensor/xio.hpp"#include "xtensor-io/ximage.hpp" xt::xarray<int> init_grid (unsigned long x_dim, unsigned long y_dim){    xt::xarray<int>::shape_type shape = { x_dim, y_dim };    xt::xarray<int> grid(shape);     grid(x_dim/2, y_dim/2) = 64000;     return grid;} int print_grid(xt::xarray<int>& grid){    // for output to the terminal uncomment next line    // only makes sense for small grid < 32x32;    // std::cout << grid << std::endl << std::endl;     // output result to an image    xt::dump_image("grid.jpg", grid);     return 0;} bool iterate_grid(xt::xarray<int>& grid, const unsigned long& x_dim, const unsigned long& y_dim){    bool changed = false;     for (unsigned long i=0; i < x_dim; ++i)    {        for (unsigned long j=0; j < y_dim; ++j)        {            if ( grid(i, j) >= 4 )            {                grid(i, j) -= 4;                changed = true;                try                {                    grid.at(i-1, j) += 1;                    grid.at(i+1, j) += 1;                    grid.at(i, j-1) += 1;                    grid.at(i, j+1) += 1;                }                catch (const std::out_of_range& oor)                {                }            }        }    }     return changed;} int main(int argc, char* argv[]){    const unsigned long x_dim { 200 };    const unsigned long y_dim { 200 };     xt::xarray<int> grid = init_grid(x_dim, y_dim);    bool changed { true };     iterate_grid(grid, x_dim, y_dim);     while (changed == true)    {        changed = iterate_grid(grid, x_dim, y_dim);    }    print_grid(grid);     return 0;}`

Compile with following CMakeList.txt:

`cmake_minimum_required(VERSION 3.1)project(abelian_sandpile) find_package(xtl REQUIRED)find_package(xtensor REQUIRED)# if xtensor was built with xsimd support:# find_package(xsimd REQUIRED)set(CMAKE_CXX_FLAGS "\${CMAKE_CXX_FLAGS} -fopenmp")include_directories(/usr/include/OpenImageIO)find_library(OIIO "OpenImageIO") add_executable(abelian_sandpile src/abelian_sandpile.cpp) target_compile_options(abelian_sandpile PRIVATE -march=native -std=c++14)target_link_libraries(abelian_sandpile xtensor \${OIIO})`

Delphi

Translation of: Python
` program Abelian_sandpile_model; {\$APPTYPE CONSOLE} {\$R *.res} uses  System.SysUtils,  Vcl.Graphics,  System.Classes; type  TGrid = array of array of Integer; function Iterate(var Grid: TGrid): Boolean;var  changed: Boolean;  i: Integer;  j: Integer;  val: Integer;  Alength: Integer;begin  Alength := length(Grid);  changed := False;   for i := 0 to High(Grid) do    for j := 0 to High(Grid[0]) do    begin      val := Grid[i, j];      if val > 3 then      begin        Grid[i, j] := Grid[i, j] - 4;         if i > 0 then          Grid[i - 1, j] := Grid[i - 1, j] + 1;         if i < Alength - 1 then          Grid[i + 1, j] := Grid[i + 1, j] + 1;         if j > 0 then          Grid[i, j - 1] := Grid[i, j - 1] + 1;         if j < Alength - 1 then          Grid[i, j + 1] := Grid[i, j + 1] + 1;        changed := True;      end;    end;  Result := changed;end; procedure Simulate(var Grid: TGrid);var  changed: Boolean;begin  while Iterate(Grid) do    ;end; procedure Zeros(var Grid: TGrid; Size: Integer);var  i, j: Integer;begin  SetLength(Grid, Size, Size);  for i := 0 to Size - 1 do    for j := 0 to Size - 1 do      Grid[i, j] := 0;end; procedure Println(Grid: TGrid);var  i, j: Integer;begin  for i := 0 to High(Grid) do  begin    Writeln;    for j := 0 to High(Grid[0]) do      Write(Format('%3d', [Grid[i, j]]));  end;  Writeln;end; function Grid2Bmp(Grid: TGrid): TBitmap;const  Colors: array[0..2] of TColor = (clRed, clLime, clBlue);var  Alength: Integer;  i: Integer;  j: Integer;begin  Alength := Length(Grid);   Result := TBitmap.Create;  Result.SetSize(Alength, Alength);   for i := 0 to Alength - 1 do    for j := 0 to Alength - 1 do    begin      Result.Canvas.Pixels[i, j] := Colors[Grid[i, j]];    end;end; procedure Grid2P6(Grid: TGrid; FileName: TFileName);var  f: text;  i, j, Alength: Integer;  ppm: TFileStream;  Header: AnsiString;const  COLORS: array[0..3] of array[0..2] of byte = //  R,   G,    B((0   ,   0,    0), (255 ,   0,    0), (0   , 255,   0), (0   ,   0, 255));begin  Alength := Length(Grid);  ppm := TFileStream.Create(FileName, fmCreate);  Header := Format('P6'#10'%d %d'#10'255'#10, [Alength, Alength]);  writeln(Header);  ppm.Write(Tbytes(Header), Length(Header));   for i := 0 to Alength - 1 do    for j := 0 to Alength - 1 do    begin      ppm.Write(COLORS[Grid[i, j]], 3);    end;  ppm.Free;end; const  DIMENSION = 10; var  Grid: TGrid;  bmp: TBitmap; begin  Zeros(Grid, DIMENSION);  Grid[4, 4] := 64;  Writeln('Before:');  Println(Grid);   Simulate(Grid);   Writeln(#10'After:');  Println(Grid);   // Output bmp  with Grid2Bmp(Grid) do  begin    SaveToFile('output.bmp');    free;  end;   // Output ppm  Grid2P6(Grid, 'output.ppm');   Readln;end. `

Forth

Works with: gforth version 0.7.3

`#! /usr/bin/gforth -d 20M\ Abelian Sandpile Model 0 assert-level ! \ command-line : parse-number  s>number? invert throw drop ;: parse-size    ." size  : " next-arg parse-number dup . cr ;: parse-height  ." height: " next-arg parse-number dup . cr ;: parse-args    cr parse-size parse-height ; parse-args constant HEIGHT constant SIZE : allot-erase   create here >r dup allot r> swap erase ;: size^2        SIZE dup * cells ;: 2cells        [ 2 cells ] literal ;: -2cells       [ 2cells negate ] literal ; size^2 allot-erase arr \ array processing: ix            swap SIZE * + cells arr + ;: center        SIZE 2/ dup ;: write-cell    ix @ u. ;: write-row     SIZE 0 ?do dup i write-cell loop drop cr ;: arr.          SIZE 0 ?do i write-row loop ; \ stack processing : stack-empty?  dup -1 = ;: stack-full?   stack-empty? invert ; \ pgm-handling : concat        { a1 l1 a2 l2 } l1 l2 + allocate throw dup dup a1 swap l1 cmove a2 swap l1 + l2 cmove l1 l2 + ;: write-pgm     ." P2" cr SIZE u. SIZE u. cr ." 3" cr arr. ;: u>s           0 <# #s #> ;: filename      s" sandpile-" SIZE u>s concat s" -" concat HEIGHT u>s concat s" .pgm" concat ;: to-pgm        filename w/o create-file throw ['] write-pgm over outfile-execute close-file throw ; \ sandpile : prep-arr      HEIGHT center ix ! ;: prep-stack    -1 HEIGHT 4 u>= if center then ;: prepare       prep-arr prep-stack ;: ensure        if else 2drop 0 2rdrop exit then ;: col>=0        dup 0>= ensure ;: col<SIZE      dup SIZE < ensure ;: row>=0        over 0>= ensure ;: row<SIZE      over SIZE < ensure ;: legal?        col>=0 col<SIZE row>=0 row<SIZE 2drop true ;: north         1. d- ;: east          1+ ;: south         1. d+ ;: west          1- ;: reduce        2dup ix dup -4 swap +! @ 4 < if 2drop then ;: increase      2dup legal? if 2dup ix dup 1 swap +! @ 4 = if 2swap else 2drop then else 2drop then ; : inc-north     2dup north increase ;: inc-east      2dup east increase ;: inc-south     2dup south increase ;: inc-west      2dup west increase ;: inc-all       inc-north inc-east inc-south inc-west 2drop ;: simulate      prepare begin stack-full? while 2dup 2>r reduce 2r> inc-all repeat drop to-pgm ." written to " filename type cr ; simulate bye`
Output:

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

Fortran

Works with: gfortran version 9.2.0

The Abelian sandpile operations are defined here.

`module abelian_sandpile_m   implicit none   private  public :: pile   type :: pile    !! usage:    !!    1) init    !!    2) run     integer, allocatable :: grid(:,:)    integer              :: n(2)   contains    procedure :: init    procedure :: run     procedure, private :: process_node    procedure, private :: inside  end type contains   logical function inside(this, i)    class(pile), intent(in) :: this    integer,     intent(in) :: i(2)     inside = ((i(1) > 0) .and. (i(1) <= this%n(1)) .and. (i(2) > 0) .and. (i(2) <= this%n(2)) )  end function   recursive subroutine process_node(this, i)    !! start process     class(pile), intent(inout) :: this    integer,     intent(in)    :: i(2)      !! node coordinates to process     integer :: i0(2,2), j(2), d, k     ! if node has more than 4 grains -> redistribute    if (this%grid(i(1),i(2)) >= 4) then      ! unit vectors: help shift only one dimension (see below)      i0 = reshape([1,0,0,1], [2,2])       ! subtract 4 grains      this%grid(i(1),i(2)) = this%grid(i(1),i(2))-4       ! add one grain to neighbor if not out of bound      do d = 1, 2               ! loop dimensions        do k = -1, 1, 2         ! loop +-1 step in direction d          j = i+k*i0(:,d)       ! j = i, but one element is shifted by +-1          if (this%inside(j)) this%grid(j(1),j(2)) = this%grid(j(1),j(2)) + 1        end do      end do       ! check neighbor nodes      do d = 1, 2               ! loop dimensions        do k = -1, 1, 2         ! loop +-1 step in direction d          j = i+k*i0(:,d)       ! j = i, but one element is shifted by +-1          if (this%inside(j)) call this%process_node(j)        end do      end do       ! check itself      call this%process_node(i)    end if  end subroutine   subroutine run(this)    !! start process     class(pile), intent(inout) :: this     ! only node that could be unstable is inital node    call this%process_node(this%n/2)  end subroutine   subroutine init(this, nx, ny, h)    class(pile), intent(out) :: this    integer,     intent(in)  :: nx, ny      !! grid dimensions    integer,     intent(in)  :: h      !! height of and grains in middle of grid     this%n = [nx, ny]    allocate (this%grid(nx,ny), source=0)    this%grid(nx/2, ny/2) = h  end subroutine end module`

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

`program main   use rgbimage_m  use abelian_sandpile_m   implicit none   integer :: nx, ny, i, j   integer :: colors(0:3,3)   type(rgbimage) :: im  type(pile) :: p   colors(0,:) = [255,255,255]  colors(1,:) = [0,0,90]  colors(2,:) = [0,0,170]  colors(3,:) = [0,0,255]   nx = 200  ny = 100   call p%init(nx, ny, 2000)  call p%run   call im%init(nx, ny)   do i = 1, nx    do j = 1, ny      call im%set_pixel(i, j, colors(p%grid(i,j),:))    end do  end do   call im%write('fig.ppm') end program`

Fōrmulæ

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

Programs in Fōrmulæ are created/edited online in its website, However they run on execution servers. By default remote servers are used, but they are limited in memory and processing power, since they are intended for demonstration and casual use. A local server can be downloaded and installed, it has no limitations (it runs in your own computer). Because of that, example programs can be fully visualized and edited, but some of them will not run if they require a moderate or heavy computation/memory resources, and no local server is being used.

F#

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

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

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

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

Go

Translation of: Rust

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

`package main import (    "fmt"    "log"    "os"    "strings") const dim = 16 // image size func check(err error) {    if err != nil {        log.Fatal(err)    }} // Outputs the result to the terminal using UTF-8 block characters.func drawPile(pile [][]uint) {    chars:= []rune(" ░▓█")    for _, row := range pile {        line := make([]rune, len(row))        for i, elem := range row {            if elem > 3 { // only possible when algorithm not yet completed.                elem = 3            }            line[i] = chars[elem]        }        fmt.Println(string(line))    }} // Creates a .ppm file in the current directory, which contains// a colored image of the pile.func writePile(pile [][]uint) {    file, err := os.Create("output.ppm")    check(err)    defer file.Close()    // Write the signature, image dimensions and maximum color value to the file.    fmt.Fprintf(file, "P3\n%d %d\n255\n", dim, dim)    bcolors := []string{"125 0 25 ", "125 80 0 ", "186 118 0 ", "224 142 0 "}    var line strings.Builder    for _, row := range pile {                for _, elem := range row {            line.WriteString(bcolors[elem])        }        file.WriteString(line.String() + "\n")        line.Reset()     }} // Main part of the algorithm, a simple, recursive implementation of the model.func handlePile(x, y uint, pile [][]uint) {    if pile[y][x] >= 4 {        pile[y][x] -= 4        // Check each neighbor, whether they have enough "sand" to collapse and if they do,        // recursively call handlePile on them.        if y > 0 {            pile[y-1][x]++            if pile[y-1][x] >= 4 {                handlePile(x, y-1, pile)            }        }        if x > 0 {            pile[y][x-1]++            if pile[y][x-1] >= 4 {                handlePile(x-1, y, pile)            }        }        if y < dim-1 {            pile[y+1][x]++            if pile[y+1][x] >= 4 {                handlePile(x, y+1, pile)            }        }        if x < dim-1 {            pile[y][x+1]++            if pile[y][x+1] >= 4 {                handlePile(x+1, y, pile)            }        }         // Uncomment this line to show every iteration of the program.        // Not recommended with large input values.        // drawPile(pile)         // Finally call the function on the current cell again,        // in case it had more than 4 particles.        handlePile(x, y, pile)    }} func main() {    // Create 2D grid and set size using the 'dim' constant.    pile := make([][]uint, dim)    for i := 0; i < dim; i++ {        pile[i] = make([]uint, dim)    }     // Place some sand particles in the center of the grid and start the algorithm.    hdim := uint(dim/2 - 1)    pile[hdim][hdim] = 16    handlePile(hdim, hdim, pile)    drawPile(pile)     // Uncomment this to save the final image to a file    // after the recursive algorithm has ended.    // writePile(pile)}`
Output:
```

░
▓░▓
░░ ░░
▓░▓
░

```

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

Using a custom monad to make the code cleaner.

`{-# LANGUAGE FlexibleContexts           #-}{-# LANGUAGE GeneralizedNewtypeDeriving #-}{-# LANGUAGE ScopedTypeVariables        #-} module Rosetta.AbelianSandpileModel.ST     ( simulate    , test    , toPGM    ) where import Control.Monad.Reader (asks, MonadReader (..), ReaderT, runReaderT)import Control.Monad.ST (runST, ST)import Control.Monad.State (evalStateT, forM_, lift, MonadState (..), StateT, modify, when)import Data.Array.ST (freeze, readArray, STUArray, thaw, writeArray)import Data.Array.Unboxed (array, assocs, bounds, UArray, (!))import Data.Word (Word32)import System.IO (hPutStr, hPutStrLn, IOMode (WriteMode), withFile)import Text.Printf (printf) type Point     = (Int, Int)type ArrayST s = STUArray s Point Word32type ArrayU    = UArray Point Word32 newtype M s a = M (ReaderT (S s) (StateT [Point] (ST s)) a)    deriving (Functor, Applicative, Monad, MonadReader (S s), MonadState [Point]) data S s = S     { bMin :: !Point    , bMax :: !Point    , arr  :: !(ArrayST s)    } runM :: M s a -> S s -> [Point]-> ST s arunM (M m) = evalStateT . runReaderT m liftST :: ST s a -> M s aliftST = M . lift . lift simulate :: ArrayU -> ArrayUsimulate a = runST \$ simulateST a simulateST :: forall s. ArrayU -> ST s ArrayUsimulateST a = do    let (p1, p2) = bounds a        s = [p | (p, c) <- assocs a, c >= 4]    b <- thaw a :: ST s (ArrayST s)    let st = S { bMin = p1               , bMax = p2               , arr  = b               }    runM simulateM st s simulateM :: forall s. M s ArrayUsimulateM = do    ps <- get    case ps of        []      -> asks arr >>= liftST . freeze        p : ps' -> do            c <- changeArr p \$ \x -> x - 4            when (c < 4) \$ put ps'            forM_ [north, east, south, west] \$ inc . (\$ p)            simulateM changeArr :: Point -> (Word32 -> Word32) -> M s Word32changeArr p f = do    a    <- asks arr    oldC <- liftST \$ readArray a p    let newC = f oldC    liftST \$ writeArray a p newC    return newC inc :: Point -> M s ()inc p = do    b <- inBounds p    when b \$ do        c <- changeArr p succ        when (c == 4) \$ modify \$ (p :) inBounds :: Point -> M s BoolinBounds p = do    st <- ask    return \$ p >= bMin st && p <= bMax st north, east, south, west :: Point -> Pointnorth (x, y) = (x, y + 1)east  (x, y) = (x + 1, y)south (x, y) = (x, y - 1)west  (x, y) = (x - 1, y) toPGM :: ArrayU -> FilePath -> IO ()toPGM a fp = withFile fp WriteMode \$ \h -> do    let ((x1, y1), (x2, y2)) = bounds a        width  = x2 - x1 + 1        height = y2 - y1 + 1    hPutStrLn h "P2"    hPutStrLn h \$ show width ++ " " ++ show height    hPutStrLn h "3"    forM_ [y1 .. y2] \$ \y -> do        forM_ [x1 .. x2] \$ \x -> do            let c = min 3 \$ a ! (x, y)            hPutStr h \$ show c ++ " "        hPutStrLn h "" initArray :: Int -> Word32 -> ArrayUinitArray size height = array     ((-size, -size), (size, size))    [((x, y), if x == 0 && y == 0 then height else 0) | x <- [-size .. size], y <- [-size .. size]] test :: Int -> Word32 -> IO ()test size height = do    printf "size = %d, height = %d\n" size height    let a  = initArray size height        b  = simulate a        fp = printf "sandpile_%d_%d.pgm" size height    toPGM b fp    putStrLn \$ "wrote image to " ++ fp`
Output:

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

J

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

Java

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

`import java.awt.*;import java.awt.event.*;import javax.swing.*; public class AbelianSandpile {    public static void main(String[] args) {        SwingUtilities.invokeLater(new Runnable() {            public void run() {                Frame frame = new Frame();                frame.setVisible(true);            }        });    }     private static class Frame extends JFrame {        private Frame() {            super("Abelian Sandpile Model");            setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);            Container contentPane = getContentPane();            JPanel controlPanel = new JPanel(new FlowLayout(FlowLayout.LEFT));            JButton start = new JButton("Restart Simulation");            start.addActionListener(e -> restartSimulation());            JButton stop = new JButton("Stop Simulation");            stop.addActionListener(e -> stopSimulation());            controlPanel.add(start);            controlPanel.add(stop);            contentPane.add(controlPanel, BorderLayout.NORTH);            contentPane.add(canvas = new Canvas(), BorderLayout.CENTER);            timer = new Timer(100, e -> canvas.runAndDraw());            timer.start();            pack();        }         private void restartSimulation() {            timer.stop();            canvas.initGrid();            timer.start();        }         private void stopSimulation() {            timer.stop();        }         private Timer timer;        private Canvas canvas;    }     private static class Canvas extends JComponent {        private Canvas() {            setBorder(BorderFactory.createEtchedBorder());            setPreferredSize(new Dimension(600, 600));        }         public void paintComponent(Graphics g) {            int width = getWidth();            int height = getHeight();            g.setColor(Color.WHITE);            g.fillRect(0, 0, width, height);            int cellWidth = width/GRID_LENGTH;            int cellHeight = height/GRID_LENGTH;            for (int i = 0; i < GRID_LENGTH; ++i) {                for (int j = 0; j < GRID_LENGTH; ++j) {                    if (grid[i][j] > 0) {                        g.setColor(COLORS[grid[i][j]]);                        g.fillRect(i * cellWidth, j * cellHeight, cellWidth, cellHeight);                    }                }            }        }         private void initGrid() {            for (int i = 0; i < GRID_LENGTH; ++i) {                for (int j = 0; j < GRID_LENGTH; ++j) {                    grid[i][j] = 0;                }            }        }         private void runAndDraw() {            for (int i = 0; i < 100; ++i)                addSand(GRID_LENGTH/2, GRID_LENGTH/2);            repaint();        }         private void addSand(int i, int j) {            int grains = grid[i][j];            if (grains < 3) {                grid[i][j]++;            }            else {                grid[i][j] = grains - 3;                if (i > 0)                    addSand(i - 1, j);                if (i < GRID_LENGTH - 1)                    addSand(i + 1, j);                if (j > 0)                    addSand(i, j - 1);                if (j < GRID_LENGTH - 1)                    addSand(i, j + 1);            }        }         private int[][] grid = new int[GRID_LENGTH][GRID_LENGTH];    }     private static final Color[] COLORS = {        Color.WHITE,        new Color(0x00, 0xbf, 0xff),        new Color(0xff, 0xd7, 0x00),        new Color(0xb0, 0x30, 0x60)    };    private static final int GRID_LENGTH = 300;}`
Output:

See: abelian_sandpile.png (offsite PNG image)

Julia

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

`module AbelSand # supports output functionality for the results of the sandpile simulations# outputs the final grid in CSV format, as well as an image file using CSV, DataFrames, Images function TrimZeros(A)    # given an array A trims any zero rows/columns from its borders    # returns a 4 tuple of integers, i1, i2, j1, j2, where the trimmed array corresponds to A[i1:i2, j1:j2]    # A can be either numeric or a boolean array     i1, j1 = 1, 1    i2, j2 = size(A)     zz = typeof(A[1, 1])(0)    # comparison of a value takes into account the type as well     # i1 is the first row which has non zero element    for i = 1:size(A, 1)        q = false        for k = 1:size(A, 2)            if A[i, k] != zz                q = true                i1 = i                break            end        end         if q == true            break        end    end     # i2 is the first from below row with non zero element    for i in size(A, 1):-1:1        q = false        for k = 1:size(A, 2)            if A[i, k] != zz                q = true                i2 = i                break            end        end         if q == true            break        end    end     # j1 is the first column with non zero element     for j = 1:size(A, 2)        q = false        for k = 1:size(A, 1)            if A[k, j] != zz                j1 = j                q = true                break            end        end         if q == true            break        end    end     # j2 is the last column with non zero element     for j in size(A, 2):-1:1        q=false        for k=1:size(A,1)            if A[k, j] != zz                j2 = j                q=true                break            end        end         if q==true            break        end    end     return i1, i2, j1, j2end function addLayerofZeros(A, extraLayer)    # adds layer of zeros from all corners to the given array A     if extraLayer <= 0        return A    end     N, M = size(A)      Z = zeros( typeof(A[1,1]), N + 2*extraLayer, M + 2*extraLayer)    Z[(extraLayer+1):(N + extraLayer ), (extraLayer+1):(M+extraLayer)] = A     return Z end function printIntoFile(A, extraLayer, strFileName, TrimSmallValues = false)    # exports a 2d matrix A into a csv file    # @extraLayer is an integers adding layer of 0-s sorrounding the output matrix     # trimming off very small values; tiny values affect the performance of CSV export    if TrimSmallValues == true        A = map(x -> if (abs(x - floor(x)) < 0.01) floor(x) else x end, A)     end     i1, i2, j1, j2  = TrimZeros( A )    A = A[i1:i2, j1:j2]     A = addLayerofZeros(A, extraLayer)     CSV.write(string(strFileName,".csv"), DataFrame(A), writeheader = false)     return A end function Array_magnifier(A, cell_mag, border_mag)    # A is the main array; @cell_mag is the magnifying size of the cell,    # @border_mag is the magnifying size of the border between lattice cells     # creates a new array where each cell of the original array A appears magnified by size = cell_mag      total_factor = cell_mag + border_mag     A1 = zeros(typeof(A[1, 1]), total_factor*size(A, 1), total_factor*size(A, 2))     for i = 1:size(A,1), j = 1:size(A,2), u = ((i-1)*total_factor+1):(i*total_factor),                                          v = ((j-1)*total_factor+1):(j*total_factor)        if(( u - (i - 1) * total_factor <= cell_mag) && (v - (j - 1) * total_factor <= cell_mag))            A1[u, v] = A[i, j]         end    end     return A1 end function saveAsGrayImage(A, fileName, cell_mag, border_mag, TrimSmallValues = false)    # given a 2d matrix A, we save it as a gray image after magnifying by the given factors    A1 = Array_magnifier(A, cell_mag, border_mag)    A1 = A1/maximum(maximum(A1))     # trimming very small values from A1 to improve performance    if TrimSmallValues == true        A1 = map(x -> if ( x < 0.01) 0.0 else round(x, digits = 2) end, A1)     end     save(string(fileName, ".png") , colorview(Gray, A1)) end function saveAsRGBImage(A, fileName, color_codes, cell_mag, border_mag)    # color_codes is a dictionary, where key is a value in A and value is an RGB triplet    # given a 2d array A, and color codes (mapping from values in A to RGB triples), save A    # into fileName as png image after applying the magnifying factors     A1 = Array_magnifier(A, cell_mag, border_mag)    color_mat = zeros(UInt8, (3, size(A1, 1), size(A1, 2)))     for i = 1:size(A1,1)        for j = 1:size(A1,2)            color_mat[:, i, j]  = get(color_codes, A1[i, j] , [0, 0, 0])        end    end     save(string(fileName, ".png") , colorview(RGB, color_mat/255)) end const N_size = 700       # the radius of the lattice Z^2, the actual size becomes (2*N+1)x(2*N+1)const dx = [1, 0, -1, 0] # for a given (x,y) in Z^2, (x + dx, y + dy) for all (dx,dy) covers the neighborhood of (x,y)const dy = [0, 1, 0, -1] struct L_coord    # represents a lattice coordinate    x::Int    y::Intend function FindCoordinate(Z::Array{L_coord,1}, a::Int, b::Int)    # in the given array Z of coordinates finds the (first) index of the tuple (a,b)    # if no match, returns -1     for i=1:length(Z)        if (Z[i].x == a) && (Z[i].y == b)            return i        end    end     return -1end function move(N)    # the main function moving the pile sand grains of size N at the origin of Z^2 until the sandpile becomes stable     Z_lat = zeros(UInt8, 2 * N_size + 1, 2 * N_size + 1)     # models the integer lattice Z^2, we will have at most 4 sands on each vertex    V_sites = falses(2 * N_size + 1, 2 * N_size + 1)         # all sites which are visited by the sandpile process, are being marked here    Odometer = zeros(UInt64, 2 * N_size + 1, 2 * N_size + 1) # stores the values of the odometer function      walking = L_coord[]    # the coordinates of sites which need to move     V_sites[N_size + 1, N_size + 1] = true     # i1, ... j2  -> show the boundaries of the box which is visited by the sandpile process    i1, i2, j1, j2 = N_size + 1, N_size + 1, N_size + 1, N_size + 1     n = N     t1 = time_ns()     while n > 0        n -= 1         Z_lat[N_size + 1, N_size + 1] += 1        if (Z_lat[N_size + 1, N_size + 1] >= 4)            push!(walking, L_coord(N_size + 1, N_size + 1))        end         while(length(walking) > 0)            w = pop!(walking)            x = w.x            y = w.y             Z_lat[x, y] -= 4            Odometer[x, y] += 4             for k = 1:4                Z_lat[x + dx[k], y + dy[k]] += 1                V_sites[x + dx[k], y + dy[k]] = true                if Z_lat[x + dx[k], y + dy[k]] >= 4                    if FindCoordinate(walking, x + dx[k] , y + dy[k]) == -1                        push!(walking, L_coord( x + dx[k], y + dy[k]))                    end                end            end             i1 = min(i1, x - 1)            i2 = max(i2, x + 1)            j1 = min(j1, y - 1)            j2 = max(j2, y + 1)        end      end #end of the main while    t2 = time_ns()     println("The final boundaries are:: ", (i2 - i1 + 1),"x",(j2 - j1 + 1), "\n")    print("time elapsed: " , (t2 - t1) / 1.0e9, "\n")     Z_lat = printIntoFile(Z_lat, 0, string("Abel_Z_", N))    Odometer = printIntoFile(Odometer, 1, string("Abel_OD_", N))     saveAsGrayImage(Z_lat, string("Abel_Z_", N), 20, 0)    color_code = Dict(1=>[255, 128, 255], 2=>[255, 0, 0],3 => [0, 128, 255])    saveAsRGBImage(Z_lat, string("Abel_Z_color_", N), color_code, 20, 0)     # for the total elapsed time, it's better to use the @time macros on the main call     return Z_lat, Odometer # these are trimmed in output module end # end of function move  end # module  using .AbelSand Z_lat, Odometer = AbelSand.move(100000) `
Output:

Lua

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

Mathematica/Wolfram Language

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

Nim

Library: nimPNG

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

` # Abelian sandpile. from math import sqrtfrom nimPNG import savePNG24from sequtils import repeatfrom strformat import fmtfrom strutils import strip, addSep, parseInt # The grid represented as a sequence of sequences of int32.type Grid = seq[seq[int32]] # Colors to use for PPM and PNG files.const Colors = [[byte 100,  40,  15],                [byte 117,  87,  30],                [byte 181, 134,  47],                [byte 245, 182,  66]] #--------------------------------------------------------------------------------------------------- func sideLength(initVal: int32): int32 =  # Return the grid side length needed for "initVal" particles.  # We make sure that the returned value is odd.  result = sqrt(initVal.toFloat / 1.75).int32 + 3  result += result and 1 xor 1 #--------------------------------------------------------------------------------------------------- func doOneStep(grid: var Grid; boundary: var array[4, int]): bool =  ## Compute one step.   result = false   for y in boundary[0]..boundary[2]:    for x in boundary[1]..boundary[3]:      if grid[y][x] >= 4:         let rem = grid[y][x] div 4        grid[y][x] = grid[y][x] mod 4         if y - 1 >= 0:          inc grid[y - 1][x], rem          if y == boundary[0]:            dec boundary[0]         if x - 1 >= 0:          inc grid[y][x - 1], rem          if x == boundary[1]:            dec boundary[1]         if y + 1 < grid.len:          inc grid[y + 1][x], rem          if y == boundary[2]:            inc boundary[2]         if x + 1 < grid.len:          inc grid[y][x + 1], rem          if x == boundary[3]:            inc boundary[3]         result = true #--------------------------------------------------------------------------------------------------- proc display(grid: Grid; initVal: int) =  ## Display the grid as an array of values.   echo fmt"Starting with {initVal} particles."  echo ""   var line = newStringOfCap(2 * grid.len - 1)  for row in grid:    for value in row:      line.addSep(" ", 0)      line.add(\$value)    echo line    line.setLen(0)  echo "" #--------------------------------------------------------------------------------------------------- proc writePpmFile(grid: Grid; name: string) =  ## Write a grid representation in a PPM file.   var file = open(name, fmWrite)  file.write(fmt"P6 {grid.len} {grid.len} 255 ")   for row in grid:    for value in row:      discard file.writeBytes(Colors[value], 0, 3)   file.close()  echo fmt"PPM image written in ""{name}""." #--------------------------------------------------------------------------------------------------- proc writePngFile(grid: Grid; name: string) =  ## Write a grid representation in a PNG file.   var pixels = newSeq[byte](3 * grid.len * grid.len)   # Build pixel list.  var idx = 0  for row in grid:    for value in row:      pixels[idx..idx+2] = Colors[value]      inc idx, 3   discard savePNG24(name, pixels, grid.len, grid.len)  echo fmt"PNG image written in ""{name}""." #--------------------------------------------------------------------------------------------------- proc askInitVal(): int32 =  # Ask user for the number of particles.   while true:    stdout.write("Number of particles? ")    try:      let input = stdin.readLine().strip().parseInt()      if input in 4..int32.high:        return input.int32      echo fmt"Value not in expected range: 4..{int32.high}"    except ValueError:      echo "Invalid input"    except EOFError:      quit(QuitSuccess) #--------------------------------------------------------------------------------------------------- # Initialize the grid.let initVal = askInitVal()let sideLen = sideLength(initVal)var grid = repeat(newSeq[int32](sideLen), sideLen)let origin = grid.len div 2var boundaries: array[4, int] = [origin, origin, origin, origin]grid[origin][origin] = initVal # Run the simulation.while doOneStep(grid, boundaries):  discard # Display grid.if grid.len <= 20:  grid.display(initVal)#grid.writePpmFile(fmt"grid_{initVal}.ppm")grid.writePngFile(fmt"grid_{initVal}.png") `
Output:
```Number of particles? 100
Starting with 100 particles.

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

PNG image written in "grid_100.png".
```

Pascal

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

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

` program Abelian2;{\$IFDEF FPC}   {\$MODE DELPHI}{\$OPTIMIZATION ON,ALL}{\$CODEALIGN proc=16}{\$ALIGN 16}{\$ELSE}  {\$APPTYPE CONSOLE}{\$ENDIF}uses  SysUtils; type  Tlimit = record             lmtLow,LmtHigh : LongWord;           end;  TRowlimits = array of Tlimit;  tOneRow  = pLongWord;  tGrid = array of LongWord; var  Grid: tGrid;  Rowlimits:TRowlimits;  s : AnsiString;  maxval,maxCoor : NativeUint; function CalcMaxCoor(maxVal : NativeUint):NativeUint;//  maxVal = 10000;maxCoor = 77-2;// maxCoor*maxCoor    *1,778;     0.009sec//  maxVal = 100000;maxCoor = 236-2;// maxCoor*maxCoor  *1.826;     0.825sec//  maxVal = 1000000;maxCoor = 732-2;// maxCoor*maxCoor *1.877;    74    secBegin  result := trunc(sqrt(maxval/1.75))+3;end; procedure clear;begin  setlength(Grid,0);  setlength(Rowlimits,0);  s := '';end; procedure InitGrid(var G:tGrid;InitVal:NativeUint);var  row,middle: nativeINt;begin//  setlength(Rowlimits,0);   setlength(G,0);  MaxCoor :=  CalcMaxCoor(InitVal);  setlength(G,sqr(maxCoor));  setlength(Rowlimits,maxCoor);  fillchar(G[0],length(G)*SizeOf(G[0]),#0);   middle := (maxCoor) div 2;  Grid[middle*maxcoor+middle] := InitVal;  For row := 1 to maxCoor do    with Rowlimits[row] do    Begin      lmtLow := middle;      lmtHigh := middle;    end;   with Rowlimits[middle] do  Begin    lmtLow := middle;    lmtHigh := middle;  end;end;procedure OutGridPPM(const G:tGrid;maxValue : NativeUint);const  color : array[0..3] of array[0..2] of Byte =             //R,G,B)            ((0,0,0),             (255,0,0),             (0,255,0),             (0,0,255));var  f :text;  pActRow: tOneRow;  col,row,sIdx,value : NativeInt;Begin  Assignfile(f,'ppm/Grid_'+IntToStr(maxValue)+'.ppm');  rewrite(f);  write(f,Format('P6 %d %d %d ',[maxCoor-1,maxCoor-1,255]));  setlength(s,(maxCoor-1)*3);  pActRow :=@G[0];  For row := maxCoor-2 downto 0 do  Begin    inc(pActRow,maxCoor);    sIdx := 1;    For col := 1 to maxCoor-1 do    Begin      value := pActRow[col];      s[sIdx]   := CHR(color[value,0]);      s[sIdx+1] := CHR(color[value,1]);      s[sIdx+2] := CHR(color[value,2]);      inc(sIdx,3);    end;    write(f,s);  end;  CloseFile(f);end; procedure OutGrid(const G:tGrid);//output of grid and test, if no sand is lostvar  pActRow: tOneRow;  col,row,sum,value : NativeUint;Begin  setlength(s,maxcoor-1);  pActRow := @G[0];  sum := 0;  For row := maxCoor-1 downto 1 do  Begin    inc(pActRow,maxcoor);    For col := 1 to maxCoor-1 do    Begin      value := pActRow[col];//      IF value>=4 then writeln(row:5,col:5,value:13);      s[col] := chr(value+48);      inc(sum,value);    end;    if maxCoor <80 then      writeln(s);  end;  writeln('columns ',maxcoor-1,' checksum ',maxVal,' ?=? ',sum);{  For row := 1 to maxCoor do    with Rowlimits[row] do      writeln(lmtLow:10,lmtHigh:10);      * }end; procedure Evolution(var G:tGrid);var  pActRow,pRowBefore,pRowAfter : tOneRow;  col,row,mul,val,done : NativeUint;begin  repeat    pRowBefore := @G[0];    pActRow    := @G[maxcoor];    pRowAfter  := @G[2*maxcoor];    done := 0;    For row := maxCoor-1 downto 1 do    Begin      with RowLimits[row] do      Begin      while (LmtLow >1) AND (pActRow[lmtLow]<> 0) do        dec(lmtLow);      while (lmtHigh < maxCoor) AND (pActRow[lmtHigh]<> 0) do        inc(lmtHigh);      For col := lmtLow to lmtHigh do      Begin        val := pActRow[col];        IF val >=4 then        Begin          mul := val DIV 4;          done := val;          inc(pRowBefore[col],mul);          inc(pActRow[col-1],mul);          pActRow[col] := val-4*Mul;          inc(pActRow[col+1],mul);          inc(pRowAfter[col],mul);        end;      end;      pRowBefore:= pActRow;      pActRow := pRowAfter;      inc(pRowAfter,maxcoor);    end;    end;  until done=0;end; procedure OneTurn(count:NativeUint);begin  Writeln(' Test abelian sandpile( ',count,' )');  MaxVal := count;  InitGrid(Grid,count);  Evolution(Grid);  OutGrid(Grid);  OutGridPPM(Grid,count);  clear;end; BEGIN  OneTurn(4);  OneTurn(16);  OneTurn(64);  OneTurn(1000);  OneTurn(10000);  OneTurn(100000);END. `
Output:
``` Test abelian sandpile( 4 )
010
101
010
columns 3 checksum 4 ?=? 4
Test abelian sandpile( 16 )
00100
02120
11011
02120
00100
columns 5 checksum 16 ?=? 16
Test abelian sandpile( 64 )
00121000
02222200
12222210
22202220
12222210
02222200
00121000
00000000
columns 8 checksum 64 ?=? 64
Test abelian sandpile( 1000 )
0000000001111111000000000
0000000130233320310000000
0000013223313133223100000
0000213222130312223120000
0002220123332333210222000
0011223233123213323221100
0033032313221223132303300
0122123203311133023212210
0322231023333333201322230
1032333332231322333332301
1231312332232322332131321
1313322133322233312233131
1330231131220221311320331
1313322133322233312233131
1231312332232322332131321
1032333332231322333332301
0322231023333333201322230
0122123203311133023212210
0033032313221223132303300
0011223233123213323221100
0002220123332333210222000
0000213222130312223120000
0000013223313133223100000
0000000130233320310000000
0000000001111111000000000
columns 25 checksum 1000 ?=? 1000
Test abelian sandpile( 10000 )
--shortened
columns 77 checksum 10000 ?=? 10000
Test abelian sandpile( 100000 )
columns 241 checksum 100000 ?=? 100000

real    0m0,815s
```

OCaml

Works with: OCaml version 4.11

In Sandpile module (sandpile.ml)

` module Make =  functor (M : sig val m : int val n : int end)  -> struct     let grid = Array.init M.m (fun _ -> Array.make M.n 0)     let print () =      for i = 0 to M.m - 1      do for j = 0 to M.n - 1         do Printf.printf "%d " grid.(i).(j)         done       ; print_newline ()       done     let unstable = Hashtbl.create 10     let add_grain x y      = grid.(x).(y) <- grid.(x).(y) + 1      ; if grid.(x).(y) >= 4 then          Hashtbl.replace unstable (x,y) () (* Use Hashtbl.replace for uniqueness *)     let topple x y      = grid.(x).(y) <- grid.(x).(y) - 4      ; if grid.(x).(y) < 4        then Hashtbl.remove unstable (x,y)      ; match (x,y) with        (* corners *)        | (0,0) -> add_grain 1 0                 ; add_grain 0 1        | (0,n) when n = M.n - 1          -> add_grain 1 n           ; add_grain 0 (n-1)        | (m,0) when m = M.m - 1          -> add_grain m 1           ; add_grain (m-1) 0        | (m,n) when m = M.m - 1 && n = M.n - 1          -> add_grain ( m ) (n-1)           ; add_grain (m-1) ( n )        (* sides *)        | (0,y) -> add_grain 1 y                 ; add_grain 0 (y+1)                 ; add_grain 0 (y-1)        | (m,y) when m = M.m - 1          -> add_grain ( m ) (y-1)           ; add_grain ( m ) (y+1)           ; add_grain (m-1) ( y )        | (x,0) -> add_grain (x+1) 0                 ; add_grain (x-1) 0                 ; add_grain ( x ) 1        | (x,n) when n = M.n - 1          -> add_grain (x-1) ( n )           ; add_grain (x+1) ( n )           ; add_grain ( x ) (n-1)        (* else *)        | (x,y) -> add_grain ( x ) (y+1)                 ; add_grain ( x ) (y-1)                 ; add_grain (x+1) ( y )                 ; add_grain (x-1) ( y )     let add_sand n x y      = for i = 1 to n        do add_grain x y        done     let avalanche ()      = while Hashtbl.length unstable > 0         do         let unstable' = Hashtbl.fold  (fun (x,y) () r -> (x,y) :: r) unstable []         in List.iter (fun (x,y) -> topple x y ) unstable'        done   end (* testing *) let ()  = let module S = Make (struct let m = 11 let n = 11 end)    in S.add_sand 500 5 5     ; S.avalanche ()     ; S.print () `

Perl

`#!/usr/bin/perl use strict; # http://www.rosettacode.org/wiki/Abelian_sandpile_modeluse warnings; my (\$high, \$wide) = split ' ', qx(stty size);my \$mask = "\0" x \$wide . ("\0" . "\177" x (\$wide - 2) . "\0") x (\$high - 5) .  "\0" x \$wide;my \$pile = \$mask =~ s/\177/ rand() < 0.02 ? chr 64 + rand 20 : "\0" /ger; for ( 1 .. 1e6 )  {  print "\e[H", \$pile =~ tr/\0-\177/ 1-~/r, "\n\$_";  my \$add = \$pile =~ tr/\1-\177/\0\0\0\200/r; # set high bit for >=4  \$add =~ /\200/ or last;  \$pile =~ tr/\4-\177/\0-\173/; # subtract 4 if >=4  for ("\0\$add", "\0" x \$wide . \$add, substr(\$add, 1), substr \$add, \$wide)    {    \$pile |= \$_;    \$pile =~ tr/\200-\377/\1-\176/; # add one to each neighbor of >=4    \$pile &= \$mask;    }  select undef, undef, undef, 0.1; # comment out for full speed  }`

Phix

Library: Phix/pGUI

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

```-- demo\rosetta\Abelian_sandpile_model.exw
include pGUI.e

Ihandle dlg, canvas
cdCanvas cddbuffer

sequence board = {{0,0,0},
{0,0,0},
{0,0,0}}

procedure drop(integer y, x)
sequence moves = {}
while true do
board[y,x] += 1
if board[y,x]>=4 then
board[y,x] -= 4
moves &= {{y,x-1},{y,x+1},{y-1,x},{y+1,x}}
end if
-- extend board if rqd (maintain a border of zeroes)
if x=1 then                             -- extend left
for i=1 to length(board) do
board[i] = prepend(board[i],0)
end for
for i=1 to length(moves) do
moves[i][2] += 1
end for
elsif x=length(board[1]) then           -- extend right
for i=1 to length(board) do
board[i] = append(board[i],0)
end for
end if
-- (copy the all-0 lines from the other end...)
if y=1 then                             -- extend up
board = prepend(board,board[\$])
for i=1 to length(moves) do
moves[i][1] += 1
end for
elsif y=length(board) then              -- extend down
board = append(board,board[1])
end if
if length(moves)=0 then exit end if
{y,x} = moves[\$]
moves = moves[1..\$-1]
end while
IupUpdate(canvas)
end procedure

function timer_cb(Ihandle /*ih*/)
integer y = floor(length(board)/2)+1,
x = floor(length(board[1])/2)+1
drop(y,x)
return IUP_DEFAULT
end function

function redraw_cb(Ihandle ih, integer /*posx*/, integer /*posy*/)
IupGLMakeCurrent(ih)
cdCanvasActivate(cddbuffer)
cdCanvasClear(cddbuffer)
for y=1 to length(board) do
for x=1 to length(board[1]) do
integer c = board[y][x]
if c!=0 then
integer colour = {CD_VIOLET,CD_RED,CD_BLUE}[c]
cdCanvasPixel(cddbuffer, x, y, colour)
end if
end for
end for
cdCanvasFlush(cddbuffer)
return IUP_DEFAULT
end function

function map_cb(Ihandle ih)
IupGLMakeCurrent(ih)
atom res = IupGetDouble(NULL, "SCREENDPI")/25.4
cddbuffer = cdCreateCanvas(CD_GL, "300x100 %g", {res})
cdCanvasSetBackground(cddbuffer, CD_PARCHMENT)
return IUP_DEFAULT
end function

procedure main()
IupOpen()
canvas = IupGLCanvas("RASTERSIZE=300x100")
IupSetCallbacks({canvas}, {"ACTION", Icallback("redraw_cb"),
"MAP_CB", Icallback("map_cb")})
dlg = IupDialog(canvas,"TITLE=\"Abelian sandpile model\"")
IupShow(dlg)
Ihandle timer = IupTimer(Icallback("timer_cb"), 10)
IupMainLoop()
IupClose()
end procedure

main()
```

Python

Python: Original, with output

` import numpy as npimport matplotlib.pyplot as plt  def iterate(grid):    changed = False    for ii, arr in enumerate(grid):        for jj, val in enumerate(arr):            if val > 3:                grid[ii, jj] -= 4                if ii > 0:                    grid[ii - 1, jj] += 1                if ii < len(grid)-1:                    grid[ii + 1, jj] += 1                if jj > 0:                    grid[ii, jj - 1] += 1                if jj < len(grid)-1:                    grid[ii, jj + 1] += 1                changed = True    return grid, changed  def simulate(grid):    while True:        grid, changed = iterate(grid)        if not changed:            return grid  if __name__ == '__main__':    start_grid = np.zeros((10, 10))    start_grid[4:5, 4:5] = 64    final_grid = simulate(start_grid.copy())    plt.figure()    plt.gray()    plt.imshow(start_grid)    plt.figure()    plt.gray()    plt.imshow(final_grid) `

Output: </n> Before:

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

After:

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

An interactive variant to the above solution:

` from os import system, namefrom time import sleep def clear():	if name == 'nt':		_ = system('cls')	else: _ = system('clear') def exit():	import sys	clear()	sys.exit() def make_area(x, y):	area = [[0]*x for _ in range(y)]	return area def make_sandpile(area, loc, height):	loc=list(n-1 for n in loc)	x, y = loc 	try: area[y][x]+=height	except IndexError: pass def run(area, by_frame=False):	def run_frame():		for y_index, group in enumerate(area):			y = y_index+1 			for x_index, height in enumerate(group):				x = x_index+1 				if height < 4: continue 				else:					make_sandpile(area, (x+1, y), 1)					make_sandpile(area, (x, y+1), 1) 					if x_index-1 >= 0:						make_sandpile(area, (x-1, y), 1)					if y_index-1 >= 0:						make_sandpile(area, (x, y-1), 1) 					make_sandpile(area, (x, y), -4) 	while any([any([pile>=4 for pile in group]) for group in area]):		if by_frame:			clear()		run_frame()		if by_frame:			show_area(area); sleep(.05) def show_area(area):	display = [' '.join([str(item) if item else ' ' for item in group]) 			   for group in area]	[print(i) for i in display] clear()if __name__ == '__main__':	area = make_area(10, 10)	print('\nBefore:')	show_area(area)	make_sandpile(area, (5, 5), 64)	run(area)	print('\nAfter:')	show_area(area) `

Output:

` Before:0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 After:0 0 0 0 0 0 0 0 0 00 0 0 1 2 1 0 0 0 00 0 2 2 2 2 2 0 0 00 1 2 2 2 2 2 1 0 00 2 2 2 0 2 2 2 0 00 1 2 2 2 2 2 1 0 00 0 2 2 2 2 2 0 0 00 0 0 1 2 1 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 `

Python: using tkinter

` ''' Python 3.6.5 code using Tkinter graphical user    interface (Canvas widget) to display final results.'''from tkinter import * # given a grid, display it on a tkinter Canvas:class Sandpile:    def __init__(self, wn, grid):        self.window = wn        self.grid = grid        self.canvas = Canvas(wn, bg='lemon chiffon')        self.canvas.pack(fill=BOTH, expand=1)                colors = {0:'dodger blue',                  1:'red',                  2:'green',                  3:'lemon chiffon'}         x = 10        y = 10        d = 5         for row in self.grid:            for value in row:                clr = colors[value]                self.canvas.create_rectangle(                    x, y, x+d, y+d,                    outline=clr,                    fill = clr)                x += 5            x = 10            y += 5 class Grid:    def __init__(self, size, center):        self.size = size        # rows/cols in (best if odd)        self.center = center    # start value at center of grid        self.grid = [[0]*self.size for i in range(self.size)]        self.grid[self.size // 2][self.size // 2] = self.center     # print the grid:    def show(self, msg):        print('  ' + msg + ':')        for row in self.grid:            print(' '.join(str(x) for x in row))               print()        return     # dissipate piles of sand as required:    def abelian(self):        while True:            found = False            for r in range(self.size):                for c in range(self.size):                    if self.grid[r][c] > 3:                        self.distribute(self.grid[r][c], r, c)                        found = True            if not found:                return     # distribute sand from a single pile to its neighbors:    def distribute(self, nbr, row, col):        qty, remain = divmod(nbr, 4)        self.grid[row][col] = remain        for r, c in [(row+1, col),                     (row-1, col),                     (row, col+1),                     (row, col-1)]:            self.grid[r][c] += qty        return     # display the grid using tkinter:    def display(self):        root = Tk()        root.title('Sandpile')        root.geometry('700x700+100+50')        sp = Sandpile(root, self.grid)        root.mainloop() # execute program for size, center value pair:# just print results for a small gridg = Grid(9,17)g.show('BEFORE')g.abelian()          # scatter the sandg.show('AFTER') # just show results in tkinter for a large grid# I wish there was a way to attach a screen shot# of the tkinter result hereg = Grid(131,25000)g.abelian()          # scatter the sandg.display()          # display result using tkinter ##  OUTPUT:####      BEFORE:##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 17 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0####      AFTER:##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 1 0 0 0 0##    0 0 0 2 1 2 0 0 0##    0 0 1 1 1 1 1 0 0##    0 0 0 2 1 2 0 0 0##    0 0 0 0 1 0 0 0 0##    0 0 0 0 0 0 0 0 0##    0 0 0 0 0 0 0 0 0 `

Raku

(formerly Perl 6)

Terminal based

Works with: Rakudo version 2019.07.1

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

`sub cleanup { print "\e[0m\e[?25h\n"; exit(0) } signal(SIGINT).tap: { cleanup(); exit(0) } unit sub MAIN (\$stack = 1000, :\$hide-progress = False ); my @color = "\e[38;2;0;0;0m█",            "\e[38;2;255;0;0m█",            "\e[38;2;255;255;0m█",            "\e[38;2;0;0;255m█",            "\e[38;2;255;255;255m█"            ; my (\$h, \$w) = qx/stty size/.words».Int;my \$buf = \$w * \$h;my @buffer = 0 xx \$buf;my \$done; @buffer[\$w * (\$h div 2) + (\$w div 2) - 1] = \$stack; print "\e[?25l\e[48;5;232m"; repeat {    \$done = True;    loop (my int \$row; \$row < \$h; \$row = \$row + 1) {        my int \$rs = \$row * \$w; # row start        my int \$re = \$rs  + \$w; # row end        loop (my int \$idx = \$rs; \$idx < \$re; \$idx = \$idx + 1) {            if @buffer[\$idx] >= 4 {                my \$grains = @buffer[\$idx] div 4;                @buffer[ \$idx - \$w ] += \$grains if \$row > 0;                @buffer[ \$idx - 1  ] += \$grains if \$idx - 1 >= \$rs;                @buffer[ \$idx + \$w ] += \$grains if \$row < \$h - 1;                @buffer[ \$idx + 1  ] += \$grains if \$idx + 1 < \$re;                @buffer[ \$idx ] %= 4;                \$done = False;            }        }    }    unless \$hide-progress {        print "\e[1;1H", @buffer.map( { @color[\$_ min 4] }).join;    }} until \$done; print "\e[1;1H", @buffer.map( { @color[\$_ min 4] }).join; cleanup;`

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

SDL2 Animation

`use NativeCall;use SDL2::Raw; my (\$width, \$height) = 900, 900; unit sub MAIN (\$stack = 10000); my int (\$w, \$h) = 160, 160; my \$buf = \$w * \$h;my @buffer = 0 xx \$buf; @buffer[\$w * (\$h div 2) + (\$w div 2) - 1] = \$stack;  SDL_Init(VIDEO); my SDL_Window \$window = SDL_CreateWindow(    "Abelian sandpile - Raku",    SDL_WINDOWPOS_CENTERED_MASK, SDL_WINDOWPOS_CENTERED_MASK,    \$width, \$height,    RESIZABLE); my SDL_Renderer \$renderer = SDL_CreateRenderer( \$window, -1, ACCELERATED +| TARGETTEXTURE ); my \$asp_texture = SDL_CreateTexture(\$renderer, %PIXELFORMAT<RGB332>, STREAMING, \$w, \$h); my \$pixdatabuf = CArray[int64].new(0, \$w, \$h, \$w); my @color = 0x00, 0xDE, 0x14, 0xAA, 0xFF; sub render {    my int \$pitch;    my int \$cursor;     # work-around to pass the pointer-pointer.    my \$pixdata = nativecast(Pointer[int64], \$pixdatabuf);    SDL_LockTexture(\$asp_texture, SDL_Rect, \$pixdata, \$pitch);     \$pixdata = nativecast(CArray[int8], Pointer.new(\$pixdatabuf[0]));     loop (my int \$row; \$row < \$h; \$row = \$row + 1) {        my int \$rs = \$row * \$w; # row start        my int \$re = \$rs  + \$w; # row end        loop (my int \$idx = \$rs; \$idx < \$re; \$idx = \$idx + 1) {            \$pixdata[\$idx] =  @buffer[\$idx] < 4 ?? @color[@buffer[\$idx]] !! @color[4];            if @buffer[\$idx] >= 4 {                my \$grains = floor @buffer[\$idx] / 4;                @buffer[ \$idx - \$w ] += \$grains if \$row > 0;                @buffer[ \$idx - 1  ] += \$grains if \$idx - 1 >= \$rs;                @buffer[ \$idx + \$w ] += \$grains if \$row < \$h - 1;                @buffer[ \$idx + 1  ] += \$grains if \$idx + 1 < \$re;                @buffer[ \$idx ] %= 4;            }        }    }     SDL_UnlockTexture(\$asp_texture);     SDL_RenderCopy(\$renderer, \$asp_texture, SDL_Rect, SDL_Rect.new(:x(0), :y(0), :w(\$width), :h(\$height)));    SDL_RenderPresent(\$renderer);} my \$event = SDL_Event.new; main: loop {     while SDL_PollEvent(\$event) {        my \$casted_event = SDL_CastEvent(\$event);         given \$casted_event {            when *.type == QUIT {                last main;            }        }    }     render();    print fps;} say ''; sub fps {    state \$fps-frames = 0;    state \$fps-now    = now;    state \$fps        = '';    \$fps-frames++;    if now - \$fps-now >= 1 {        \$fps = [~] "\b" x 40, ' ' x 20, "\b" x 20 ,            sprintf "FPS: %5.2f  ", (\$fps-frames / (now - \$fps-now)).round(.01);        \$fps-frames = 0;        \$fps-now = now;    }    \$fps}`

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

Rust

`// This is the main algorithm.//// It loops over the current state of the sandpile and updates it on-the-fly.fn advance(field: &mut Vec<Vec<usize>>, boundary: &mut [usize; 4]) -> bool{     // This variable is used to check whether we changed anything in the array. If no, the loop terminates.    let mut done = false;     for y in boundary[0]..boundary[2]    {        for x in boundary[1]..boundary[3]        {            if field[y][x] >= 4            {                // This part was heavily inspired by the Pascal version. We subtract 4 as many times as we can                // and distribute it to the neighbors. Also, in case we have outgrown the current boundary, we                // update it to once again contain the entire sandpile.                 // The amount that gets added to the neighbors is the amount here divided by four and (implicitly) floored.                // The remaining sand is just current modulo 4.                let rem: usize = field[y][x] / 4;                field[y][x] %= 4;                 // The isize casts are necessary because usize can not go below 0.                // Also, the reason why x and y are compared to boundary[2]-1 and boundary[3]-1 is because for loops in                // Rust are upper bound exclusive. This means a loop like 0..5 will only go over 0,1,2,3 and 4.                if y as isize - 1 >= 0 {field[y-1][x] += rem; if y == boundary[0] {boundary[0]-=1;}}                if x as isize - 1 >= 0 {field[y][x-1] += rem; if x == boundary[1] {boundary[1]-=1;}}                if y+1 < field.len() {field[y+1][x] += rem; if x == boundary[2]-1 {boundary[2]+=1;}}                if x+1 < field.len() {field[y][x+1] += rem; if y == boundary[3]-1 {boundary[3]+=1;}}                 done = true;            }        }    }     done} // This function can be used to display the sandpile in the console window.//// Each row is mapped onto chars and those characters are then collected into a string.// These are then printed to the console.//// Eg.: [0,1,1,2,3,0] -> [' ','░','░','▒','▓',' ']-> " ░░▒▓ "fn display(field: &Vec<Vec<usize>>){    for row in field    {        let char_row = {            row.iter().map(|c| {match c {                0 => ' ',                1 => '░',                2 => '▒',                3 => '▓',                _ => '█'            }}).collect::<String>()        };        println!("{}", char_row);    }} // This function writes the end result to a file called "output.ppm".//// PPM is a very simple image format, however, it entirely uncompressed which leads to huge image sizes.// Even so, for demonstrative purposes it's perfectly good to use. For something more robust, look into PNG libraries.//// Read more about the format here: http://netpbm.sourceforge.net/doc/ppm.htmlfn write_pile(pile: &Vec<Vec<usize>>) {    use std::fs::File;    use std::io::Write;     // We first create the file (or erase its contents if it already existed).    let mut file = File::create("./output.ppm").unwrap();     // Then we add the image signature, which is "P3 <newline>[width of image] [height of image]<newline>[maximum value of color]<newline>".    write!(file, "P3\n{} {}\n255\n", pile.len(), pile.len()).unwrap();     for row in pile {        // For each row, we create a new string which has more or less enough capacity to hold the entire row.        // This is for performance purposes, but shouldn't really matter much.        let mut line = String::with_capacity(row.len() * 14);         // We map each value in the field to a color.        // These are just simple RGB values, 0 being the background, the rest being the "sand" itself.        for elem in row {            line.push_str(match elem {                0 => "100 40 15 ",                1 => "117 87 30 ",                2 => "181 134 47 ",                3 => "245 182 66 ",                _ => unreachable!(),            });        }         // Finally we write this string into the file.        write!(file, "{}\n", line).unwrap();    }} fn main() {    // This is how big the final image will be. Currently the end result would be a 16x16 picture.    let field_size = 16;    let mut playfield = vec![vec![0; field_size]; field_size];     // We put the initial sand in the exact middle of the field.    // This isn't necessary per se, but it ensures that sand can fully topple.    //    // The boundary is initially just the single tile which has the sand in it, however, as the algorithm    // progresses, this will grow larger too.    let mut boundary = [field_size/2-1, field_size/2-1, field_size/2, field_size/2];    playfield[field_size/2 - 1][field_size/2 - 1] = 16;     // This is the main loop. We update the field until it returns false, signalling that the pile reached its    // final state.    while advance(&mut playfield, &mut boundary) {};     // Once this happens, we simply display the result. Uncomment the line below to write it to a file.    // Calling display with large field sizes is not recommended as it can easily become too large for the console.    display(&playfield);    //write_pile(&playfield);}`

Output:

```

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

```

VBA

`Sub SetupPile(a As Integer, b As Integer)Application.ScreenUpdating = FalseFor i = 1 To aFor j = 1 To bCells(i, j).value = ""Cells(i, j).Select With Selection.Borders(xlEdgeLeft)    .LineStyle = xlContinuous    .Weight = xlMediumEnd WithWith Selection.Borders(xlEdgeTop)    .LineStyle = xlContinuous    .Weight = xlMediumEnd WithWith Selection.Borders(xlEdgeBottom)    .LineStyle = xlContinuous    .Weight = xlMediumEnd WithWith Selection.Borders(xlEdgeRight)    .LineStyle = xlContinuous    .Weight = xlMediumEnd With With Selection    .HorizontalAlignment = xlCenter    .VerticalAlignment = xlCenterEnd With Next jNext iApplication.ScreenUpdating = TrueEnd Sub  Sub Abelian_Sandpile()Dim PileWidth As IntegerDim PileHeight As IntegerDim FieldArray() As Integer Debug.Print "Start:" & Now() 'Set Size of Playing FieldPileWidth = 25PileHeight = 25 ReDim FieldArray(PileWidth - 1, PileHeight - 1) 'Paint Basic GridSetupPile PileWidth, PileHeight 'Drop sand amount into middle of playing fieldSandDropAmount = 1000'Get around excel's incorrect roundingSandDropColumn = Round((PileWidth / 2) + 0.001, 0)SandDropRow = Round((PileHeight / 2) + 0.001, 0) Cells(SandDropRow, SandDropColumn) = SandDropAmountFieldArray(SandDropRow - 1, SandDropColumn - 1) = SandDropAmount Continue = False 'Check if Pile is already stabilized at the startFor i = 1 To PileWidth 'ColFor j = 1 To PileHeight 'RowIf FieldArray(j - 1, i - 1) > 3 Then Continue = TrueNext jNext i 'While not stabilizedWhile ContinueFor i = 1 To PileWidthFor j = 1 To PileHeight    If FieldArray(j - 1, i - 1) > 3 Then    'Reduce by 4    FieldArray(j - 1, i - 1) = FieldArray(j - 1, i - 1) - 4    'Increase Neighbours    't    If j >= 2 Then FieldArray(j - 2, i - 1) = FieldArray(j - 2, i - 1) + 1    'r    If i < PileWidth Then FieldArray(j - 1, i) = FieldArray(j - 1, i) + 1    'b    If j < PileHeight Then FieldArray(j, i - 1) = FieldArray(j, i - 1) + 1    'l    If i >= 2 Then FieldArray(j - 1, i - 2) = FieldArray(j - 1, i - 2) + 1    'Next round    GoTo Nextone    End IfNext jNext i Nextone: 'Check if now stabilizedContinue = FalseFor i = 1 To PileWidthFor j = 1 To PileHeight'Paint every step if needed'Cells(j, i) = FieldArray(j - 1, i - 1) If FieldArray(j - 1, i - 1) > 3 Then Continue = TrueNext jNext i Wend 'Print out final stepFor i = 1 To PileWidthFor j = 1 To PileHeightCells(j, i) = FieldArray(j - 1, i - 1)Next jNext i 'Make field square and remove 0Cells.SelectSelection.ColumnWidth = 2Selection.RowHeight = 13.5Selection.Replace What:="0", Replacement:="", LookAt:=xlPart, SearchOrder:=xlByRows, MatchCase:=False, SearchFormat:=False, ReplaceFormat:=FalseRange("A1").Select Range(Cells(1, 1), Cells(PileHeight, PileWidth)).Select 'Conditional FormatSelection.FormatConditions.AddColorScale ColorScaleType:=3Selection.FormatConditions(Selection.FormatConditions.Count).SetFirstPrioritySelection.FormatConditions(1).ColorScaleCriteria(1).Type = xlConditionValueLowestValueWith Selection.FormatConditions(1).ColorScaleCriteria(1).FormatColor    .Color = 8109667    .TintAndShade = 0End WithSelection.FormatConditions(1).ColorScaleCriteria(2).Type = xlConditionValuePercentileSelection.FormatConditions(1).ColorScaleCriteria(2).value = 50With Selection.FormatConditions(1).ColorScaleCriteria(2).FormatColor    .Color = 8711167    .TintAndShade = 0End WithSelection.FormatConditions(1).ColorScaleCriteria(3).Type = xlConditionValueHighestValueWith Selection.FormatConditions(1).ColorScaleCriteria(3).FormatColor    .Color = 7039480    .TintAndShade = 0End WithRange("A1").Select Debug.Print "W,H,A:" & PileWidth & "," & PileHeight & "," & SandDropAmountDebug.Print "End:" & Now() End Sub`

Output:

```On Excel Page
```

Wren

Library: Wren-fmt
`import "/fmt" for Fmt class Sandpile {    // 'a' is a list of integers in row order    construct new(a) {        var count = a.count        _rows = count.sqrt.floor        if (_rows * _rows != count) Fiber.abort("The matrix of values must be square.")        _a = a        _neighbors = List.filled(count, 0)        for (i in 0...count) {            _neighbors[i] = []            if (i % _rows > 0)     _neighbors[i].add(i-1)            if ((i + 1)%_rows > 0) _neighbors[i].add(i+1)            if (i - _rows >= 0)    _neighbors[i].add(i-_rows)            if (i + _rows < count) _neighbors[i].add(i+_rows)        }    }     isStable { _a.all { |i| i <= 3 } }     // topples until stable    topple() {        while (!isStable) {            for (i in 0..._a.count) {                if (_a[i] > 3) {                    _a[i] = _a[i] - 4                    for (j in _neighbors[i]) _a[j] = _a[j] + 1                }            }        }    }     toString {        var s = ""        for (i in 0..._rows) {            for (j in 0..._rows) s = s + Fmt.swrite("\$2d ", _a[_rows*i + j])            s = s + "\n"        }        return s    }} var printAcross = Fn.new { |str1, str2|    var r1 = str1.split("\n")    var r2 = str2.split("\n")    var rows = r1.count - 1    var cr = (rows/2).floor    for (i in 0...rows) {        var symbol = (i == cr) ? "->" : "  "        Fmt.print("\$s \$s \$s", r1[i], symbol, r2[i])    }    System.print()} var a1 = List.filled(25, 0)a1[12] = 4var a2 = List.filled(25, 0)a2[12] = 6var a3 = List.filled(25, 0)a3[12] = 16var a4 = List.filled(100, 0)a4[55] = 64for (a in [a1, a2, a3, a4]) {    var s = Sandpile.new(a)    var str1 = s.toString    s.topple()    var str2 = s.toString    printAcross.call(str1, str2)}`
Output:
``` 0  0  0  0  0      0  0  0  0  0
0  0  0  0  0      0  0  1  0  0
0  0  4  0  0  ->  0  1  0  1  0
0  0  0  0  0      0  0  1  0  0
0  0  0  0  0      0  0  0  0  0

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

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

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