Abelian sandpile model/Identity

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

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

E.g. s1 =

    1 2 0
    2 1 1
    0 1 3

Or s2 =

    2 1 3
    1 0 1
    0 1 0

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

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

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

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

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

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

The final result is the stable sandpile on the right.

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

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

Show confirming output here, with your example.

References


C++

<lang cpp>#include <algorithm>

  1. include <array>
  2. include <cassert>
  3. include <initializer_list>
  4. include <iostream>

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

class abelian_sandpile {

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

public:

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

private:

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

};

abelian_sandpile::abelian_sandpile() {

   cells_.fill(0);

}

abelian_sandpile::abelian_sandpile(std::initializer_list<int> init) {

   assert(init.size() == sp_cells);
   std::copy(init.begin(), init.end(), cells_.begin());

}

abelian_sandpile& abelian_sandpile::operator+=(const abelian_sandpile& other) {

   for (size_t i = 0; i < sp_cells; ++i)
       cells_[i] += other.cells_[i];
   stabilize();
   return *this;

}

bool abelian_sandpile::operator==(const abelian_sandpile& other) {

   return cells_ == other.cells_;

}

bool abelian_sandpile::is_stable() const {

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

}

void abelian_sandpile::topple() {

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

}

void abelian_sandpile::stabilize() {

   while (!is_stable())
       topple();

}

abelian_sandpile operator+(const abelian_sandpile& a, const abelian_sandpile& b) {

   abelian_sandpile c(a);
   c += b;
   return c;

}

std::ostream& operator<<(std::ostream& out, const abelian_sandpile& as) {

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

}

int main() {

   std::cout << std::boolalpha;
   std::cout << "Avalanche:\n";
   abelian_sandpile sp{4,3,3, 3,1,2, 0,2,3};
   while (!sp.is_stable()) {
       std::cout << sp << "stable? " << sp.is_stable() << "\n\n";
       sp.topple();
   }
   std::cout << sp << "stable? " << sp.is_stable() << "\n\n";
   std::cout << "Commutativity:\n";
   abelian_sandpile s1{1,2,0, 2,1,1, 0,1,3};
   abelian_sandpile s2{2,1,3, 1,0,1, 0,1,0};
   abelian_sandpile sum1(s1 + s2);
   abelian_sandpile sum2(s2 + s1);
   std::cout << "s1 + s2 equals s2 + s1? " << (sum1 == sum2) << "\n\n";
   std::cout << "s1 + s2 = \n" << sum1;
   std::cout << "\ns2 + s1 = \n" << sum2;
   std::cout << '\n';
   std::cout << "Identity:\n";
   abelian_sandpile s3{3,3,3, 3,3,3, 3,3,3};
   abelian_sandpile s3_id{2,1,2, 1,0,1, 2,1,2};
   abelian_sandpile sum3(s3 + s3_id);
   abelian_sandpile sum4(s3_id + s3_id);
   std::cout << "s3 + s3_id equals s3? " << (sum3 == s3) << '\n';
   std::cout << "s3_id + s3_id equals s3_id? " << (sum4 == s3_id) << "\n\n";
   std::cout << "s3 + s3_id = \n" << sum3;
   std::cout << "\ns3_id + s3_id = \n" << sum4;
   return 0;

}</lang>

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

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

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

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

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

Commutativity:
s1 + s2 equals s2 + s1? true

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

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

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

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

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

F#

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

Output:
[[3; 3; 3]
 [3; 1; 2]
 [0; 2; 3]]

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

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

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

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

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

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

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

Factor

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

Works with: Factor version 0.99 2020-07-03

<lang factor>USING: arrays grouping io kernel math math.vectors prettyprint qw sequences ;

CONSTANT: neighbors {

   { 1 3 } { 0 2 4 } { 1 5 } { 0 4 6 } { 1 3 5 7 }
   { 2 4 8 } { 3 7 } { 4 6 8 } { 5 7 }

}

! Sandpile words

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

! Output words

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

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

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

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

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

"s3_id + s3_id = s3_id" print nl id id .s+</lang>

Output:
Avalanche:

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

s1 + s2 = s2 + s1

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

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

s3 + s3_id = s3

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

s3_id + s3_id = s3_id

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

Go

Translation of: Wren

<lang go>package main

import (

   "fmt"
   "strconv"
   "strings"

)

type sandpile struct{ a [9]int }

