# Quicksand (piles)

In this challenge, you take a positive integer as input, which represents the height of a sand pile, located at (0,0) on an infinite square grid. For example, if our input is 123, the sand grid looks initially like this:

$$\\begin{matrix} \ddots & \vdots & \vdots & \vdots & \cdot^{\cdot^\cdot} \\ \cdots & 0 & 0 & 0 & \cdots\\ \cdots & 0 & 123 & 0 & \cdots \\ \cdots & 0 & 0 & 0 & \cdots \\ \cdot^{\cdot^\cdot} & \vdots & \vdots & \vdots & \ddots \end{matrix}\$$

Now, piles of sand are unstable, so they topple if their height is 4 or greater. When sand piles topple, they send sand equally in all four directions, but due to the quantum nature of sand, the amount sent is an integer. Or in other words, if there is a pile of height n, it sends sand as follows:

$$\\begin{matrix} \ddots & \vdots & \vdots & \vdots & \vdots & \vdots & \cdot^{\cdot^\cdot} \\ \cdots & 0 & 0 & 0 & 0 & 0 &\cdots\\ \cdots & 0 & 0 & \lfloor\frac{n}{4}\rfloor & 0 & 0 & \cdots\\ \cdots & 0 & \lfloor\frac{n}{4}\rfloor & n\mod 4 & \lfloor\frac{n}{4}\rfloor & 0 & \cdots\\ \cdots & 0 & 0 & \lfloor\frac{n}{4}\rfloor & 0 & 0 & \cdots\\ \cdots & 0 & 0 & 0 & 0 & 0 &\cdots\\ \cdot^{\cdot^\cdot} & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{matrix}\$$

This means that sand piles evolve step by step. For example, if we have this initial position, here's how the sand piles evolve (I've removed the dots for clarity).

$$\\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 6 & 0 \\ 0 & 2 & 7 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}\$$

$$\\begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 1 & 3 & 1 \\ 0 & 3 & 4 & 1 \\ 0 & 0 & 1 & 0 \end{matrix}\$$

$$\\begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 4 & 0 & 2 \\ 0 & 0 & 2 & 0 \end{matrix}\$$

$$\\begin{matrix} 0 & 0 & 2 & 0 \\ 0 & 3 & 0 & 2 \\ 1 & 0 & 2 & 2 \\ 0 & 1 & 2 & 0 \end{matrix}\$$

Now, since all the sand piles have less than 4 sand, this sand grid is stable.

# Rules

Your task is to take an positive integer $$\i\$$ as input, and output the eventually stable grid that the following initial position evolves to:

$$\\begin{matrix} \ddots & \vdots & \vdots & \vdots & \cdot^{\cdot^\cdot} \\ \cdots & 0 & 0 & 0 & \cdots\\ \cdots & 0 & i & 0 & \cdots \\ \cdots & 0 & 0 & 0 & \cdots \\ \cdot^{\cdot^\cdot} & \vdots & \vdots & \vdots & \ddots \end{matrix}\$$

I will execute your code with the following command:

time ./your_executable i > /dev/null

The highest i where the time command reports a real time of <=1 min will be your score. I will run the programs on an AMD Ryzen 7 1800X linux system with 32 GB of ram. Your code should also compile in 10 minutes or less, and should be less than 60kb. This is simply to prevent calculating the output at comptime or including the output in the source code.

Your program will choose some odd integer l and write to stdout l*l numbers ("0","1","2","3" meaning bytes 48, 49, 50 and 51) and nothing else (not even a trailing newline). These numbers represent the final state of the sand pile. l must be large enough to include the whole pile.

Your code doesn't have to work for all inputs i (for example, you can only accept powers of two, or just a single specific number). Please explicitly state if this is the case in your answer. Also, if your code works significantly faster for some inputs, please also tell this in the answer, as I can't test all values of i and instead have to use a binary search or something similar.

# Fractals

If you plot the outputs, you get really amazing fractals. Here's an example for a sand pile of height $$\3\times10^7\$$:

Credit colt_browning CC BY 4.0

What's amazing to me is that even though the rules are very simple and discrete, complex and organic behavior emerges both at a small and large scale.

Feel free to include a picture of your sand piles.

These are called Abelian sandpiles. Numberphile has a video about them.

Score Language Author
10000000 C (clang -march=native -Ofast) Polichinelle
8388608 Bash + xz + node (compressed) Community
4800000 Rust (rustc -Copt-level=3 -Ctarget-cpu=native) alephalpha
1572864 Python + numpy ovs
1376256 C (clang -march=native -Ofast) AnttiP
1350000 Rust (RUSTFLAGS="-Ctarget-cpu=native" cargo ...) Aiden4
410000 C (clang -march=native -O3) astroide
• Amazing challenge!! Commented Feb 6, 2022 at 14:55
• "What's amazing to me is that even though the rules are very simple and discrete, complex and organic behavior emerges both at a small and large scale." - this is because this a form of a cellular automaton, which are famous for their simple rules and complex results (e.g. Game of Life) Commented Feb 6, 2022 at 14:56
• This is not gonna win any prices for speed, but here is some animated sand dropping
– ovs
Commented Feb 6, 2022 at 15:14
• Why can't I output a square grid?
– Neil
Commented Feb 6, 2022 at 15:57
• Wesley Pegden (mentioned in the Numberphile video) has a nice gallery of zoomable images with up to $2^{30}$ grains of sand
– ovs
Commented Feb 7, 2022 at 15:09

# C (clang), 8895 bytes, score 14000000 on TIO (revised again)

This is a multithreaded version of my earlier revision, which is still available below. Please note:

• The number of threads, which is set by -DNUM_THREADS=..., should be a power of two. I used two threads for the run on TIO. The default is eight.
• The code uses spin locks. You can avoid this by compiling with the option -DUSE_COND_VAR, but the program becomes slower.
• The program uses C11 threads, so you may have to compile with the library -lpthread.
• For the run on TIO, I compiled with the flags -Ofast -march=native.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <assert.h>
#include <stdatomic.h>

#ifndef ALIGN
#define ALIGN 4
#endif

#ifndef GROUP_SIZE
#define GROUP_SIZE 50
#endif

// Number of threads should be a power of 2

#endif

#define DIGS_PER_BM (8u * sizeof(bitmask) / 3u)
#define LSH (3 * (DIGS_PER_BM - 1) - 2)
#define RSH (3 * (DIGS_PER_BM - 1) + 2)

#define SIZE(n) ((size_t)(n) * ((size_t)n + 1) * DIGS_PER_BM * ALIGN * ALIGN / 2)

