# 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. It's history, mathematical definition and properties can be found under it's 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 it's 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
```

## 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; 	if(argc!=3){		printf("Usage: %s <Sand pile side> <Center pile height>",argv);		return 0;	} 	sandPileEdge = atoi(argv);	centerPileHeight = atoi(argv); 	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) + strlen(argv) + 23)*sizeof(char)); 	strcpy(fileName,"Final_Sand_Pile_");	strcat(fileName,argv);	strcat(fileName,"_");	strcat(fileName,argv);	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 = (sandPile[i][j] + i)%256;			colour = (sandPile[i][j] + j)%256;			colour = (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")

target_compile_options(abelian_sandpile PRIVATE -march=native -std=c++14)

```

## Fōrmulæ

Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text (more info). 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 transportation effects more than visualization and edition.

The option to show Fōrmulæ programs and their results is showing images. Unfortunately images cannot be uploaded in Rosetta Code.

## 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: 
sandpile with 50000 grains of sand: ./sandpile.fs 201 50000: 
sandpile with 500000 grains of sand: ./sandpile.fs 601 500000: 

## 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: 
sandpile with 10000 grains of sand: test 40 10000: 
sandpile with 100000 grains of sand: test 150 100000: 
sandpile with 1000000 grains of sand: test 400 1000000: 

## 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:

## 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,length(G)*SizeOf(G),#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;  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;  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;    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
```

## 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  }`

## Perl 6

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 {                ++@buffer[ \$idx - \$w ] if \$row > 0;                ++@buffer[ \$idx - 1  ] if \$idx - 1 > \$rs;                ++@buffer[ \$idx + \$w ] if \$row < \$h - 1;                ++@buffer[ \$idx + 1  ] if \$idx + 1 < \$buf;                @buffer[ \$idx ] -= 4;                \$done = False;            }        }    }    unless \$hide-progress {        print join '', @buffer.map( { @color[\$_ min 4] })    }} until \$done; print join '', @buffer.map( { @color[\$_ min 4] }); cleanup;`

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

## Phix

Library: 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.exwinclude pGUI.e Ihandle dlg, canvascdCanvas 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] += 1            end for        elsif x=length(board) 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            end for        elsif y=length(board) then              -- extend down            board = append(board,board)        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)/2)+1    drop(y,x)    return IUP_DEFAULTend 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) 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_DEFAULTend 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_DEFAULTend 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\"")    IupCloseOnEscape(dlg)    IupShow(dlg)    Ihandle timer = IupTimer(Icallback("timer_cb"), 10)    IupMainLoop()    IupClose()end procedure main()`

## Python

` 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.]] `

## Rust

`// Set image size.const DIM: usize = 16; // This function outputs the result to the console using UTF-8 block characters.fn draw_pile(pile: &Vec<Vec<usize>>) {    for row in pile {        let mut line = String::with_capacity(row.len());        for elem in row {            line.push(match elem {                0 => ' ',                1 => '░',                2 => '▒',                3 => '▓',                _ => '█'            });        }         println!("{}", line);    }} // This function creates a file called "output.ppm" in the directory the program was run, which contains// a colored image of the pile.fn write_pile(pile: &Vec<Vec<usize>>) {    use std::fs::File; // Used for opening the file.    use std::io::Write; // Used for writing to the file.     // Learn more about PPM here: http://netpbm.sourceforge.net/doc/ppm.html    let mut file = File::create("./output.ppm").unwrap();     // We write the signature, image dimensions and maximum color value to the file.    let _ = write!(file, "P3\n {} {}\n255\n", DIM, DIM).unwrap();     for row in pile {        let mut line = String::with_capacity(row.len()*6);        for elem in row {            line.push_str(match elem {                0 => "125 0 25 ", // Background color for cells that have no "sand" in them.                 // Depending on how many particles of sand is there in the cell we use a different shade of yellow.                1 => "125 80 0 ",                2 => "186 118 0 ",                3 => "224 142 0 ",                 // It is impossible to have more than 3 particles of sand in one cell after the program has run,                // however, Rust demands that all branches have to be considered in a match statement, so we                // explicitly tell the compiler, that this is an unreachable branch.                _ => unreachable!()             });        }         let _ = write!(file, "{}", line).unwrap();    }} // This is the main part of the algorithm, a simple, recursive implementation of the model.fn handle_pile(x: usize, y: usize, pile: &mut Vec<Vec<usize>>) {    if pile[y][x] >= 4 {        pile[y][x] -= 4;         // We check each neighbor, whether they have enough "sand" to collapse and if they do,         // we recursively call handle_pile on them.        if y > 0 {            pile[y-1][x] += 1;            if pile[y-1][x] >= 4 {handle_pile(x, y-1, pile)}}         if x > 0 {            pile[y][x-1] += 1;            if pile[y][x-1] >= 4 {handle_pile(x-1, y, pile)}}         if y < DIM-1 {            pile[y+1][x] += 1;            if pile[y+1][x] >= 4 {handle_pile(x, y+1, pile)}}         if x < DIM-1 {            pile[y][x+1] += 1;            if pile[y][x+1] >= 4 {handle_pile(x+1, y, pile)}}         // Uncomment this line to show every iteration of the program. Not recommended with large input values.        //draw_pile(&pile);         // Finally we call the function on the current cell again, in case it had more than 4 particles.        handle_pile(x,y,pile);    }}  fn main() {    use std::thread::Builder; // Used to spawn a new thread.     /* Rust by default uses a 2Mb stack, which gets quickly filled (resulting in a stack overflow) if we use any value larger than     * about 30,000 as our input value. To circumvent this, we spawn a thread with 32Mbs of stack memory, which can easily handle     * hundreds of thousands of sand particles. I tested the program using 256,000, but it should theoretically work with larger     * values too.      */     let _ = Builder::new().stack_size(33554432).spawn(|| {        // This is our 2D grid. It's size can be set using the DIM constant found at the top of the code.        let mut pile: Vec<Vec<usize>> = vec![vec![0;DIM]; DIM];         // We place this much sand in the center of the grid.        pile[DIM/2 - 1][DIM/2 - 1] = 16;         // We start the algorithm on the pile we just created.        handle_pile(DIM/2 - 1, DIM/2 - 1, &mut pile);          draw_pile(&pile)         // Uncomment this to save the image to a file after the recursive algorithm has ended.        //write_pile(&pile)    }).unwrap().join();}`

Output:

```

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

```