var neighbors = [][]int{

   {1, 3}, {0, 2, 4}, {1, 5}, {0, 4, 6}, {1, 3, 5, 7}, {2, 4, 8}, {3, 7}, {4, 6, 8}, {5, 7},

}

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

func (s *sandpile) plus(other *sandpile) *sandpile {

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

}

func (s *sandpile) isStable() bool {

   for _, e := range s.a {
       if e > 3 {
           return false
       }
   }
   return true

}

// just topples once so we can observe intermediate results func (s *sandpile) topple() {

   for i := 0; i < 9; i++ {
       if s.a[i] > 3 {
           s.a[i] -= 4
           for _, j := range neighbors[i] {
               s.a[j]++
           }
           return
       }
   }

}

func (s *sandpile) String() string {

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

}

func main() {

   fmt.Println("Avalanche of topplings:\n")
   s4 := newSandpile([9]int{4, 3, 3, 3, 1, 2, 0, 2, 3})
   fmt.Println(s4)
   for !s4.isStable() {
       s4.topple()
       fmt.Println(s4)
   }
   fmt.Println("Commutative additions:\n")
   s1 := newSandpile([9]int{1, 2, 0, 2, 1, 1, 0, 1, 3})
   s2 := newSandpile([9]int{2, 1, 3, 1, 0, 1, 0, 1, 0})
   s3_a := s1.plus(s2)
   for !s3_a.isStable() {
       s3_a.topple()
   }
   s3_b := s2.plus(s1)
   for !s3_b.isStable() {
       s3_b.topple()
   }
   fmt.Printf("%s\nplus\n\n%s\nequals\n\n%s\n", s1, s2, s3_a)
   fmt.Printf("and\n\n%s\nplus\n\n%s\nalso equals\n\n%s\n", s2, s1, s3_b)
   fmt.Println("Addition of identity sandpile:\n")
   s3 := newSandpile([9]int{3, 3, 3, 3, 3, 3, 3, 3, 3})
   s3_id := newSandpile([9]int{2, 1, 2, 1, 0, 1, 2, 1, 2})
   s4 = s3.plus(s3_id)
   for !s4.isStable() {
       s4.topple()
   }
   fmt.Printf("%s\nplus\n\n%s\nequals\n\n%s\n", s3, s3_id, s4)
   fmt.Println("Addition of identities:\n")
   s5 := s3_id.plus(s3_id)
   for !s5.isStable() {
       s5.topple()
   }
   fmt.Printf("%s\nplus\n\n%s\nequals\n\n%s", s3_id, s3_id, s5)

}</lang>

Output:
Avalanche of topplings:

4 3 3 
3 1 2 
0 2 3 

0 4 3 
4 1 2 
0 2 3 

1 0 4 
4 2 2 
0 2 3 

1 1 0 
4 2 3 
0 2 3 

2 1 0 
0 3 3 
1 2 3 

Commutative additions:

1 2 0 
2 1 1 
0 1 3 

plus

2 1 3 
1 0 1 
0 1 0 

equals

3 3 3 
3 1 2 
0 2 3 

and

2 1 3 
1 0 1 
0 1 0 

plus

1 2 0 
2 1 1 
0 1 3 

also equals

3 3 3 
3 1 2 
0 2 3 

Addition of identity sandpile:

3 3 3 
3 3 3 
3 3 3 

plus

2 1 2 
1 0 1 
2 1 2 

equals

3 3 3 
3 3 3 
3 3 3 

Addition of identities:

2 1 2 
1 0 1 
2 1 2 

plus

2 1 2 
1 0 1 
2 1 2 

equals

2 1 2 
1 0 1 
2 1 2 

Haskell