{
size_t na = n * ALIGN;

if (thd_num == 0)
{
*dest =   (*src & ORIG_MASK) +
+ (*src << 1 & XFER_MASK)
+ (*src >> 5 & XFER_MASK) + (src[1] << LSH & XFER_MASK)
+ 2 * (src[na] >> 2 & XFER_MASK)
+ 3 * (*src >> 5 & 1);           // Origin

s = src + 1;
d = dest + 1;
for (size_t ct = na - 1; ct-- != 0; s++, d++)
+ (*s << 1 & XFER_MASK) + (s[-1] >> RSH)
+ (*s >> 5 & XFER_MASK) + (s[1] << LSH & XFER_MASK)
+ 2 * (s[na] >> 2 & XFER_MASK);
} else
{
s = src + na;
d = dest + na;
}

for (size_t k = n; k != 0; k--, s -= ALIGN, d -= ALIGN)
{
size_t kau = k * ALIGN;
size_t ka = (k - 1) * ALIGN;

if (((n - k) & (NUM_THREADS - 1)) != thd_num)
{
s += kau * (ALIGN * DIGS_PER_BM - (k == n)) + ALIGN;
d += kau * (ALIGN * DIGS_PER_BM - (k == n)) + ALIGN;
continue;
}

for (size_t ct = kau * (ALIGN * DIGS_PER_BM - 1 - (k == n));
ct-- != 0; s1++, d1++)
+ (s1[-kau] >> 2 & XFER_MASK)
+ (*s1 << 1 & XFER_MASK) + (s1[-1] >> RSH)
+ (*s1 >> 5 & XFER_MASK) + (s1[1] << LSH & XFER_MASK)
+ (s1[kau] >> 2 & XFER_MASK);

for (size_t ct = ALIGN - 1; ct-- != 0; s1++, d1++)
*d1 = 0;

+ (s1[-kau] >> 2 & XFER_MASK);

if (ka != 0)
{
*d1++ += (s1++)[1] << LSH & XFER_MASK;

for (size_t ct = ka; ct-- != 1; s1++, d1++)
+ (s1[-kau] >> 2 & XFER_MASK)
+ (*s1 << 1 & XFER_MASK) + (s1[-1] >> RSH)
+ (*s1 >> 5 & XFER_MASK) + (s1[1] << LSH & XFER_MASK)
+ (s1[ka] >> 2 & XFER_MASK);

+ (s1[-kau] >> 2 & XFER_MASK)
+ (*s1 << 1 & XFER_MASK) + (s1[-1] >> RSH)
+ (*s1 >> 5 & XFER_MASK)
+ (s1[ka] >> 2 & XFER_MASK);
}

// Fix up main and left of main diagonals

for (size_t yd = 0; yd != ALIGN; yd++, s++, d++)
{
if (yd != 0)
d[-1] = 0;

for (size_t yf = 0; yf != DIGS_PER_BM; yf++)
{
if (k == n && yd == 0 && yf == 0)

s += kau;
d += kau;
}
if (k != 1 || yd != ALIGN - 1)
d[-kau] += s[1 - kau] << LSH & XFER_MASK;
}
}
}

int get_pixel(const bitmask *src, size_t n, size_t x, size_t y)
{
if (y > x)
{
size_t temp = x;
x = y;
y = temp;
}
size_t y_block = y / (ALIGN * DIGS_PER_BM);
size_t y_remdr = y % (ALIGN * DIGS_PER_BM);
x -= y - y_remdr;
size_t block_start = n * y_block - y_block * (y_block - 1) / 2;
return 3 & src[block_start * ALIGN * ALIGN * DIGS_PER_BM
+ y_remdr * (n - y_block) * ALIGN
+ x / DIGS_PER_BM] >> (3 * (x % DIGS_PER_BM));
}

void print(const bitmask *src, size_t n)
{
long long np = n * ALIGN * DIGS_PER_BM - 1;
for (long long y = -np; y <= np; y++)
{
for (long long x = -np; x <= np; x++)
putchar('0' + get_pixel(src, n, llabs(x), llabs(y)));
}
}

typedef struct
{
size_t n;
unsigned long long weight;
} Grid;

{
size_t ct = SIZE(p_g->n);
for (bitmask *src = p_g->grid; ct-- != 0; src++)
*src <<= 1;
p_g->weight <<= 1;
{
(*p_g->grid)++;
p_g->weight++;
}
}

void expand(Grid *p_g)
{
size_t n = p_g->n;
size_t new_n = ((size_t)(ceil(sqrt(p_g->weight) / 2.7) + 4)
+ DIGS_PER_BM * ALIGN - 1)
/ (DIGS_PER_BM * ALIGN);
if (n < new_n)
{
size_t byte_size = SIZE(new_n) * sizeof(bitmask);
assert(new_grid != NULL);
memset(new_grid, 0, byte_size);
for (size_t k = n, l = new_n; k != 0; k--, l--)
{
for (size_t ct2 = ALIGN * DIGS_PER_BM; ct2-- != 0; )
{
memcpy(dest, src, k * ALIGN * sizeof(bitmask));
src += k * ALIGN;
dest += l * ALIGN;
}
}
free(p_g->grid);
p_g->n = new_n;
p_g->grid = new_grid;
}
}

void initial(Grid *p_g)
{
size_t byte_size = SIZE(1) * sizeof(bitmask);
p_g->n = 1;
p_g->weight = 0;
assert(p_g->grid != NULL);
memset(p_g->grid, 0, byte_size);
}

void copy_empty(const Grid *p_g, Grid *p_h)
{
size_t byte_size = SIZE(p_g->n) * sizeof(bitmask);
p_h->n = p_g->n;
assert(p_h->grid != NULL);
}

void free_grid(Grid *p_g)
{
free(p_g->grid);
}

bool is_done(const Grid *p_g)
{
size_t ct = SIZE(p_g->n);
for (bitmask *src = p_g->grid; ct-- != 0; src++)
if ((*src & HIGH_MASK) != 0)
return false;

return true;
}

Grid grid_1;
Grid grid_2;
bool done_flag = false;

#ifndef USE_COND_VAR

atomic_ullong barr_ct;

void barrier(unsigned long long *p_barr_ct2)
{
if (atomic_fetch_add(&barr_ct, 1) >= *p_barr_ct2 - 1)
return;
thrd_yield();
}

#else

unsigned long long barr_ct = 0;
cnd_t cond;
mtx_t mutex;

void barrier(unsigned long long *p_barr_ct2)
{
int rv = mtx_lock(&mutex);
assert(rv == thrd_success);

if (++barr_ct >= *p_barr_ct2)
{
rv = mtx_unlock(&mutex);
assert(rv == thrd_success);
assert(rv == thrd_success);
return;
}

for (;;)
{
rv = cnd_wait(&cond, &mutex);
assert(rv == thrd_success);
if (barr_ct >= *p_barr_ct2)
{
rv = mtx_unlock(&mutex);
assert(rv == thrd_success);
return;
}
}
/*NOTREACHED*/
}

#endif

void process(unsigned thd_num, unsigned long long *p_barr_ct2)
{
for (int i = 0; i < GROUP_SIZE; i++)
{
one_round(grid_1.grid, grid_2.grid, grid_1.n, thd_num);
barrier(p_barr_ct2);
one_round(grid_2.grid, grid_1.grid, grid_1.n, thd_num);
barrier(p_barr_ct2);
}
}

{
unsigned long long barr_ct2 = 0;
unsigned thd_num = *(unsigned *)p_info;

for (;;)
{
barrier(&barr_ct2);
if (done_flag)
break;
process(thd_num, &barr_ct2);
}

return 0;
}

int main(int argc, char **argv)
{
assert(argc >= 2);
size_t starting_pile = atol(argv[1]);
int rv;

initial(&grid_1);
while (mask << 1 <= starting_pile)

#ifndef USE_COND_VAR

atomic_init(&barr_ct, 0);

#else

rv = cnd_init(&cond);
assert(rv == thrd_success);
rv = mtx_init(&mutex, mtx_plain);
assert(rv == thrd_success);

#endif

for (unsigned i = 0; i < NUM_THREADS - 1; i++)
{
thd_num[i] = i + 1;
assert(rv == thrd_success);
}

unsigned long long barr_ct2 = 0;

{
double_(&grid_1, (mask & starting_pile) != 0);
expand(&grid_1);
copy_empty(&grid_1, &grid_2);
for (;;)
{
barrier(&barr_ct2);
process(0, &barr_ct2);
if (is_done(&grid_1))
break;
}
free_grid(&grid_2);
}

done_flag = true;
barrier(&barr_ct2);
for (unsigned i = 0; i < NUM_THREADS - 1; i++)
{
rv = thrd_join(thd_ids[i], (int *)NULL);
assert(rv == thrd_success);
}

#ifdef USE_COND_VAR

cnd_destroy(&cond);
mtx_destroy(&mutex);

#endif

print(grid_1.grid, grid_1.n);

free_grid(&grid_1);

return 0;
}


Try it online!

# C (gcc), 6348 bytes, score 10500000 on TIO (revised)

This is a revision of my earlier answers to take advantage of the eightfold symmetry of the square. Since I used the doubling trick, the sandpile heights never exceed seven, so they are stored as three-bit octal digits, 21 per unsigned long long.

For TIO I compiled with the flags -O3 -march=native.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <assert.h>

#ifndef ALIGN
#define ALIGN 4
#endif

#ifndef GROUP_SIZE
#define GROUP_SIZE 50
#endif

