# 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!!
– null
Feb 6 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) Feb 6 at 14:56
• This is not gonna win any prices for speed, but here is some animated sand dropping
– ovs
Feb 6 at 15:14
• Why can't I output a square grid?
– Neil
Feb 6 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
Feb 7 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
Feb 13 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. Feb 7 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 Feb 7 at 7:45 • faster? – att Feb 8 at 18:30 • @att thanks a lot, edited. Again this seems to be most significant on TIO, but still a nice improvement! – ovs Feb 9 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 Feb 9 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
Feb 7 at 22:45
• Mathemagica. (filler)
– null
Feb 8 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
Feb 9 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? Feb 9 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
Feb 9 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 + Nalgebra, 7392 bytes

Should be plenty fast, although there is still a ton of room for optimization (among other issues, it still calculates the full grid). The strategy I took was based on two observations:

1. The maximum value of the sandpile never increases

2. The vast majority of the iterations occur when the maximum value is less than 256

Based on those two things, I came up with a two step algorithm:

1. use a fairly naive algorithm to bring the maximum value to less than 256
2. store one cell per byte in a 256 bit vector and go nuts with simd

The code will panic if you try to run it on a machine without avx2 support, and on debug mode will pretty print the grid rather than output all the elements in a row. Build cargo rustc --release -- -Ctarget-cpu=native and run ./target/release/sandpiles <arg>. Fair warning, the first build does take a few minutes while cargo builds the dependencies. Also, if you want to use rust-analyzer with this code, you're going to want to turn off the mismatched-arg-count and unresolved-macro-call errors.

src/main.rs

use nalgebra::*;
use std::env::args;
fn main() {
assert!(is_x86_feature_detected!("avx2"));
let arg = args()
.nth(1)
.expect("No input provided")
.parse::<u32>()
.expect("Not parseable as an unsigned integer");
let s = sandpile(arg);
#[cfg(debug_assertions)]
{
println!("{}", trim(unpack_simd(s)));
}
#[cfg(not(debug_assertions))]
{
let s = trim(unpack_simd(s));
write_mat(&s);
}
}
fn max_size(n: u32) -> usize {
let val = ((n as f64).sqrt() / 2.7).ceil() as usize + 3;
if val % 2 == 0 {
val + 1
} else {
val
}
}
fn sandpile(n: u32) -> DMatrix<SimdInt> {
let mut matrix = dmatrix![n];
reserve(&mut matrix, max_size(n));
matrix = step(matrix, 256);
let mat = step_simd(pack_simd(matrix));
mat
}
#[inline]
fn step(mut matrix: DMatrix<u32>, max: u32) -> DMatrix<u32> {
while matrix.max() >= max {
let size = matrix.ncols();
let mut collapsed = (&matrix).map(|a| a % 4);
let piles = matrix / 4;
let mut slice = collapsed.slice_range_mut(1.., ..);
slice += piles.slice_range(..size - 1, ..);
slice = collapsed.slice_range_mut(..size - 1, ..);
slice += piles.slice_range(1.., ..);
slice = collapsed.slice_range_mut(.., 1..);
slice += piles.slice_range(.., ..size - 1);
slice = collapsed.slice_range_mut(.., ..size - 1);
slice += piles.slice_range(.., 1..);

matrix = collapsed;
}
matrix
}

fn trim(mut matrix: DMatrix<u32>) -> DMatrix<u32> {
let trim_count = matrix.row_iter().take_while(|r| r.max() == 0).count();
matrix = matrix.remove_rows(0, trim_count);
let trim_count = matrix.row_iter().skip_while(|r| r.max() > 0).count();

let rows = matrix.nrows();
matrix = matrix.remove_rows(rows - trim_count, trim_count);
let trim_count = matrix.column_iter().take_while(|c| c.max() == 0).count();
matrix = matrix.remove_columns(0, trim_count);
let trim_count = matrix.column_iter().skip_while(|c| c.max() > 0).count();
let cols = matrix.ncols();
matrix = matrix.remove_columns(cols - trim_count, trim_count);
matrix
}