<lang haskell>{-# LANGUAGE TupleSections #-}

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


TEST ---------------------------

main :: IO () main = do

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

SAND PILES ------------------------

addSand :: Int -> Int -> Int addSand xs ys =

 (head . cascade . chunksOf (length xs)) $ zipWith (+) (concat xs) (concat ys)

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

 where
   w = length xs

convergence :: (a -> a -> Bool) -> [a] -> [a] convergence p = go

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

tumble :: Int -> [Int] -> [Int] tumble w xs = maybe xs go $ findIndex (w <) xs

 where
   go i = zipWith f [0 ..] xs
     where
       neighbours = indexNeighbours w i
       f j x
         | j `elem` neighbours = succ x
         | i == j = x - succ w
         | otherwise = x

indexNeighbours :: Int -> Int -> [Int] indexNeighbours w = go

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

DISPLAY --------------------------

showCascade :: [(String, Int)] -> String showCascade pairs =

 unlines $
 fmap unwords $
 transpose $
 fmap
   (\(pfx, xs) ->
       unwords <$> transpose (centered pfx : transpose (fmap (fmap show) xs)))
   pairs

centered :: String -> [String] centered s = [pad, s, pad <> replicate r ' ']

 where
   lng = length s
   pad = replicate lng ' '
   (q, r) = quotRem (2 + lng) 2</lang>
Cascade:
 4 3 3    0 4 3    1 0 4    1 1 0    2 1 0
 3 1 2 -> 4 1 2 -> 4 2 2 -> 4 2 3 -> 0 3 3
 0 2 3    0 2 3    0 2 3    0 2 3    1 2 3

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

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

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

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

J

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

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

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

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

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

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

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

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

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

   raku_id verify_identity raku
1

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

Julia

<lang julia>import Base.+, Base.print

struct Sandpile

   pile::Matrix{UInt8}

end

function Sandpile(s::String)

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

end

const HMAX = 3

function avalanche!(s::Sandpile, lim=HMAX)

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

end

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

function print(io::IO, s::Sandpile)

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

end

const s1 = Sandpile("""

   1 2 0
   2 1 1
   0 1 3""")

const s2 = Sandpile("""

   2 1 3
   1 0 1
   0 1 0""")

const s3 = Sandpile("""

   3 3 3
   3 3 3
   3 3 3""")

const s3_id = Sandpile("""

   2 1 2
   1 0 1
   2 1 2""")

const s3a = Sandpile("""

   4 3 3
   3 1 2
   0 2 3""")

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

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

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

</lang>

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


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

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


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


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

Lua

Uses Abelian sandpile model here, then extends.. <lang Lua>sandpile.__index = sandpile sandpile.new = function(self, vals)

 local inst = setmetatable({},sandpile)
 inst.cell, inst.dim = {}, #vals
 for r = 1, inst.dim do
   inst.cell[r] = {}
   for c = 1, inst.dim do
     inst.cell[r][c] = vals[r][c]
   end
 end
 return inst

end sandpile.add = function(self, other)

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

end

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

Output:
s1 =
1 2 0
2 1 1
0 1 3

s2 =
2 1 3
1 0 1
0 1 0

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

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

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

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

s3 =
3 3 3
3 3 3
3 3 3

s3_id =
2 1 2
1 0 1
2 1 2

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

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

Nim

<lang Nim> import sequtils import strutils

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

  1. ---------------------------------------------------------------------------------------------------

iterator neighbors(i, j: int): tuple[a, b: int] =

 ## Yield the indexes of the neighbours of cell at indexes (i, j).
 if i > 0:
   yield (i - 1, j)
 if i < 2:
   yield (i + 1, j)
 if j > 0:
   yield (i, j - 1)
 if j < 2:
   yield (i, j + 1)
  1. ---------------------------------------------------------------------------------------------------

proc print(s: openArray[SandPile]) =

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

proc printSum(s1, s2, s3: SandPile) =

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

func isStable(sandPile: SandPile): bool =

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

proc topple(sandPile: var SandPile) =

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

proc stabilize(sandPile: var SandPile) =

 ## Stabilize a sandpile.
 while not sandPile.isStable():
   sandPile.topple()
  1. ---------------------------------------------------------------------------------------------------

proc `+`(s1, s2: SandPile): SandPile =

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

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

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

 s.topple()
 list.add(s)

list.print() echo Separator

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

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

echo "s3_id + s3_id = s3_id\n" printSum(s3_id, s3_id, s3_id + s3_id) </lang>

Output:
Avalanche

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

-----

s1 + s2 == s2 + s1

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

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

-----

s3 + s3_id == s3

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

-----

s3_id + s3_id = s3_id

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

Phix

<lang Phix>constant s1 = {"1 2 0",

              "2 1 1",
              "0 1 3"},
        s2 = {"2 1 3",
              "1 0 1",
              "0 1 0"},
        s3 = {"3 3 3",
              "3 3 3",
              "3 3 3"},
     s3_id = {"2 1 2",
              "1 0 1",
              "2 1 2"},
        s4 = {"4 3 3",
              "3 1 2",
              "0 2 3"}

function add(sequence s, t)

   for i=1 to 3 do
       for j=1 to 5 by 2 do
           s[i][j] += t[i][j]-'0'
       end for
   end for
   return s

end function

function topple(sequence s, integer one=0)

   for i=1 to 3 do
       for j=1 to 5 by 2 do
           if s[i][j]>'3' then
               s[i][j] -= 4
               if i>1 then s[i-1][j] += 1 end if
               if i<3 then s[i+1][j] += 1 end if
               if j>1 then s[i][j-2] += 1 end if
               if j<5 then s[i][j+2] += 1 end if
               if one=1 then return s end if
               one = -1
           end if
       end for
   end for
   return iff(one=1?{}:iff(one=-1?topple(s):s))

end function

procedure shout(sequence s)

   sequence r = repeat("",5)
   for i=1 to length(s) do
       sequence si = s[i]
       if string(si) then
           string ti = repeat(' ',length(si))
           r[1] &= ti
           r[2] &= si
           r[3] &= ti
       else
           for j=1 to 3 do
               r[j] &= si[j]
           end for
       end if
   end for
   puts(1,join(r,"\n"))

end procedure

puts(1,"1. Show avalanche\n\n") sequence s = s4,

        res = {"    ",s}

while true do

   s = topple(s,1)
   if s={} then exit end if
   res &= {" ==> ",s}

end while shout(res)

puts(1,"2. Prove s1 + s2 = s2 + s1\n\n") shout({" ",s1," + ",s2," = ",topple(add(s1,s2))}) shout({" ",s2," + ",s1," = ",topple(add(s2,s1))})

puts(1,"3. Show that s3 + s3_id == s3\n\n") shout({" ",s3," + ",s3_id," = ",topple(add(s3,s3_id))})

puts(1,"4. Show that s3_id + s3_id == s3_id\n\n") shout({" ",s3_id," + ",s3_id," = ",topple(add(s3_id,s3_id))})</lang>

Output:
1. Show avalanche

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

2. Prove s1 + s2 = s2 + s1

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

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

3. Show that s3 + s3_id == s3

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

4. Show that s3_id + s3_id == s3_id

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

Python

Object Oriented

<lang python>from itertools import product from collections import defaultdict


class Sandpile():

   def __init__(self, gridtext):
       array = [int(x) for x in gridtext.strip().split()]
       self.grid = defaultdict(int,
                               {(i //3, i % 3): x 
                                for i, x in enumerate(array)})
   _border = set((r, c) 
                 for r, c in product(range(-1, 4), repeat=2) 
                 if not 0 <= r <= 2 or not 0 <= c <= 2
                 )
   _cell_coords = list(product(range(3), repeat=2))
   
   def topple(self):
       g = self.grid
       for r, c in self._cell_coords:
           if g[(r, c)] >= 4:
               g[(r - 1, c)] += 1
               g[(r + 1, c)] += 1
               g[(r, c - 1)] += 1
               g[(r, c + 1)] += 1
               g[(r, c)] -= 4
               return True
       return False
   
   def stabilise(self):
       while self.topple():
           pass
       # Remove extraneous grid border
       g = self.grid
       for row_col in self._border.intersection(g.keys()):
           del g[row_col]
       return self
   
   __pos__ = stabilise     # +s == s.stabilise()
   
   def __eq__(self, other):
       g = self.grid
       return all(g[row_col] == other.grid[row_col]
                  for row_col in self._cell_coords)
   def __add__(self, other):
       g = self.grid
       ans = Sandpile("")
       for row_col in self._cell_coords:
           ans.grid[row_col] = g[row_col] + other.grid[row_col]
       return ans.stabilise()
      
   def __str__(self):
       g, txt = self.grid, []
       for row in range(3):
           txt.append(' '.join(str(g[(row, col)]) 
                               for col in range(3)))
       return '\n'.join(txt)
   
   def __repr__(self):
       return f'{self.__class__.__name__}(""""\n{self.__str__()}""")'
       

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

   1 2 0
   2 1 1
   0 1 3

""") s2 = Sandpile("""

   2 1 3
   1 0 1
   0 1 0 

""") s3 = Sandpile("3 3 3 3 3 3 3 3 3") s3_id = Sandpile("2 1 2 1 0 1 2 1 2") </lang>

Command line session to complete task.
In [2]: unstable
Out[2]: 
Sandpile(""""
4 3 3
3 1 2
0 2 3""")

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

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

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

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

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

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

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

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

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

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

In [13]: 

Functional

<lang python>Abelian Sandpile – Identity

from operator import add, eq


  1. -------------------------- TEST --------------------------
  2. main :: IO ()

def main():

   Tests of cascades and additions
   s0 = [[4, 3, 3], [3, 1, 2], [0, 2, 3]]
   s1 = [[1, 2, 0], [2, 1, 1], [0, 1, 3]]
   s2 = [[2, 1, 3], [1, 0, 1], [0, 1, 0]]
   s3 = [[3, 3, 3], [3, 3, 3], [3, 3, 3]]
   s3_id = [[2, 1, 2], [1, 0, 1], [2, 1, 2]]
   series = list(cascadeSeries(s0))
   for expr in [
           'Cascade:',
           showSandPiles(
               [(' ', series[0])] + [
                   (':', xs) for xs in series[1:]
               ]
           ),
           ,
           f's1 + s2 == s2 + s1 -> {addSand(s1)(s2) == addSand(s2)(s1)}',
           showSandPiles([
               (' ', s1),
               ('+', s2),
               ('=', addSand(s1)(s2))
           ]),
           ,
           showSandPiles([
               (' ', s2),
               ('+', s1),
               ('=', addSand(s2)(s1))
           ]),
           ,
           f's3 + s3_id == s3 -> {addSand(s3)(s3_id) == s3}',
           showSandPiles([
               (' ', s3),
               ('+', s3_id),
               ('=', addSand(s3)(s3_id))
           ]),
           ,
           f's3_id + s3_id == s3_id -> {addSand(s3_id)(s3_id) == s3_id}',
           showSandPiles([
               (' ', s3_id),
               ('+', s3_id),
               ('=', addSand(s3_id)(s3_id))
           ]),
   ]:
       print(expr)


  1. ----------------------- SANDPILES ------------------------
  1. addSand :: Int -> Int -> Int

def addSand(xs):

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


  1. cascadeSeries :: Int -> [[[Int]]]

def cascadeSeries(rows):

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


  1. convergence :: (a -> a -> Bool) -> [a] -> [a]

def convergence(p):

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


  1. nextState Int -> Int -> [Int] -> [Int]

def nextState(w):

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


  1. indexNeighbours :: Int -> Int -> [Int]

def indexNeighbours(w):

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


  1. ------------------------ DISPLAY -------------------------
  1. showSandPiles :: [(String, Int)] -> String

def showSandPiles(pairs):

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


  1. ------------------------ GENERIC -------------------------
  1. chunksOf :: Int -> [a] -> a

def chunksOf(n):

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


  1. concat :: a -> [a]

def concat(xs):

   The concatenation of all
      elements in a list.
   
   return [x for lst in xs for x in lst]


  1. findIndex :: (a -> Bool) -> [a] -> Maybe Int

def findIndex(p):

   Just the first index at which an
      element in xs matches p,
      or Nothing if no elements match.
   
   def go(xs):
       return next(
           (i for (i, x) in enumerate(xs) if p(x)),
           None
       )
   return go


  1. iterate :: (a -> a) -> a -> Gen [a]

def iterate(f):

   An infinite list of repeated
      applications of f to x.
   
   def go(x):
       v = x
       while True:
           yield v
           v = f(v)
   return go


  1. maybe :: b -> (a -> b) -> Maybe a -> b

def maybe(v):

   Either the default value v, if x is None,
      or the application of f to x.
   
   def go(f):
       def g(x):
           return v if None is x else f(x)
       return g
   return go


  1. MAIN ---

if __name__ == '__main__':

   main()</lang>
Output:
Cascade:
  4 3 3   0 4 3   1 0 4   1 1 0   2 1 0
  3 1 2 : 4 1 2 : 4 2 2 : 4 2 3 : 0 3 3
  0 2 3   0 2 3   0 2 3   0 2 3   1 2 3

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

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

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

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

Raku

Works with: Rakudo version 2020.05

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

<lang perl6>class ASP {

   has $.h = 3;
   has $.w = 3;
   has @.pile = 0 xx $!w * $!h;
   method topple {
       my $buf = $!w * $!h;
       my $done;
       repeat {
           $done = True;
           loop (my int $row; $row < $!h; $row = $row + 1) {
               my int $rs = $row * $!w; # row start
               my int $re = $rs  + $!w; # row end
               loop (my int $idx = $rs; $idx < $re; $idx = $idx + 1) {
                   if self.pile[$idx] >= 4 {
                       my $grain = self.pile[$idx] div 4;
                       self.pile[ $idx - $!w ] += $grain if $row > 0;
                       self.pile[ $idx - 1  ]  += $grain if $idx - 1 >= $rs;
                       self.pile[ $idx + $!w ] += $grain if $row < $!h - 1;
                       self.pile[ $idx + 1  ]  += $grain if $idx + 1 < $re;
                       self.pile[ $idx ] %= 4;
                       $done = False;
                   }
               }
           }
       } until $done;
       self.pile;
   }

}

  1. some handy display layout modules

use Terminal::Boxer:ver<0.2+>; use Text::Center;

for 3, (4,3,3,3,1,2,0,2,3), (2,1,2,1,0,1,2,1,2), # 3 square task example

   3, (2,1,2,1,0,1,2,1,2), (2,1,2,1,0,1,2,1,2), # 3 square identity
   5, (4,1,0,5,1,9,3,6,1,0,8,1,2,5,3,3,0,1,7,5,4,2,2,4,0), (2,3,2,3,2,3,2,1,2,3,2,1,0,1,2,3,2,1,2,3,2,3,2,3,2) # 5 square test
 -> $size, $pile, $identity {
   my $asp = ASP.new(:h($size), :w($size));
   $asp.pile = |$pile;
   my @display;
   my %p = :col($size), :3cw, :indent("\t");
   @display.push: rs-box |%p, |$identity;
   @display.push: rs-box |%p, $asp.pile;
   @display.push: rs-box |%p, $asp.topple;
   $asp.pile Z+= $identity.list;
   @display.push: rs-box |%p, $asp.pile;
   @display.push: rs-box |%p, $asp.topple;
   put %p<indent> ~ qww<identity 'test pile' toppled 'plus identity' toppled>».&center($size * 4 + 1).join: %p<indent>;
   .put for [Z~] @display».lines;
   put ;

}</lang>

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

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

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

Wren

Library: Wren-fmt

<lang ecmascript>import "/fmt" for Fmt

class Sandpile {

   static init() {
       __neighbors = [
           [1, 3], [0, 2, 4], [1, 5], [0, 4, 6], [1, 3, 5, 7], [2, 4, 8], [3, 7], [4, 6, 8], [5, 7]
       ]
   }
   // 'a' is a list of 9 integers in row order
   construct new(a) {
       _a = a
   }
   a { _a }
   +(other) {
       var b = List.filled(9, 0)
       for (i in 0..8) b[i] = _a[i] + other.a[i]
       return Sandpile.new(b)
   }
   isStable { _a.all { |i| i <= 3 } }
   // just topples once so we can observe intermediate results
   topple() {
       for (i in 0..8) {
           if (_a[i] > 3) {
               _a[i] = _a[i] - 4
               for (j in __neighbors[i]) _a[j] = _a[j] + 1
               return
           }
       }
   }
   toString {
       var s = ""
       for (i in 0..2) {
           for (j in 0..2) s = s + "%(a[3*i + j]) "
           s = s + "\n"
       }
       return s
   }

}

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

   s4.topple()
   System.print(s4)

}

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

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

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

Output:
Avalanche of topplings:

4 3 3 
3 1 2 
0 2 3 

0 4 3 
4 1 2 
0 2 3 

1 0 4 
4 2 2 
0 2 3 

1 1 0 
4 2 3 
0 2 3 

2 1 0 
0 3 3 
1 2 3 

Commutative additions:

1 2 0 
2 1 1 
0 1 3 

plus

2 1 3 
1 0 1 
0 1 0 

equals

3 3 3 
3 1 2 
0 2 3 

and

2 1 3 
1 0 1 
0 1 0 

plus

1 2 0 
2 1 1 
0 1 3 

also equals

3 3 3 
3 1 2 
0 2 3 

Addition of identity sandpile:

3 3 3 
3 3 3 
3 3 3 

plus

2 1 2 
1 0 1 
2 1 2 

equals

3 3 3 
3 3 3 
3 3 3 

Addition of identities:

2 1 2 
1 0 1 
2 1 2 

plus

2 1 2 
1 0 1 
2 1 2 

equals

2 1 2 
1 0 1 
2 1 2 

Rust

<lang Rust>#[derive(Clone)] struct Box {

   piles: [[u8; 3]; 3],

}

impl Box {

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

}

fn main() {

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

} </lang>

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