#define DIGS_PER_BM (8u * sizeof(bitmask) / 3u)
#define LSH (3 * (DIGS_PER_BM - 1) - 2)
#define RSH (3 * (DIGS_PER_BM - 1) + 2)

#define SIZE(n) ((size_t)(n) * ((size_t)n + 1) * DIGS_PER_BM * ALIGN * ALIGN / 2)

{
size_t na = n * ALIGN;

*dest =   (*src & ORIG_MASK) +
+ (*src << 1 & XFER_MASK)
+ (*src >> 5 & XFER_MASK) + (src[1] << LSH & XFER_MASK)
+ 2 * (src[na] >> 2 & XFER_MASK)
+ 3 * (*src >> 5 & 1);           // Origin

const bitmask *s = src + 1;
bitmask *d = dest + 1;
for (size_t ct = na - 1; ct-- != 0; s++, d++)
+ (*s << 1 & XFER_MASK) + (s[-1] >> RSH)
+ (*s >> 5 & XFER_MASK) + (s[1] << LSH & XFER_MASK)
+ 2 * (s[na] >> 2 & XFER_MASK);

for (size_t k = n; k != 0; k--, s -= ALIGN, d -= ALIGN)
{
size_t kau = k * ALIGN;
size_t ka = (k - 1) * ALIGN;

for (size_t ct = kau * (ALIGN * DIGS_PER_BM - 1 - (k == n));
ct-- != 0; s1++, d1++)
+ (s1[-kau] >> 2 & XFER_MASK)
+ (*s1 << 1 & XFER_MASK) + (s1[-1] >> RSH)
+ (*s1 >> 5 & XFER_MASK) + (s1[1] << LSH & XFER_MASK)
+ (s1[kau] >> 2 & XFER_MASK);

for (size_t ct = ALIGN - 1; ct-- != 0; s1++, d1++)
*d1 = 0;

+ (s1[-kau] >> 2 & XFER_MASK);

if (ka != 0)
{
*d1++ += (s1++)[1] << LSH & XFER_MASK;

for (size_t ct = ka; ct-- != 1; s1++, d1++)
+ (s1[-kau] >> 2 & XFER_MASK)
+ (*s1 << 1 & XFER_MASK) + (s1[-1] >> RSH)
+ (*s1 >> 5 & XFER_MASK) + (s1[1] << LSH & XFER_MASK)
+ (s1[ka] >> 2 & XFER_MASK);

+ (s1[-kau] >> 2 & XFER_MASK)
+ (*s1 << 1 & XFER_MASK) + (s1[-1] >> RSH)
+ (*s1 >> 5 & XFER_MASK)
+ (s1[ka] >> 2 & XFER_MASK);
}

// Fix up main and left of main diagonals

for (size_t yd = 0; yd != ALIGN; yd++, s++, d++)
{
if (yd != 0)
d[-1] = 0;

for (size_t yf = 0; yf != DIGS_PER_BM; yf++)
{
if (k == n && yd == 0 && yf == 0)

s += kau;
d += kau;
}
if (k != 1 || yd != ALIGN - 1)
d[-kau] += s[1 - kau] << LSH & XFER_MASK;
}
}
}

int get_pixel(const bitmask *src, size_t n, size_t x, size_t y)
{
if (y > x)
{
size_t temp = x;
x = y;
y = temp;
}
size_t y_block = y / (ALIGN * DIGS_PER_BM);
size_t y_remdr = y % (ALIGN * DIGS_PER_BM);
x -= y - y_remdr;
size_t block_start = n * y_block - y_block * (y_block - 1) / 2;
return 3 & src[block_start * ALIGN * ALIGN * DIGS_PER_BM
+ y_remdr * (n - y_block) * ALIGN
+ x / DIGS_PER_BM] >> (3 * (x % DIGS_PER_BM));
}

void print(const bitmask *src, size_t n)
{
long long np = n * ALIGN * DIGS_PER_BM - 1;
for (long long y = -np; y <= np; y++)
{
for (long long x = -np; x <= np; x++)
putchar('0' + get_pixel(src, n, llabs(x), llabs(y)));
}
}

typedef struct
{
size_t n;
unsigned long long weight;
} Grid;

{
size_t ct = SIZE(p_g->n);
for (bitmask *src = p_g->grid; ct-- != 0; src++)
*src <<= 1;
p_g->weight <<= 1;
{
(*p_g->grid)++;
p_g->weight++;
}
}

void expand(Grid *p_g)
{
size_t n = p_g->n;
size_t new_n = ((size_t)(ceil(sqrt(p_g->weight) / 2.7) + 4)
+ DIGS_PER_BM * ALIGN - 1)
/ (DIGS_PER_BM * ALIGN);
if (n < new_n)
{
size_t byte_size = SIZE(new_n) * sizeof(bitmask);
assert(new_grid != NULL);
memset(new_grid, 0, byte_size);
for (size_t k = n, l = new_n; k != 0; k--, l--)
{
for (size_t ct2 = ALIGN * DIGS_PER_BM; ct2-- != 0; )
{
memcpy(dest, src, k * ALIGN * sizeof(bitmask));
src += k * ALIGN;
dest += l * ALIGN;
}
}
free(p_g->grid);
p_g->n = new_n;
p_g->grid = new_grid;
}
}

void initial(Grid *p_g)
{
size_t byte_size = SIZE(1) * sizeof(bitmask);
p_g->n = 1;
p_g->weight = 0;
assert(p_g->grid != NULL);
memset(p_g->grid, 0, byte_size);
}

void copy_empty(const Grid *p_g, Grid *p_h)
{
size_t byte_size = SIZE(p_g->n) * sizeof(bitmask);
p_h->n = p_g->n;
assert(p_h->grid != NULL);
}

void free_grid(Grid *p_g)
{
free(p_g->grid);
}

bool is_done(const Grid *p_g)
{
size_t ct = SIZE(p_g->n);
for (bitmask *src = p_g->grid; ct-- != 0; src++)
if ((*src & HIGH_MASK) != 0)
return false;

return true;
}

int main(int argc, char **argv)
{
assert(argc >= 2);
size_t starting_pile = atol(argv[1]);
Grid grid_1, grid_2;

initial(&grid_1);
while (mask << 1 <= starting_pile)

{
double_(&grid_1, (mask & starting_pile) != 0);
expand(&grid_1);
copy_empty(&grid_1, &grid_2);
for (;;)
{
for (int i = 0; i < GROUP_SIZE; i++)
{
one_round(grid_1.grid, grid_2.grid, grid_1.n);
one_round(grid_2.grid, grid_1.grid, grid_1.n);
}
if (is_done(&grid_1))
break;
}
free_grid(&grid_2);
}

print(grid_1.grid, grid_1.n);

free_grid(&grid_1);

return 0;
}


Try it online!

• Congs first answer better than hardcode (if it work, not confirmed)
– l4m2
Commented Feb 13, 2022 at 3:29

# Python 3 + numpy, scores ~450000 on TIO

import math
import sys
import numpy as np

def step(g):
a = g>>2
g &= 3
g[1:] += a[:-1]
g[:-1] += a[1:]
g[:, 1:] += a[:, :-1]
g[:, :-1] += a[:, 1:]
g[0] += a[1]
g[:, 0] += a[:, 1]

def nd(g, n):
for y in range(0, n):
for x in range(0, n):
if g[y, x] > 3:
return True
return False

def grd(i, size):
a = i>>2
if a:
grid = grd(a, size)<<2
while nd(grid, size):
for _ in range(a>>(a.bit_length()>>1)):
step(grid)
else:
grid = np.zeros((size, size), dtype=np.ubyte)
grid[0,0] = i&3
return grid

def main(i):
size = math.ceil(math.sqrt(i)/2.7)+3
grid = grd(i, size)
for y in [*range(size-1, 0, -1), *range(size)]:
print(*grid[y][:0:-1], *grid[y], sep='', end='')

if __name__ == '__main__':
i = int(sys.argv[1])
main(i)