#[inline]
fn reserve(mat: &mut DMatrix<u32>, additional: usize) {
if !matches!(
Some(0)
) {
let increase = additional * 2;
let size = mat.ncols();
let mut expanded = DMatrix::zeros(size + increase, size + increase);
let mut slice =
slice += &*mat;
*mat = expanded;
}
}
fn write_mat(mat: &DMatrix<u32>) {
print!("{}", mat.iter().map(u32::to_string).collect::<String>())
}
use std::arch::x86_64::*;
#[repr(C)]
#[derive(Copy, Clone)]
union SimdInt {
a: [u8; 32],
s: __m256i,
}
use std::fmt;
impl fmt::Debug for SimdInt {
fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
f.debug_list().entries(unsafe { self.a }).finish()
}
}
impl PartialEq for SimdInt {
#[inline(always)]
fn eq(&self, other: &Self) -> bool {
unsafe { self.a == other.a }
}
}
impl SimdInt {
#[inline(always)]
unsafe fn store(self, ptr: *mut Self) {
_mm256_store_si256(ptr.cast(), self.s)
}
#[inline(always)]
fn any_ge_four(self) -> bool {
unsafe { self.a.into_iter().any(|n| n >= 4) }
}
#[inline(always)]
fn mod_four(self) -> Self {
let mask = SimdInt { a: [0x03; 32] };
}
#[inline(always)]
fn div_four(self) -> Self {
let mask = SimdInt { a: [!0x03; 32] };
unsafe {
SimdInt {
}
}
}
}
impl BitAnd for SimdInt {
type Output = Self;
#[inline(always)]
fn bitand(self, rhs: Self) -> Self::Output {
SimdInt {
s: unsafe { _mm256_and_si256(self.s, rhs.s) },
}
}
}
impl From<SimdInt> for [u8; 32] {
fn from(a: SimdInt) -> Self {
unsafe { a.a }
}
}

type Output = Self;
#[inline(always)]
fn add(self, other: Self) -> Self {
SimdInt {
s: unsafe { _mm256_add_epi8(self.s, other.s) },
}
}
}
#[inline(always)]
fn add_assign(&mut self, rhs: Self) {
unsafe { (*self + rhs).store(self) }
}
}
fn pack_simd(mut mat: DMatrix<u32>) -> DMatrix<SimdInt> {
debug_assert!(mat.max() < 256);
mat = mat.transpose();
let cols = mat.ncols();
let padding = 32 - cols % 32;

let vec = mat
.row_iter()
.map(|r| {
let mut packed = DVector::from_element(r.len() / 32, SimdInt { a: [0; 32] });
for i in 0..(cols + padding) / 32 {
let mut arr = [0u8; 32];
let j = 32 * i;
for k in 0..32 {
arr[k] = r[j + k] as u8;
}
packed[i] = SimdInt { a: arr };
}
packed
})
.collect::<Vec<_>>();
DMatrix::from_columns(&vec)
}
fn unpack_simd(mat: DMatrix<SimdInt>) -> DMatrix<u32> {
let final_row_len = mat.nrows() * 32;
let final_col_len = mat.ncols();
let mut out = DMatrix::zeros(final_col_len, final_row_len);
for (c, mut r) in mat.column_iter().zip(out.row_iter_mut()) {
for i in 0..mat.nrows() {
let j = i * 32;
let packed: [u8; 32] = c[i].into();
for k in 0..32 {
r[j + k] = packed[k] as u32;
}
}
}
out
}
fn step_simd(mut mat: DMatrix<SimdInt>) -> DMatrix<SimdInt> {
let size = mat.ncols();
while mat.iter().copied().any(SimdInt::any_ge_four) {
for _ in 0..256 {
let mut piles = mat.map(SimdInt::mod_four);
let spill = mat.map(SimdInt::div_four);
let mut slice = piles.slice_range_mut(.., 1..);
slice += spill.slice_range(.., ..size - 1);
slice = piles.slice_range_mut(.., ..size - 1);
slice += spill.slice_range(.., 1..);
let rows = piles.nrows();
let slice = piles.as_mut_slice();
for (left, right) in slice.chunks_exact_mut(rows).zip(spill.column_iter()) {
let (_, center, _) = unsafe { left.align_to_mut::<u8>() };
center.rotate_right(1);
let mut slice = unsafe {
VectorSliceMut::<SimdInt, Dynamic>::from_slice_unchecked(left, 0, rows)
};
slice += right;
let left = slice.as_mut_slice();
let (_, center, _) = unsafe { left.align_to_mut::<u8>() };
center.rotate_left(2);
let mut slice = unsafe {
VectorSliceMut::<SimdInt, Dynamic>::from_slice_unchecked(left, 0, rows)
};
slice += right;
let left = slice.as_mut_slice();
let (_, center, _) = unsafe { left.align_to_mut::<u8>() };
center.rotate_right(1);
}
mat = piles;
}
}

mat
}



Cargo.toml:

[package]
name = "sandpiles"
version = "0.1.0"
edition = "2021"

[dependencies]
nalgebra = "0.30"
[profile.release]
lto = "fat"


Playground

• Do you have a proof that the max won’t increase to ≥ 256 after decreasing to < 256? Feb 10 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. Feb 11 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 Feb 11 at 7:34

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