Try it online!

Takes the initial height of the sand pile as a command line argument.

This stores and calculates only the bottom right quarter of the grid which can be used to construct the full output in the end. Some inspiration is taken from @alephalpha's comment and a better bound of the required space is taken from their answer.

Many thanks to @att, who found a way to keep the values on the grid at all times below 255 (or below 7 to be exact), which allows using bytes as a integer datatype, which speeds up both memory access and computations.
In rough terms this is achieved by recursively spreading higher powers of four first, then adding the lower ones.

The timing on TIO is a bit strange, 450000 finishes in 15 seconds, but 500000 times out. Local timings:

$time python3 run.py 1300000 > /dev/null python3 run.py 1300000 > /dev/null 59.94s user 0.32s system 99% cpu 1:00.33 total  • I tried this, which reaches 300000 on TIO. But using numba will make it slower. I don't know why. Commented Feb 7, 2022 at 5:38 • @alephalpha Interesting, this does seem to perform much better on TIO than my first revision, but on my laptop there is only a small difference in performance. And when I adapt a similar approach for only calculating a quarter of the grid (using 2d arrays again for simplicity) I get the same performance on TIO but double (~900000) locally. I guess this is related to memory speed / available cache? Numba's JIT can't do much here because the step function only contains very few calls to numpy's C functions, which Python is already good enough at. – ovs Commented Feb 7, 2022 at 7:45 • faster? – att Commented Feb 8, 2022 at 18:30 • @att thanks a lot, edited. Again this seems to be most significant on TIO, but still a nice improvement! – ovs Commented Feb 9, 2022 at 9:23 • I believe the maximum possible stored value will be 15: 3<<2=12 can be generated from the recursive step, and a number with residue 3 mod 4, surrounded by four 12s (or higher) will become 15 on the next step. – att Commented Feb 9, 2022 at 19:12 # Charcoal + tr, score circa 2560 on TIO python3 /opt/charcoal/charcoal.py -c ＵＢ0ＦＮ«≔⟦Ｅχ⁰⟧υＦυ«Ｊ⊟κ⊟κ≔﹪⊕ΣＫＫ⁴θＩθ¿¬θＦ⁴«Ｍ✳⁺³⊗λ⊞υ⟦ⅉⅈ -i$1 | tr -d \\n


Try it online! Verbose Charcoal source: Try it online! Takes i as input and outputs a square grid. Explanation: Simply increments the centre cell i times, then when it overflows adds each orthogonally adjacent cell to the list of cells to increment.

# Python 3, score circa 5120 on TIO

import collections
import sys
import numpy
c = collections.Counter()
c[0j] = int(sys.argv[1])
[(m, n)] = c.most_common(1)
while n > 3:
c[m] = n % 4
for d in [1, 1j, -1, -1j]:
c[m + d] += n // 4
[(m, n)] = c.most_common(1)
l = max(int(m.real) for m in c)
print(end = ''.join(str(c[i + j * 1j]) for j in range(-l, l + 1) for i in range(-l, l + 1)))


Try it online! Explanation: Keeps track of the sand using a Counter indexed by Gaussian integers. The centre cell starts off with the program's argument and then the cell with the largest value gets overflowed each time until no cell has a value greater than 3.

# C (gcc), 3282 bytes

Not optimized at all, but it should still be quite fast, because it's C.

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

#define swap(a, b, TYPE) \
do {                 \
TYPE tmp = a;    \
a = b;           \
b = tmp;         \
} while (0)

long long int* data_a = NULL;
long long int* data_b = NULL;
int size_x = 3;
int size_y = 3;

void copy_a_into_b(int size) {
for (int i = 0; i < size; i++) {
data_b[i] = data_a[i];
}
}

void larger() {
int new_size_x = size_x + 2;
int new_size_y = size_y + 2;
copy_a_into_b(size_x * size_y);
data_a = (long long int*)realloc(
data_a, new_size_x * new_size_y * sizeof(long long int));
for (int i = 0; i < new_size_x * new_size_y; i++) {
data_a[i] = 0;
}
int offset_x = 1;
int offset_y = 1;
for (int i = 0; i < size_x; i++) {
for (int j = 0; j < size_y; j++) {
data_a[(i + offset_x) * new_size_y + j + offset_y] =
data_b[i * size_y + j];
}
}
data_b = (long long int*)realloc(
data_b, new_size_x * new_size_y * sizeof(long long int));
copy_a_into_b(new_size_x * new_size_y);
size_x = new_size_x;
size_y = new_size_y;
}

int is_stable() {
for (int i = 0; i < size_x; i++) {
for (int j = 0; j < size_y; j++) {
if (data_a[i * size_y + j] > 3) {
return 0;
}
}
}
return 1;
}

void zero(long long int* target, int size) {
for (int i = 0; i < size; i++) {
target[i] = 0;
}
}

for (int i = 0; i < size; i++) {
data_b[i] += data_a[i];
}
}

void update() {
zero(data_b, size_x * size_y);
int should_grow = 0;
for (int i = 0; i < size_x; i++) {
for (int j = 0; j < size_y; j++) {
if (data_a[i * size_y + j] > 3) {
if (i <= 1 || j <= 1 || i >= size_x - 2 || j >= size_y - 2) {
should_grow = 1;
}
int quarter = data_a[i * size_y + j] >> 2;
data_a[i * size_y + j] = data_a[i * size_y + j] & 3;
// data_b[i * size_y + j] = 0;
data_b[(i - 1) * size_y + j] += quarter;
data_b[(i + 1) * size_y + j] += quarter;
data_b[i * size_y + j - 1] += quarter;
data_b[i * size_y + j + 1] += quarter;
}
}
}
swap(data_a, data_b, long long int*);
copy_a_into_b(size_x * size_y);
if (should_grow) {
larger();
}
}

int main(int argc, char* argv[]) {
int height = atoi(argv[1]);
int size = (int)sqrt((double)height) + 1;
if (size % 2 == 0) {
size++;
}
size_x = size;
size_y = size;
data_a = (long long int*)malloc(sizeof(long long int) * size_x * size_y);
data_b = (long long int*)malloc(sizeof(long long int) * size_x * size_y);
zero(data_a, size_x * size_y);
zero(data_b, size_x * size_y);
copy_a_into_b(size_x * size_y);
int x = size / 2;
int y = size / 2;
data_a[x * size_y + y] = height;
int q = 0;
while (!is_stable()) {
q++;
update();
}
for (int i = 0; i < size_x * size_y; i++) {
printf("%d", (int)data_a[i]);
}
fflush(stdout);
return 0;
}


# Wolfram Language (Mathematica), 284 bytes, score ~150000 on TIO

Fixed a bug thanks to @att.

With[{size = 2 Ceiling[Sqrt@#/2.7] + 3},
WriteString[$Output, ##] & @@ FixedPoint[ BitAnd[#, 3] + With[{q = BitShiftRight[#, 2]}, Sum[RotateLeft[q, d], {d, {1, -1, size, -size}}]] &, CenterArray[#, size^2]]] &@ToExpression[$ScriptCommandLine[[2]]]


Try it online!

• Since OP disallows a trailing newline, you'd need WriteString or BinaryWrite instead of Print
– att
Commented Feb 7, 2022 at 22:45
• Mathemagica. (filler) Commented Feb 8, 2022 at 15:10

# MATLAB, 310 bytes, score ~2^16-2^17

function piles = sand(h)

piles = h;

while any(piles > 3,'all')

piles = conv2([0 1/4 0; 1/4 0 1/4; 0 1/4 0],4*floor(piles/4)) + conv2([0 0 0; 0 1 0; 0 0 0],mod(piles,4).*(piles>=4)) + conv2([0 0 0; 0 1 0; 0 0 0],piles.*(piles<4));
piles(all(~piles,2),:) = [];
piles(:,all(~piles,1)) = [];
end


Try it online!

I did this on a relatively average PC, so I don't know if better hardware could get this done any faster. something on the order of 2^16 was the best I could get in less than a minute, 2^17 took almost 2 minutes.

# C (clang), score 1376256

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define u8 unsigned char

int main(int argc, char* argv[]) {
if (argc != 2) {
fprintf(stderr, "Please provide i (and only i)\n");
return -1;
}
int height = atoi(argv[1]);
if (height < 1) {
return -1;
}

int maxm = 0;
while (height >> ++maxm){};

// Whatever... These just have to be big enough
int size = height + 1000;
int width = height + 1000;

u8* tmpbuf = malloc(width);
u8* mainbuf = calloc(size, 1);

// Null ptrs are zero, deal with it
u8** bufptrs = calloc(__builtin_popcount(height) - 1, sizeof(u8*));
int bufptrsidx = 0;

mainbuf[0] = 1;

int index = height;
int rcount = 0;
int grandmax = 1;
while (1) {
int stuff_happening = 1;
int lmax = grandmax;
while (stuff_happening) {
stuff_happening = 0;
int p = 0;
int l = 0;
int next_lmax = 0;
while (l < lmax) {
int lw = l + 1;
u8 carry = 0;
for (int j = 0; j < lw; j++) {
u8 k = mainbuf[p+j] / 4;
mainbuf[p+j]%= 4;
mainbuf[p+j]+= carry;
if (j + 1 == lw) {
mainbuf[p+j]+= carry;
} else {
mainbuf[p+1-lw+j]+= k;
mainbuf[p+j]+= tmpbuf[j];
}
if (j + 2 == lw) {
mainbuf[p+1-lw+j]+= k;
if (j == 0) {
mainbuf[p+1-lw+j]+= 2*k;
}
}
tmpbuf[j] = k;
if (j > 0) {
mainbuf[p+j-1]+= k;
}
if (j == 1) {
mainbuf[p+j-1]+= k;
}
carry = k;
}
p+= lw;
l+= 1;
if (l == lmax) {
lmax+= 1;
}
next_lmax = l + 1;
stuff_happening|= 1;
}
}
if (lmax > grandmax) {grandmax = lmax;}
lmax = next_lmax;
}
int tmp = index & 1;
if (index/= 2) { // Double until largets power of two reached
// If this is needed in the future, make a copy
if (tmp) {
bufptrs[bufptrsidx++] = memcpy(malloc(width), mainbuf, width);
}
for (int i = 0;i<size;++i){mainbuf[i]*= 2;}
} else if (rcount < bufptrsidx) { // When that is done, add old arrays starting from 1
for (int i = 0;i<size;++i){mainbuf[i]+= bufptrs[rcount][i];}
rcount+= 1;
} else { // Exit
break;
}
}

for (int i = 0;i < bufptrsidx;i++){free(bufptrs[i]);}
free(bufptrs);
free(tmpbuf);

for (int y = 0; y < grandmax*2-1; y++) {
int yp = grandmax - y - 1;
if (yp < 0){yp*=-1;}
for (int x = 0; x < grandmax*2-1; x++) {
int xp = grandmax - x - 1;
if (xp < 0){xp*=-1;}
int r,c;
if (xp > yp) {
r = xp;
c = yp;
} else {
r = yp;
c = xp;
}
int ro = r*(r+1)/2;
u8 ch = mainbuf[ro+c];
putchar(ch + 48);
}
// Newlines for pretty print
//putchar(10);
}
free(mainbuf);
}


Try it online!

This answer demonstrates some techniques that may be useful.

First, we only do the calculations on one eighth of the grid, since there is eight-fold symmetry.

The second technique we use is the doubling trick (first noticed by @att). If we wish to calculate the resulting sand pile for input $$\2i\$$, we can first calculate the resulting sand pile for input $$\i\$$, double every value, and iterate until we are done. Similarily, we can calculate the resulting sandpile for input $$\a+b\$$ by calculating $$\a\$$ and $$\b\$$ individually and then just adding them together and iterating.

Combined with the observation that the maximum height of a sandpile doesn't increase (assuming you do "global" iterations), we can use just a byte to store the height of a sandpile, because the height of a sandpile never exceeds seven.

Currently the two main easiest ways to make this answer faster are multi-threading and use of SSE instructions.

For scoring, I used this manually unrolled version, since it was a bit faster when compiled with clang:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define u8 unsigned char

int main(int argc, char* argv[]) {
if (argc != 2) {
fprintf(stderr, "Please provide i (and only i)\n");
return -1;
}
int height = atoi(argv[1]);
if (height < 1) {
return -1;
}

int maxm = 0;
while (height >> ++maxm){};

// Whatever... These just have to be big enough
int size = height + 1000;
int width = height + 1000;

u8* tmpbuf = malloc(width);
u8* mainbuf = calloc(size, 1);

// Null ptrs are zero, deal with it
u8** bufptrs = calloc(__builtin_popcount(height) - 1, sizeof(u8*));
int bufptrsidx = 0;

mainbuf[0] = 1;

int index = height;
int rcount = 0;
int grandmax = 1;
while (1) {
int stuff_happening = 1;
int lmax = grandmax;
while (stuff_happening) {
stuff_happening = 0;
int p = 0;
int l = 0;
int next_lmax = 0;
while (l < lmax && l < 4) {
int lw = l + 1;
u8 carry = 0;
for (int j = 0; j < lw; j++) {
u8 k = mainbuf[p+j] / 4;
mainbuf[p+j]%= 4;
mainbuf[p+j]+= carry;
if (j + 1 == lw) {
mainbuf[p+j]+= carry;
} else {
mainbuf[p+1-lw+j]+= k;
mainbuf[p+j]+= tmpbuf[j];
}
if (j + 2 == lw) {
mainbuf[p+1-lw+j]+= k;
if (j == 0) {
mainbuf[p+1-lw+j]+= 2*k;
}
}
tmpbuf[j] = k;
if (j > 0) {
mainbuf[p+j-1]+= k;
}
if (j == 1) {
mainbuf[p+j-1]+= k;
}
carry = k;
}
p+= lw;
l+= 1;
if (l == lmax) {
lmax+= 1;
}
next_lmax = l + 1;
stuff_happening|= 1;
}
}
while (l < lmax && l >= 4) {
// At least 5
int lw = l + 1;
u8 carry = 0;
for (int j = 0; j < 2; j++) {
u8 k = mainbuf[p+j] / 4;
mainbuf[p+j]%= 4;
mainbuf[p+j]+= carry;
if (j + 1 == lw) {
mainbuf[p+j]+= carry;
} else {
mainbuf[p+1-lw+j]+= k;
mainbuf[p+j]+= tmpbuf[j];
}
if (j + 2 == lw) {
mainbuf[p+1-lw+j]+= k;
if (j == 0) {
mainbuf[p+1-lw+j]+= 2*k;
}
}
tmpbuf[j] = k;
if (j > 0) {
mainbuf[p+j-1]+= k;
}
if (j == 1) {
mainbuf[p+j-1]+= k;
}
carry = k;
}
for (int j = 2; j < lw - 2; j++) {
u8 k = mainbuf[p+j] / 4;
mainbuf[p+j]%= 4;
mainbuf[p+j]+= carry;
if (j + 1 == lw) {
mainbuf[p+j]+= carry;
} else {
mainbuf[p+1-lw+j]+= k;
mainbuf[p+j]+= tmpbuf[j];
}
if (j + 2 == lw) {
mainbuf[p+1-lw+j]+= k;
if (j == 0) {
mainbuf[p+1-lw+j]+= 2*k;
}
}
tmpbuf[j] = k;
if (j > 0) {
mainbuf[p+j-1]+= k;
}
if (j == 1) {
mainbuf[p+j-1]+= k;
}
carry = k;
}
for (int j = lw - 2; j < lw; j++) {
u8 k = mainbuf[p+j] / 4;
mainbuf[p+j]%= 4;
mainbuf[p+j]+= carry;
if (j + 1 == lw) {
mainbuf[p+j]+= carry;
} else {
mainbuf[p+1-lw+j]+= k;
mainbuf[p+j]+= tmpbuf[j];
}
if (j + 2 == lw) {
mainbuf[p+1-lw+j]+= k;
if (j == 0) {
mainbuf[p+1-lw+j]+= 2*k;
}
}
tmpbuf[j] = k;
if (j > 0) {
mainbuf[p+j-1]+= k;
}
if (j == 1) {
mainbuf[p+j-1]+= k;
}
carry = k;
}
p+= lw;
l+= 1;
if (l == lmax) {
lmax+= 1;
}
next_lmax = l + 1;
stuff_happening|= 1;
}
}
if (lmax > grandmax) {grandmax = lmax;}
lmax = next_lmax;
}
int tmp = index & 1;
if (index/= 2) { // Double until largets power of two reached
// If this is needed in the future, make a copy
if (tmp) {
bufptrs[bufptrsidx++] = memcpy(malloc(width), mainbuf, width);
}
for (int i = 0;i<size;++i){mainbuf[i]*= 2;}
} else if (rcount < bufptrsidx) { // When that is done, add old arrays starting from 1
for (int i = 0;i<size;++i){mainbuf[i]+= bufptrs[rcount][i];}
rcount+= 1;
} else { // Exit
break;
}
}

for (int i = 0;i < bufptrsidx;i++){free(bufptrs[i]);}
free(bufptrs);
free(tmpbuf);

for (int y = 0; y < grandmax*2-1; y++) {
int yp = grandmax - y - 1;
if (yp < 0){yp*=-1;}
for (int x = 0; x < grandmax*2-1; x++) {
int xp = grandmax - x - 1;
if (xp < 0){xp*=-1;}
int r,c;
if (xp > yp) {
r = xp;
c = yp;
} else {
r = yp;
c = xp;
}
int ro = r*(r+1)/2;
u8 ch = mainbuf[ro+c];
putchar(ch + 48);
}
// Newlines for pretty print
//putchar(10);
}
free(mainbuf);
}
$$$$

• Quadrupling (rather than doubling) is convenient because any residue mod 4 will always be located at the center due to symmetry, making the add step much faster.
– att
Commented Feb 9, 2022 at 19:22
• @att Sorry, I don't quite understand. Could you try to explain again? At the start the piles are obviously divisible by four, but after the first iteration, the edges will contain piles whose heights are not divisible by four, right? Also, how does divisibility make the add step any faster? Commented Feb 9, 2022 at 19:32
• I think we might have slightly different conceptions of the "trick", but in any case: The original idea I had in mind was that we can first recursively settle $\lfloor\frac i4\rfloor$. 4 is natural due to the fourfold symmetry; any residue mod 4 must lie at the origin, as the sum of all values outside the origin must be divisible by 4. Adding the residue back in would then be just one operation. Thinking again, mod 2 should work similarly; you'd need to |= the residue instead of =.
– att
Commented Feb 9, 2022 at 22:33

# Java (JDK), 1418 bytes

import java.util.Scanner;
import java.lang.Math;
public class Main {
public static void main(String[] args) {
Scanner scanner = new Scanner(System.in);
int middleSize = scanner.nextInt();
int sideLength = 2 * (int) (Math.sqrt(middleSize) / 2) + 3;

int[][] grid = new int[sideLength][sideLength];
int[][] tempGrid = new int[sideLength][sideLength];

for (int i = 0; i < sideLength; i++) {
for (int j = 0; j < sideLength; j++) {
if (i == (int) (sideLength / 2) && j == (int) (sideLength / 2)) {
grid[i][j] = middleSize;
tempGrid[i][j] = middleSize;
} else {
grid[i][j] = 0;
tempGrid[i][j] = 0;
}
}
}

while (true) {
for (int i = 0; i < sideLength; i++) {
for (int j = 0; j < sideLength; j++) {
tempGrid[i][j] = grid[i][j] % 4;
if (i > 0) tempGrid[i][j] += (int) grid[i-1][j] / 4;
if (i < sideLength - 1) tempGrid[i][j] += (int) grid[i+1][j] / 4;
if (j > 0) tempGrid[i][j] += (int) grid[i][j-1] / 4;
if (j < sideLength - 1) tempGrid[i][j] += (int) grid[i][j+1] / 4;
}
}
boolean stable = true;
for (int i = 0; i < sideLength; i++) {
for (int j = 0; j < sideLength; j++) {
grid[i][j] = tempGrid[i][j];
if (grid[i][j] > 3) stable = false;
}
}
if (stable) break;
}

for (int i = 0; i < sideLength; i++) {
for (int j = 0; j < sideLength; j++) {
System.out.print(grid[i][j]);
}
}
}
}


Try it online!

I somehow managed to put together a functional answer to this challenge, in spite of my very limited Java knowledge. (So limited in fact that I had to look up the names of libraries which I knew about but whose names had been forgotten due to me barely ever using Java.) It seems to work correctly, and the score is likely around 150000 because it took 55 seconds with an input of 150000.

Builds a grid with a size roughly the same as the input ("size" meaning area) and keeps tossing the sand around until no more sand can be tossed around. Not the most optimized, but 150000 is not a bad score.

BTW, I chose Java because it's a compiled language and usually compiled languages are faster.

# Rust, 1903 bytes, score ~1000000 on TIO, ~3500000 on my computer

A simple port of @ovs's Python answer, without @att's improvement. It is fast on TIO, but much slower than @ovs's answer on my computer.

Now only compute one eighth of the grid, and borrow the trick in @Aiden's answer: first bring the maximum value to less than 256, and then switch to u8.

macro_rules! sandpile {
($grid:ident,$size:ident, $boundary:expr) => {{ let mut quater =$grid.clone();
while $grid.iter().any(|&x| x >$boundary) {
quater.clear();
quater.extend($grid.iter().map(|&x| x >> 2));$grid.iter_mut().for_each(|x| *x &= 3);
$grid[0] += 4 * quater[1]; for i in 1..$size - 1 {
for j in i + 1..$size - 1 {$grid[i * $size + j] += quater[(i + 1) *$size + j]
+ quater[i * $size + (j + 1)] + quater[(i - 1) *$size + j]
+ quater[i * $size + (j - 1)]; }$grid[i] += 2 * quater[i + $size] + quater[i + 1] + quater[i - 1];$grid[i * $size + i] += 2 * quater[(i - 1) *$size + i] + 2 * quater[i * $size + (i + 1)]; } }$grid
}};
}

fn sandpile32(n: u32, size: usize) -> Vec<u32> {
let mut grid = vec![0; size * size];
grid[0] = n;

sandpile!(grid, size, 255)
}

fn sandpile8(grid32: &[u32], size: usize) -> Vec<u8> {
let mut grid = Vec::with_capacity(grid32.len());
grid.extend(grid32.iter().map(|&x| x as u8));

sandpile!(grid, size, 3)
}

fn print(grid: &[u8], size: usize) {
let isize = size as isize;
for i in 1 - isize..isize {
for j in 1 - isize..isize {
let (i, j) = (i.abs() as usize, j.abs() as usize);
let index = std::cmp::min(i, j) * size + std::cmp::max(i, j);
print!("{}", grid[index]);
}
}
}

fn main() {
let n = std::env::args()
.nth(1)
.and_then(|x| x.parse::<u32>().ok())
.expect("Give me a number");
let size = ((n as f32).sqrt() / 2.7) as usize + 3;
let grid32 = sandpile32(n, size);
let grid = sandpile8(&grid32, size);
print(&grid, size);
}


Try it online!

# Rust, 2044 bytes, score ~750000 on TIO, ~6500000 on my computer

This one is a port of @Polichinelle's C answer, and borrow the trick in @Aiden's answer: first bring the maximum value to less than 256, and then switch to u8.

It is slow on TIO. But it become much faster when compiled with -Ctarget-cpu=native on my computer.

macro_rules! sandpile {
($grid:ident,$size:ident, $boundary:expr) => {{ let mut buf =$grid.clone();
while $grid.iter().any(|&x| x >$boundary) {
std::mem::swap(&mut $grid, &mut buf);$grid[0] = (buf[0] & 3) + (buf[1] & !3);

for i in 1..$size - 1 {$grid[i] = (buf[i] & 3)
+ 2 * (buf[i + $size] >> 2) + (buf[i + 1] >> 2) + (buf[i - 1] >> 2); } for i in$size..$size * ($size - 1) {
$grid[i] = (buf[i] & 3) + (buf[i + 1] >> 2) + (buf[i +$size] >> 2)
+ (buf[i - 1] >> 2)
+ (buf[i - $size] >> 2); } for i in 1..$size - 1 {
$grid[i *$size] = (buf[i * $size] & 3) + 2 * (buf[i *$size + 1] >> 2)
+ (buf[(i + 1) * $size] >> 2) + (buf[(i - 1) *$size] >> 2);
}

for i in 1..$size {$grid[i * $size - 1] = 0; } }$grid
}};
}

fn sandpile32(n: u32, size: usize) -> Vec<u32> {
let mut grid = vec![0; size * size];
grid[0] = n;

sandpile!(grid, size, 255)
}

fn sandpile8(grid32: &[u32], size: usize) -> Vec<u8> {
let mut grid = Vec::with_capacity(grid32.len());
grid.extend(grid32.iter().map(|&x| x as u8));

sandpile!(grid, size, 3)
}

fn print(grid: &[u8], size: usize) {
let isize = size as isize;
for i in 1 - isize..isize {
for j in 1 - isize..isize {
let (i, j) = (i.abs() as usize, j.abs() as usize);
print!("{}", grid[i * size + j]);
}
}
}

fn main() {
let n = std::env::args()
.nth(1)
.and_then(|x| x.parse::<u32>().ok())
.expect("Give me a number");
let size = ((n as f32).sqrt() / 2.7) as usize + 3;
let grid32 = sandpile32(n, size);
let grid = sandpile8(&grid32, size);
print(&grid, size);
}


Try it online!

# F# 6, 3521 bytes, score ~ 33000 on my computer

Requires .NET 6. TIO doesn't support it, so I can't provide a sample run of it.

The algorithm is total junk, but I wanted to try and implement it in as pure a functional way as possible - no mutable keyword, all value/record types, no changing items in collections. Turns out that it's not a good idea, but it is a lot of fun.

namespace Sand

open FSharp.Collections.ParallelSeq

type Coordinates =
{
Row: int
Column: int
}

{
Coordinates: Coordinates
Value: int
}

module SandPile =
let maxHeight = 4

|> PSeq.groupBy (fun x -> x.Coordinates)
|> PSeq.map (fun (_, changes) ->
changes |> Seq.reduce (fun total next -> { total with Value = total.Value + next.Value })
)

let distribution = originalAdjustment.Value / maxHeight
seq {
yield { Adjustment.Coordinates = { originalAdjustment.Coordinates with Row = originalAdjustment.Coordinates.Row - 1 }; Value = distribution  }
yield { Adjustment.Coordinates = { originalAdjustment.Coordinates with Column = originalAdjustment.Coordinates.Column - 1 }; Value = distribution  }
yield { Adjustment.Coordinates = { originalAdjustment.Coordinates with Row = originalAdjustment.Coordinates.Row + 1 }; Value = distribution  }
yield { Adjustment.Coordinates = { originalAdjustment.Coordinates with Column = originalAdjustment.Coordinates.Column + 1 }; Value = distribution  }
}
else
)

|> PSeq.tryFind (fun x -> x.Value >= maxHeight)
|> Option.isNone

let rows = adjustments |> PSeq.map (fun x -> x.Coordinates.Row) |> Seq.sort |> Seq.toArray
let columns = adjustments |> PSeq.map (fun x -> x.Coordinates.Column) |> Seq.sort |> Seq.toArray

let lowRow = rows |> Array.head
let highRow = rows |> Array.last
let lowColumn = columns |> Array.head
let highColumn = columns |> Array.last

let gridDict = adjustments |> PSeq.groupBy (fun x -> x.Coordinates.Row, x.Coordinates.Column) |> Map.ofSeq

seq {
for row = lowRow to highRow do
for column = lowColumn to highColumn do
yield
Map.tryFind
(row, column)
gridDict
|> Option.map (fun x -> x |> Seq.head |> fun x -> x.Value |> string)
|> Option.defaultValue "0"
}
|> String.concat " "
|> printf "%s"

let loop = (calculateDeltas >> applyAdjustments >> PSeq.toArray)

module Program =

[<EntryPoint>]
let main argv =
argv[0]
|> int
|> fun initialSeed -> { Adjustment.Coordinates = { Row = 0; Column = 0 }; Value = initialSeed }
|> Array.singleton
|> Seq.unfold (fun surface ->
let isStable = SandPile.isStable surface
if isStable then None
else
let newSurface = SandPile.loop surface
Some (newSurface, newSurface)
)
|> Seq.last
|> SandPile.printOut

0


# Rust, 10667 bytes

Significantly faster than my last answer, and comfortably handles 4000000 on my machine. The idea is still the same as my previous answer, by first reducing everything to under 255 and then operating on bytes. I dropped nalgebra and platform intrinsics this time and added a few nightly features, mainly for portable SIMD types. Additionally, only a quarter of the grid is calculated now. Building this one with just rustc instead of cargo is fine but ensure -Ctarget-cpu=native is passed to rustc.

#![feature(portable_simd)]
#![feature(allocator_api)]
#![feature(array_chunks)]
use std::mem;
use std::ops::{Index, IndexMut};
use std::ptr::NonNull;
use std::{env::args, fmt::Debug, simd::*};
fn main() {
let arg = args()
.nth(1)
.expect("No input provided")
.trim()
.parse::<u32>()
.expect("Not parsable as an unsigned integer");
let mut sandpile = Mat::sandpile(arg);
sandpile.topple::<255>();
let mut sandpile = sandpile.map(|n| n as u8);
sandpile.topple::<3>();
sandpile.trim();
sandpile.unquarter();
sandpile.print();
}
fn max_size(n: u32) -> usize {
((n as f64).sqrt() / 2.7).ceil() as usize + 3
}

#[derive(Clone)]
struct Mat<T> {
data: Vec<T, AlignedAlloc<SIMD_ALIGN>>,
row_size: usize,
}

impl<T> Index<(usize, usize)> for Mat<T> {
type Output = T;
fn index(&self, (x, y): (usize, usize)) -> &T {
&self.data[x + y * self.row_size]
}
}
impl<T> IndexMut<(usize, usize)> for Mat<T> {
fn index_mut(&mut self, (x, y): (usize, usize)) -> &mut T {
&mut self.data[x + y * self.row_size]
}
}

impl<T> Mat<T> {
fn row(&self, y: usize) -> &[T] {
&self.data[y * self.row_size..(y + 1) * self.row_size]
}
fn row_mut(&mut self, y: usize) -> &mut [T] {
&mut self.data[y * self.row_size..(y + 1) * self.row_size]
}
fn row_count(&self) -> usize {
self.data.len() / self.row_size
}
fn row_windows(&self) -> RowWindows<'_, T> {
RowWindows {
data: &self.data,
row_size: self.row_size,
}
}
fn map<U: Default + Copy, F: Fn(T) -> U>(self, f: F) -> Mat<U> {
let mut buf = Vec::with_capacity_in(self.data.len(), AlignedAlloc);
buf.extend(self.data.into_iter().map(f));
let mut mat = Mat {
data: buf,
row_size: self.row_size,
};
mat
}
}

macro_rules! make_topple {
($int:ty,$simd:ty, $lanes:literal) => { fn topple<const MAX:$int>(&mut self) {
if !self.data.iter().any(|&x| x > MAX) {
return;
}
let mut buf = self.clone();
let four = <$simd>::splat(4); loop { let mut max = <$simd>::splat(0);
for ((up, center, down), buf) in self
.row_windows()
.zip(buf.data.chunks_exact_mut(self.row_size).skip(1))
{
for (((((cur, left), right), up), down), buf) in center[1..]
.array_chunks::<$lanes>() .zip(center.array_chunks::<$lanes>())
.zip(center[2..].array_chunks::<$lanes>()) .zip(up[1..].array_chunks::<$lanes>())
.zip(down[1..].array_chunks::<$lanes>()) .zip(buf[1..].array_chunks_mut::<$lanes>())
{
let cur = <$simd>::from_array(*cur); let left = <$simd>::from_array(*left);
let right = <$simd>::from_array(*right); let up = <$simd>::from_array(*up);
let down = <$simd>::from_array(*down); let sum = left / four + right / four + up / four + down / four + cur % four; max = max.simd_max(sum); *buf = sum.to_array(); } } self.row(2)[1..] .array_chunks::<$lanes>()
.zip(buf.row_mut(1)[1..].array_chunks_mut::<$lanes>()) .for_each(|(cur, buf)| { let cur = <$simd>::from_array(*cur);
let tmp = <\$simd>::from_array(*buf);
let sum = tmp + cur / four;
max = max.simd_max(sum);
*buf = sum.to_array();
});
let col_max = (1..self.row_count() - 1)
.map(|y| {
buf[(1, y)] += self[(2, y)] / 4;
buf[(1, y)]
})
.max()
.unwrap();
mem::swap(self, &mut buf);
let max = max.reduce_max().max(col_max);
if max <= MAX {
break;
}
}
}
#[allow(dead_code)] // u32 version is unused
fn trim(&mut self) {
let mut num_rows = 0;
for row in self.data.chunks_exact(self.row_size).skip(1) {
if row.iter().all(|&x| x == 0) {
break;
}
num_rows += 1;
}
let mut buf = Vec::with_capacity_in(num_rows * num_rows, AlignedAlloc);
buf.extend(std::iter::repeat(0).take(num_rows * num_rows));
let mut new = Mat {
data: buf,
row_size: num_rows,
};
for x in 0..num_rows {
for y in 0..num_rows {
new[(x, y)] = self[(x + 1, y + 1)];
}
}
*self = new;
}
};
}

impl<T> Mat<T>
where
T: Debug + Default,
{
fn print(&mut self) {
self.data.truncate(self.row_size * self.row_size);
for n in self.data.iter() {
print!("{n:?}");
}
}
}
impl<T: Default + Copy> Mat<T> {
let round_mask = SIMD_ALIGN / mem::align_of::<T>() - 1;
let padding = SIMD_ALIGN / mem::align_of::<T>();
let size = self.row_size;
let aligned_size = if size & round_mask == 0 {
size
} else {
};
let mut data = Vec::with_capacity_in(
AlignedAlloc::<SIMD_ALIGN>,
);
data.extend(
std::iter::repeat_with::<T, _>(Default::default)
);
let mut new = Mat {
data,
};
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[(x, y)] = self[(x, y)];
}
}
*self = new;
}
fn unquarter(&mut self) {
let mut new = Mat {
data: Vec::with_capacity_in(
self.row_size * self.row_count() * 4,
AlignedAlloc::<SIMD_ALIGN>,
),
row_size: self.row_size * 2 - 1,
};

new.data.extend(
std::iter::repeat_with::<T, _>(Default::default)
.take((self.row_size * 2 - 1).pow(2) - 1 + self.row_size * 2 - 1),
);
let size = new.row_size;
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[((size - 1) / 2 + x, (size - 1) / 2 + y)] = self[(x, y)];
}
}
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[(x, y)] = self[(self.row_size - 1 - x, self.row_count() - 1 - y)];
}
}
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[((size - 1) / 2 + x, y)] = self[(x, self.row_count() - 1 - y)];
}
}
for x in 0..self.row_size {
for y in 0..self.row_count() {
new[(x, (size - 1) / 2 + y)] = self[(self.row_size - 1 - x, y)];
}
}
*self = new;
}
}
impl Mat<u32> {
fn sandpile(n: u32) -> Self {
let size = max_size(n);
let mut data = Vec::with_capacity_in(size * (size), AlignedAlloc::<SIMD_ALIGN>);
data.extend(std::iter::repeat(0).take(size * size));
let mut mat = Mat {
data,
row_size: size,
};
mat[(1, 1)] = n;
mat
}
make_topple!(u32, u32x8, 8);
}
impl Mat<u8> {
make_topple!(u8, u8x32, 32);
}

pub struct RowWindows<'a, T> {
data: &'a [T],
row_size: usize,
}
impl<'a, T> Iterator for RowWindows<'a, T> {
type Item = (&'a [T], &'a [T], &'a [T]);
fn next(&mut self) -> Option<Self::Item> {
if self.data.len() < self.row_size * 3 {
return None;
}
let (first, rest) = self.data.split_at(self.row_size);
self.data = rest;
let (second, rest) = rest.split_at(self.row_size);
let (third, _rest) = rest.split_at(self.row_size);

Some((first, second, third))
}
}

impl<T: Debug> Debug for Mat<T> {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
for i in 0..self.data.len() / self.row_size {
for j in 0..self.row_size {
write!(f, "{:?} ", self.data[i * self.row_size + j])?;
}
writeln!(f)?;
}
Ok(())
}
}

const SIMD_ALIGN: usize = mem::align_of::<u32x8>();

#[derive(Clone, Copy, Debug)]
struct AlignedAlloc<const ALIGN: usize>;
use std::alloc::{AllocError, Allocator, Global, Layout};
unsafe impl<const ALIGN: usize> Allocator for AlignedAlloc<ALIGN> {
fn allocate(&self, layout: Layout) -> Result<NonNull<[u8]>, AllocError> {
Global.allocate(layout.align_to(ALIGN).unwrap())
}

unsafe fn deallocate(&self, ptr: NonNull<u8>, layout: Layout) {
Global.deallocate(ptr, layout.align_to(ALIGN).unwrap())
}

fn allocate_zeroed(&self, layout: Layout) -> Result<NonNull<[u8]>, AllocError> {
Global.allocate_zeroed(layout.align_to(ALIGN).unwrap())
}

unsafe fn grow(
&self,
ptr: NonNull<u8>,
old_layout: Layout,
new_layout: Layout,
) -> Result<NonNull<[u8]>, AllocError> {
Global.grow(
ptr,
old_layout.align_to(ALIGN).unwrap(),
new_layout.align_to(ALIGN).unwrap(),
)
}

unsafe fn grow_zeroed(
&self,
ptr: NonNull<u8>,
old_layout: Layout,
new_layout: Layout,
) -> Result<NonNull<[u8]>, AllocError> {
Global.grow_zeroed(
ptr,
old_layout.align_to(ALIGN).unwrap(),
new_layout.align_to(ALIGN).unwrap(),
)
}

unsafe fn shrink(
&self,
ptr: NonNull<u8>,
old_layout: Layout,
new_layout: Layout,
) -> Result<NonNull<[u8]>, AllocError> {
Global.shrink(
ptr,
old_layout.align_to(ALIGN).unwrap(),
new_layout.align_to(ALIGN).unwrap(),
)
}
}


Playground

I tried parallelizing this with rayon but the code got lost somewhere in synchronization primitives and took roughly ten times longer. I might return to this later and see if I can do better.

This code is fully cross-platform, and in fact, I cross-compiled it for 32-bit arm at one point. Because this uses several nightly features it likely won't work on the latest nightly forever. If it doesn't, here's my current rustc --version:

rustc 1.66.0-nightly (57f097ea2 2022-10-01)

• Do you have a proof that the max won’t increase to ≥ 256 after decreasing to < 256? Commented Feb 10, 2022 at 23:53
• @AndersKaseorg expand out [[0,255,0],[255,255,255],[0,255,0]]` once to see why I'm pretty sure this is the case. Other answers also seem to be relying on this fact. Commented Feb 11, 2022 at 5:22
• @AndersKaseorg Assuming you do "global" iterations (meaning in the same way as suggested in the post), then the maximum height of a sandpile (divided by four) is always decreasing. You can think if it's possible for a cell to reach a value of 256 or higher, given that it's neighbours are of value <256 Commented Feb 11, 2022 at 7:34