[Question]
  [
In this [fastest-code](/questions/tagged/fastest-code "show questions tagged 'fastest-code'") 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\$:
[![Sand fractal 3*10^7](https://upload.wikimedia.org/wikipedia/commons/e/ee/Sandpile_on_infinite_grid%2C_3e7_grains.png)](https://upload.wikimedia.org/wikipedia/commons/e/ee/Sandpile_on_infinite_grid%2C_3e7_grains.png)
[Credit colt\_browning CC BY 4.0](https://commons.wikimedia.org/wiki/File:Sandpile_on_infinite_grid,_3e7_grains.png)
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.
# Further reading
These are called [Abelian sandpiles](https://en.wikipedia.org/wiki/Abelian_sandpile_model). Numberphile has a [video](https://www.youtube.com/watch?v=1MtEUErz7Gg) about them.
# Leaderboard
| Score | Language | Author |
| --- | --- | --- |
| 10000000 | C (clang -march=native -Ofast) | [Polichinelle](https://codegolf.stackexchange.com/a/242612/84290) |
| 8388608 | Bash + xz + node (compressed) | [Community](https://codegolf.stackexchange.com/a/242526/84290) |
| 4800000 | Rust (rustc -Copt-level=3 -Ctarget-cpu=native) | [alephalpha](https://codegolf.stackexchange.com/a/242598/84290) |
| 1572864 | Python + numpy | [ovs](https://codegolf.stackexchange.com/a/242442/84290) |
| 1376256 | C (clang -march=native -Ofast) | [AnttiP](https://codegolf.stackexchange.com/a/242610/84290) |
| 1350000 | Rust (RUSTFLAGS="-Ctarget-cpu=native" cargo ...) | [Aiden4](https://codegolf.stackexchange.com/a/242695/84290) |
| 410000 | C (clang -march=native -O3) | [astroide](https://codegolf.stackexchange.com/a/242446/84290) |
 
  
[Answer]
# [C (clang)](http://clang.llvm.org/), 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>
#include <threads.h>
typedef unsigned long long bitmask;
#ifndef ALIGN
#define ALIGN 4
#endif
#ifndef GROUP_SIZE
#define GROUP_SIZE 50
#endif
// Number of threads should be a power of 2
#ifndef NUM_THREADS
#define NUM_THREADS 8
#endif
#define ORIG_MASK 0333333333333333333333ull
#define XFER_MASK 0111111111111111111111ull
#define HIGH_MASK 0444444444444444444444ull
#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)
void one_round(const bitmask *src, bitmask *dest, size_t n, unsigned thd_num)
{
  size_t na = n * ALIGN;
  const bitmask *s;
  bitmask *d;
  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++)
      *d =    (*s & ORIG_MASK)
            + (*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;
    }
    const bitmask *s1 = s;
    bitmask *d1 = d;
    for (size_t ct = kau * (ALIGN * DIGS_PER_BM - 1 - (k == n));
         ct-- != 0; s1++, d1++)
      *d1 =   (*s1 & ORIG_MASK)
            + (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;
    *d1 =  (*s1 & ORIG_MASK)
         + (s1[-kau] >> 2 & XFER_MASK);
    if (ka != 0)
    {
      *d1++ += (s1++)[1] << LSH & XFER_MASK;
      for (size_t ct = ka; ct-- != 1; s1++, d1++)
        *d1 =   (*s1 & ORIG_MASK)
              + (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);
      *d1 =   (*s1 & ORIG_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)
          continue; // Skip top row
        bitmask diag_mask = 3ull << (yf * 3);
        *d += s[-kau] >> 2 & XFER_MASK & diag_mask;
        *d += *s >> 5 & XFER_MASK & diag_mask;
        *d &= ~(diag_mask >> 3);
        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;
  bitmask *grid;
} Grid;
void double_(Grid *p_g, bool add_one)
{
  size_t ct = SIZE(p_g->n);
  for (bitmask *src = p_g->grid; ct-- != 0; src++)
    *src <<= 1;
  p_g->weight <<= 1;
  if (add_one)
  {
    (*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);
    bitmask *new_grid =
           (bitmask *)aligned_alloc(ALIGN * sizeof(bitmask), byte_size);
    assert(new_grid != NULL);
    memset(new_grid, 0, byte_size);
    bitmask *src = p_g->grid;
    bitmask *dest = new_grid;
    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;
  p_g->grid = (bitmask *)aligned_alloc(ALIGN * sizeof(bitmask), byte_size);
  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;
  p_h->grid = (bitmask *)aligned_alloc(ALIGN * sizeof(bitmask), byte_size);
  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;
}
// Common data between threads
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)
{
  *p_barr_ct2 += NUM_THREADS;
  if (atomic_fetch_add(&barr_ct, 1) >= *p_barr_ct2 - 1)
    return;
  while (atomic_load(&barr_ct) < *p_barr_ct2)
    thrd_yield();
}
#else
unsigned long long barr_ct = 0;
cnd_t cond;
mtx_t mutex;
void barrier(unsigned long long *p_barr_ct2)
{
  *p_barr_ct2 += NUM_THREADS;
  int rv = mtx_lock(&mutex);
  assert(rv == thrd_success);
  if (++barr_ct >= *p_barr_ct2)
  {
    rv = mtx_unlock(&mutex);
    assert(rv == thrd_success);
    rv = cnd_broadcast(&cond);
    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);
  }
}
int helper_thread(void *p_info)
{
  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(sizeof(bitmask) == 8);
  assert(argc >= 2);
  size_t starting_pile = atol(argv[1]);
  size_t mask = 1;
  int rv;
  initial(&grid_1);
  while (mask << 1 <= starting_pile)
    mask <<= 1;
#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
  thrd_t thd_ids[NUM_THREADS - 1];
  unsigned thd_num[NUM_THREADS - 1];
  for (unsigned i = 0; i < NUM_THREADS - 1; i++)
  {
    thd_num[i] = i + 1;
    rv = thrd_create(&thd_ids[i], helper_thread, &thd_num[i]);
    assert(rv == thrd_success);
  }
  unsigned long long barr_ct2 = 0;
  for (; mask != 0; mask >>= 1)
  {
    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!](https://tio.run/##zVltb9vIEf6uX7FBUB2pl1h0cmgARQbS2GcHTeyDfWmLGgFBkStpa2qpklQsNef@9Koz@84XyU56B1SALXJ3dnb22WdmZ1bxME4jPt/tnjMep@uEkjdFmbDsxeKkU2lK2bTeljM@r7Yto3LRGDnNsrTaGBUFzcuGYFRmSxZXm8tFTqOkwMZOuV3RhM7ImhdszmlC0ozP5b8pK5dRcTfuwNAZR6G3H96fX3aewyPjVL6RV53nlCdsZqXOr68@/RzevP/7mRG1TeTHkRlwdEQu18spzUk2I8ooUiyydZqQKSURWWX3svPYar/89DH85eL67O3pjVHvtJHX1h7Ve3X9/jz8@Pbmz2T0su2zTlMj@7efzq6VbND2cWUv3p9fKNlXbR9X9vT9@U34M@j@00fivV6THinYv2g28xTIPjkiL9e@kf9wc0G8lyDmuSOHJPDh37GVu94v10c5I4jYe9wnnocTh6WPLz37ykE@wBZXT0/tsf4@Eiq/ZCwhGadhnq154sUZL0pNFtIr8nhg3xJalAMi5yB8YGlWLpKQr5d@52uHmP6ITAjXswHvCKkrH0ObVS5E2Ix4ShuZTMjIhzZUSuTsoJIQD80iXUsFgEeI6E9fibx5QwKQMzTwW6VOTsiPFSnsgY7b4DNqwM07oOMYcUdpHn1GVceHhMXeupMG/tgRABe6ytmc8Y4YV8BqURb2ciwaEmgQKJiWWZYTtekkRnQAdeDLGF6GQ/IMAByTot8fkKTf18b0EoEiGlIBsQlOE0ABze0wEEsFurYNasfzW9BsxxJX/EBoWlDDCYsQjxoQyaYHxNKF6Q5RGsOXROduOARKk@FE8hSQMs@We3potIbBd5bTlT7o8u6kt7qkl5z2PA5dEBm6xHMDHIr7aIpxIRzxVcFSkP5EzAqoaOethgaYEdyE@wiyYxQC8d1DwU1LxtdUvj/IRdR9N0DopYR1YWxM1Kob1DxoTOAaNLascHkcCCIHLpMDHRCCg1QugtshTP@4gwpN7awPDtM@2MP74GnER8F2C/fBKXFsOPsekEZKjYLsAGIH4XIYDYR/pgK0JWwPZ0bmeWiI3752paSVInYxQdtinrrn37Lr373vv8XO273fj/b/PdE7T16PE1DgsPuJbch6RZYR4yTikK3SWYkJomhIWDTPeJQWTf5vE8Fo/H6mHAFfkCzVw07TEvkqhUfW2ESs07pGbY6ZmmOGw5xghU2WkF@NOuETIoCRbleYCOPF48xkMiaq6QiLKNzcsRUpsxXJs/uOEdJBFWEIxdOEYHKLm@aByh55afkhDnXwumLP7sOz0TOujWk5sfeKdyfk3561CMZVjNDHlR2T1FoeOi5Y6OXk11/dfRQHortFYj1iaXhEiLe2iKKVP3SAXoyXZE7LcMU2NG1NaW0Oq5425mkrs1hBGXJCNo0soKTLFWzGRk66gcetfNzCI3bKtMPIb8NpmsW4f1vIudvOP@EYRjqnyyQX0n/YL73BLGULiCh5R4OYLSzKKC9VBq4tGJonOIdtY4D1yjFqyGm5zjmkqV1Mq25dTfXqoWJSLaSJKKAXAlNxO7PJjdqGbMAOR6vgsayHNoCFCwFg8KBqlxUU2eXBTZYbakthvnIrk2YqMtYZox2CWzvkK/B98gbG4oOMAV9tcLLSGy290dIbGzFW6zJeRLn3w@gHWLKlqTAZCJmm0bTwNr5@2voyIRLMtiV@UebruKwUXCjVUvnfUzZflJVaa54zSNMeyLn4ljgm2Xqa0tDDNtJbhXMo/LIsJVGShFAfVmo7cViLIhTkhifcN5C5OwAyolvMVklS8ljjoeq0iURdiEt7bSN6ojFCI@71jGq/35f@54yWTQ@GJHSzggPGrq1aqWpDueNFnN6H2GFL7Jgy2KV/5qXnTCRc58Uf8eR85Tc5feDTby3OK@Gv@TmqXg2oUb6GiZM30vBG0JpuSxris944Kda8uqil9CiGIJOJa5XdZj9KBd/CKAX3NuGqpnRg51cTyAsuz6gHYlx@@vBB9S7psqC2d0BGTQ17mVarSeStgVbVLJtFPQjOpoTqpWE6HFZziWrSemwy8V41SYAuw/dmsgALjFdbT12ooOffOfGohp5TD4lit1GCiqNSFL0TWEet66Fj/89ySj3rOI7bcLN82yb3vQKd9SjGWcmidI9LNegW7KGambzh/yIxq1ryv9JOkc6qdFmnOGc6G6TTK4@z1TaEY77cqlPHiZj6cXEYDRU190GykJDYmCSafmsQFg0Q9AqRJmLL67vb4A@MEIcEK8IE4nMdkN/p0BAXKuoS0Nzb@pX0XmUyM6gfqMhQVQOcmnQsT1JxZf0uWy4zKDaiMiJTWt5TyvXldacj1oGGhMBO@wKpklg0rjicpdEcbNYT6YvtTzdn4bury9PwL2@vOx15bx9C@i6u4qM8D@NSH7z4ymjutRzdgKESPpZQOg3o7M5Nkjkp5VQzCllGCMem11UDBpjnnUwqKsxpI9FBHfcLllKjJs0iq8GH06ViEY4EsJJwy2iaeJIPz8UNXaftJwg5UHp2zBPkRMYhrizLDTwv1yXd/A6gQDGQf4FJcRZMQb2umMl1BuyfyKUU6zimReGb2@h@X9tdBc8esUb7mtf1H57BDEYwpjlgHUdF6XURlacONxtnLzrH45ptqP4@YkrzgHyTgYjBfgTsoXYQhcencZeiSznId3qXV7/AVr67ODvtHUl6yd@DTPafoRKv/lPEgDyFOQIuJAiTxT4Dhttft@DdTfLtjyQyJLyQp4QMCe5L8ALyCX2fq9IRRWfHgHGb0pqe71JqSuAFTVc0D2U48wRcgADjs0yufr@LHuvTtw4rtPcs1j1fqhu3Uk@b160tGQllQqeO11OwUd016C01O9mtLc9G89FYrxWvjMRORvkcsimssEivBy9f5FoV/eq/0wEXX7uBAEcjy4/dklwUwIzPoVBL8QiH2Jii5Jfb4LMrp25pAht1ZBBRmVJX7qPvRFkxQtzCQbFTmUbiovqFzj1HCyEqVOM0TrAfYQRTsdgJA1LKBJjDXmlcWo4SDj0QDasU8H48gmpnVcdEKWjEkuK29hPI5zautQoJlhlJx21rwjXf1SoZXvgx@/uZWKGwLQYCltTragvZ50HVgYCHVsuTYqdg6qNOZlxH7rZMdtT12kSez3INukBXPBoo@nRrxJFpkLRPlb0u84ibvxpV8kF7qOPKNr7vc2frsKO6q2pv19mhtsNWuI7b2wpFpp6uSQJJN92Sadw@o76TJJYM/8ggmrhMEKGl5zsl6qNbD@7a5q3ohFir5dnW8UP0KdOqD07HfeQlV@upIDOVOm6BP67Fyd1uF7waic9/YoSx2A2vZpBv7IanDi6T490wXcLfStJ@N/wrlBfwn27KPILvFU0iXrJ4N1xGebyY8KhkX@huWJTJJA6C/wI "C (clang) – Try It Online")
# [C (gcc)](https://gcc.gnu.org/), 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>
typedef unsigned long long bitmask;
#ifndef ALIGN
#define ALIGN 4
#endif
#ifndef GROUP_SIZE
#define GROUP_SIZE 50
#endif
#define ORIG_MASK 0333333333333333333333ull
#define XFER_MASK 0111111111111111111111ull
#define HIGH_MASK 0444444444444444444444ull
#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)
void one_round(const bitmask *src, bitmask *dest, size_t n)
{
  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++)
    *d =    (*s & ORIG_MASK)
          + (*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;
    const bitmask *s1 = s;
    bitmask *d1 = d;
    for (size_t ct = kau * (ALIGN * DIGS_PER_BM - 1 - (k == n));
         ct-- != 0; s1++, d1++)
      *d1 =   (*s1 & ORIG_MASK)
            + (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;
    *d1 =  (*s1 & ORIG_MASK)
         + (s1[-kau] >> 2 & XFER_MASK);
    if (ka != 0)
    {
      *d1++ += (s1++)[1] << LSH & XFER_MASK;
      for (size_t ct = ka; ct-- != 1; s1++, d1++)
        *d1 =   (*s1 & ORIG_MASK)
              + (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);
      *d1 =   (*s1 & ORIG_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)
          continue; // Skip top row
        bitmask diag_mask = 3ull << (yf * 3);
        *d += s[-kau] >> 2 & XFER_MASK & diag_mask;
        *d += *s >> 5 & XFER_MASK & diag_mask;
        *d &= ~(diag_mask >> 3);
        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;
  bitmask *grid;
} Grid;
void double_(Grid *p_g, bool add_one)
{
  size_t ct = SIZE(p_g->n);
  for (bitmask *src = p_g->grid; ct-- != 0; src++)
    *src <<= 1;
  p_g->weight <<= 1;
  if (add_one)
  {
    (*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);
    bitmask *new_grid =
           (bitmask *)aligned_alloc(ALIGN * sizeof(bitmask), byte_size);
    assert(new_grid != NULL);
    memset(new_grid, 0, byte_size);
    bitmask *src = p_g->grid;
    bitmask *dest = new_grid;
    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;
  p_g->grid = (bitmask *)aligned_alloc(ALIGN * sizeof(bitmask), byte_size);
  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;
  p_h->grid = (bitmask *)aligned_alloc(ALIGN * sizeof(bitmask), byte_size);
  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(sizeof(bitmask) == 8);
  assert(argc >= 2);
  size_t starting_pile = atol(argv[1]);
  size_t mask = 1;
  Grid grid_1, grid_2;
  initial(&grid_1);
  while (mask << 1 <= starting_pile)
    mask <<= 1;
  for (; mask != 0; mask >>= 1)
  {
    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!](https://tio.run/##zVh7b9s2EP/fn4JFMU@yo8ZyWqyA4gAd1qbB0qZIUGxYUAiyRNuCZcqT5MRe6n30eXd8iNQr67YOmIDEFHkP3u@OxzuFzjwMD4enMQuTTUTJaV5EcfpscdarTCXxtD6XxWxenVsFxaLBOU3TpDoZ5DnNCpzrFbs1jeiMbFgezxmNSJKyufg3jYtVkC@9HnDOGBK9urw4f997CsOYUfFGnveeUhbFM011fn318YN/c/HL65JUT5EXI80gV6@uL879d69ufiSjk7ZnkyQl7c9vXl9LWrftMWnfXpy/lbTP2x6T9oeL8xv/A8j@/h2xXm7IgOTxbzSdWRIFmxyTk41d0l/evCXWCZBZJqdDXBv@jTXddTfdEOlKQgTHYjaxLFTsFza@DPQrA3oXZ0w5A@kE9XvMRd6lcURSRv0s3bDIClOWF8qbZJBn4ZF@i2heHBGhgzC799Aj5VtAJoQp2RAGRJDDLCEWyiF97TuwBwjUM5QEp6fEBarSa3YLzdkZeVGhwRVYuHU/IT8i3SlhjBAhLQs@oaBxNyl3gqnQtT29To6PyVUWz2OGdtYhA5ORETzgwaoGD@Y5InJhlmZE@ouEiBNACK724MVxyJMJGXkkHw6PSDQcis1xEQLOCph1mJpAcpBuHZebDVHWZGnH9ctRbceUx4Fp6BLt9OBH2Ld0HIgn4kxE2ICt5RiVPHBFijXYAPNSh1hlDZaspTgqZgw2veOiewSzdg1ORpKh4RdUDCaqw1M7mvAHiidglm17GhjTiS73oqvcSKRC7ki3y5PcMbl764D6x6NVutDtcLvb4feSrd317l/7Xm2xfYddcAocG5HeAdJIipGQPYLYo3BJKfEMnBVwlYLzQSsbDslwgjJgC@22SyGtIaKNcduM@VKf/x2v/2O/fw3Pa993o/2/D/TeF9uDBHthFeT@N/GWbNZkFcSMBAzKIDorSDoTE1EczFMWJHkz/ncRj2j8fSIPAr5gsFQzvQpLjFdBPNKbjbid@mjUdMykjhmyGckKp3RAPpTi@JngCYz0@3yLwM@HMz40UYJkWsRsQz1E4WYZr0mRrkmW3vdKIpVUEQafjyYEizJ0mgUiB@RExwe/0@DU5R3eh3Epx6vxtFxaneT9Cfnd0jsCvsomcpQHG9A8UW1m3zPBwlNOPn82/civHtNF3B5uGl4R/K0toyjh@x6EV8wKMqeFv463NGktxVS1VY625Wgn6jEeMuSMbBsXaEFXa3DGVijdwnAnhjsY4qLXE3tRAv1pkobovx3Uim33Hz8YJXVGV1HGqb/ppt7iBb8DRCS9IYFr8/MiyApZS6odOOUI7mE96WKdPUYJGS02GYOqrY@l160pqV71VrZUS2k8CyhDQBXTmsuyoo1lC/swpPI4FnX8FrAwIQAM9rLmXkNHVjzqZOFQ3WOxtVljN0uRsqrULOhah63h7JNT4MWByAEPOjlp6q2i3irqrc4Y600RLoLM@nb0LZisw5RvGQIySYJpbm1tNdrZoiDika17R@hEN2FRaR2QqqWlvKfxfFFUauh5FkOZtifn/FfgGKWbaUJ9C@fIYO3PoWGBJpYEUeRDX1PpUvhlzZsnoHPOmF1CZnoAaPgy11YpUrKwrMVFwzIRqHNysV89iSex3IRC3BqUou3hUJw/g1tM7csgods1XDDatmrPpTbKjFPE6L2PC7o1DGkMXvo1KyxDET86z77Dm/O53YzpR55ha1NZSX/N57ja0kouW8HEyKnYeCNpTXcF9XGsHCfImi13raRHMgSZTMxdaTfbQcLjzQ8SON5luqoJPdL6pQLxNcQqxUNgvP94eSlXV3SVU716REZNCZ2RVutJRPOsRHmNYoK3UnDYJFG9q0ocp1pLVIvWcVmJD6pFAiyV8d4sFsDAcL2z5IcAPPlLIx/V0DP6Id4QN7o3flXylngCdtSW9j39f5ZRaumDYxwbVpqv54TfK9DpExWzuIiDpONINcLN7Qi1Unnj/PPCrLqTfxt2Mui0SDPqZMyVi42gU5aH6XrnwzVf7OStY2RMNVw8jobMml2QLAQkOifxqa8NwqIBgrIQw4S7vO7dRvwAB78k4tyPID/XAfmPLg1MdepbWPm90a6U97KSmUH/QHmFKifg1qSeqhGxz7BwEGRzOIJ4LZPBAF7uxMYlVPWPklDNvzShRG5yNiFjs47jVVPM5nC7J@j3oEgTpLyDrtCkk6U9D3@OGxruu0fid8z3rs5aX6xx9vsFyrU4O@/j4Lqs6BRIyHWuQMHtiVkBqyzkJ@LeERlKlQJ9tRehpl9TIAAXCUNesOYOiXlSSlFiMLaNROx5LRkW3RKLFiyGW01/zIb3tvZLf3oVmp6JMyy0mS/uM2Yk1BpbjbKDzexiVOQry/XtPc1osPRq2VccKxME3g6LIrZLa6/J69pmVEOm3B8OB3f0YoTPH@EsCeb5wfkJMgP8p9siC@AXKscAGs/w4ORFNAld9@BcnRycVZCFiwkLiviOHpxk9Sc "C (gcc) – Try It Online")
[Answer]
# [Python 3](https://docs.python.org/3/) + [numpy](https://numpy.org/), 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!](https://tio.run/##dVNNb6MwEL3zK@bUmNZhIWxVKQoc9xfsDSHLXRzHUnBc43RD/zw7xlBIq50L43nz8d7YmN6dLjofBtWai3XQcneKJr/ru9nV19b0wDvQJooacYTOCUNkvI8AjUMBsix340HCQwF5cKtsX8NTAbzab7N6ink3BBGdYxSWVArr7HD6hJaadO6ySk1XiXVgqhsiKeiJ6vFioQelwXItBUkXZEZv/0W9qSMO6incaighv8e8WeGuVsNvexXR6vyLnzsR@EjbEEWhUx9itT41rw8H8KWttKrxy8UaPtUcDrtP@O9JncWoEPPueq4VsUURL0vCk1fl2Flo6U4kLsssjr/rCPeLXeMREkj/Gyttkg9hLx0hfvA0nkLjeiMKBK@vvROh3ldUKcXrQakP@Xo1HgqbabnSRE1kfDNM9u8x@SPUmYxe92YdpvzYJS/xUx59WdG81vubrh6Ddg9tM3wjFLYZ8lyF43oRZ6zSjjyOjPu62qf@@WH2FMARwhSbDQWhG/zGUYRXxpjmrWAMigI2jHkljG1CU@U1Y0v8nRJu5Ts@2MBw0jsMw8/nFO0f "Python 3 – 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](https://codegolf.stackexchange.com/a/242477/64121).
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
```
[Answer]
# Charcoal + tr, score circa 2560 on TIO
```
python3 /opt/charcoal/charcoal.py -c ＵＢ0ＦＮ«≔⟦Ｅχ⁰⟧υＦυ«Ｊ⊟κ⊟κ≔﹪⊕ΣＫＫ⁴θＩθ¿¬θＦ⁴«Ｍ✳⁺³⊗λ⊞υ⟦ⅉⅈ -i $1 | tr -d \\n
```
[Try it online!](https://tio.run/##S0oszvj/v6CyJCM/z1hBP7@gRD85I7EoOT8xB87QK6hU0E1WeL9n6/s9iwze71n2fs@6Q6sfdU55NB/IXHq@/VHjhkfzl59vBUqdbz20@v2eVY@65p/bBSE6p7zfCeRPPbf4/R6g1OpHjVvO7Xi/Z@W5HYf2H1oDYi4DCoF0rX00Z/Ojxl2HNj/qmn5u96OueedbgTY8au181NqhoJupoGKoUKNQUqSgm6IQE5P3//9/MzMA "Bash – Try It Online") Verbose Charcoal source: [Try it online!](https://tio.run/##NY2xasMwEIZn@SlEphOo4NCQpVPbLA04CNKhJWRQbCURlnWOpMtS8uyuKlM4cZz4/@9rrzq0qN007U16021/CUi@g0W9EC/VGQOHDz9S2tFwMgGE4D8Ve43RXjwcGj3Cspa8FkfJKedZKVAJsS0N4yeCwhF6Ifm8/0L//QY7cpj5bTCD8cl0sKcBlDF9FuXKKr9baahgfYJ3HRPcZoY9c9hhOXmxrmYra/BuYGODaZNFD8pRhGfJN0gnlw0ukwuAKYpXIMkP35A9XyCO5f9R5XlM03o9Pd3dLw "Charcoal – 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.
[Answer]
# [Python 3](https://docs.python.org/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!](https://tio.run/##fY7BTgMhEIbvPMVcTMFuqaRNDyZ68THIxmxYbCHMsGGpuk@PsDXGg/FAAvN9P/9MS75EOpTicIopg4khWJNdpJl9j@bl50pXnBZm4Om3J1/ilbJNXDCjH3xfqaPMa0wO6fyuVS@Y5tgBicaMxDjnVxMRI3El2MfFBQsEz3B4ZABGY9MI7uBYn28xwVg/BK06UL6DnWrH981dbdjC2MO2Rfb7NfNfW6hDHD55WxFlskMQawe2DiPYlBqxNFZvs5E@OuJzTtxoV4s83Ncl@lvEt0ga6Gz5LnQQKlc34v4iQpRSTqcv "Python 3 – 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.
[Answer]
# [C (gcc)](https://gcc.gnu.org/), 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;
    }
}
void add_a_into_b(int size) {
    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;
            }
        }
    }
    add_a_into_b(size_x * size_y);
    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;
}
```
TIO link:
[Try it online!](https://tio.run/##zVdLc5swEL7zK7bJpAPGbuzkSOJbb5m2h/bQSTOMeBllsGWDiENS//W6K0AgYUgybg7VwZZ3v12t9in7k4Xv7/endOUneRDCVcYDyj7Fc0MjJdTTaUvCY0ExToMwoqsQsi1Zm2QM3hi@//z22YJfBuAKGDxDd1UssQQU@HIN10AcnUWQ5jm9Uh6yUMjpsHawjWkSgjm1DCNhqwWUH3TFRxAQTlyh8suPmxunl@s1XKRBRp9C9xFJl8rvovptPDAagM/WhUtcZDLXMyXGgufSmIilUBIpykwd/Loq@bizbQkqPVSefUvvEFdZiXunupCxq89KSLoIU1PKCb2rcOs2RtYbGy6cQ0AhAUUL0I2vxUc1zKowjc9M3V1WGpIkYb6p3wGDr9g0Us@vFLNIV2TV5/S5akBTv/dI5b2p9Jr0AIuiLOSlg2ZOl1q01KFYuY8H5zXQ@wp6L6Fo2r0OVcwzKXpeWmPpvrFRR8Ms8B5Gt1xkhjTxETJ1irQX3rUx894cM@/YmOn5M6CjxjZZ2sIURqEyCkekfBmKzM048ZLQfK2i/j1KNAJTJpLuZJjDZRctVhryPF3JjNMD0Q1JjZ05TTU/hSnrxAe4qHA@hqP7SKWgUwnyRBIE79qq7Bd6Vb5GVhu28rIy2fobTWlQzPIkcBcp27YX@J9CLiTwfGwa8Pu3UF/vKMyb/juBi4o5bzoukvq0lemvXXnmHIB2h0bgBTc5SXmYKuOie4G57PI9ragLHtTyUUy6rpLz84Fu1Masp3Oh3yYwszoSmEP1VV4StI8S1K0Tpx8jZr8k1l/sWqH1Z3v5UpITUxZGp1u/fUiLtFQSSc01@WhQClTkz5LQVVkpyPbH4MckHYn9w@2d@r6IQ7qIuXiYcUbNkj@7UwsWrRBzRgyFbJNy0wxYjv3aqgQt4T3FRIE@w/K4xkRRjRQM21Ynt/ao6QyKljT4OFlWc653eEnn9T91vPfS1rY8Mn4VMtgV3xR9DIR0FZyrr7/igFrX@aOa4@K9UUe6Fd20xVw/qD8o01gN3kZGTizZ99VYDjfwxoqDTr5OUSAyT86Ck3GVXs2o0XVHSZ7FJv47YTmvOe1c3u33@9lUrD9@lJBFtp8ky/3ka0Qy/hc "C (gcc) – Try It Online")
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 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!](https://tio.run/##NY5Na8MwDIbv@xWClF3mfTRl7DRIF7ZTYaMe9GA8MLXaCGY7cRToVvLbMzvbBEIfPHpfOcMNOsO0N9O0I27UuadvhEcooUb6JH9UsotcFbflzYOGK1iNAi4AYBeJUXLMxOJ14HZgAUWh4RKqaibghU5o3wJ5VvMMT8Rrb1UhYJWlfpdJavbtkmkCZEMH3tKx4cyV@s9uDjk4tQ1sGDd4YNUJsFrA2aZcCrhOmZ9PXS7jqNMz/9c1esa4jtF8Zd0MfJQ6E9V7eD61EfueglcLuY/Uch2cM95uyKNSmZumaXl/l@IH "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# [MATLAB](https://www.gnu.org/software/octave/), 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!](https://tio.run/##hY7LDoMgEEX3fMXsBEsqWlZS/RHjgviIJhSaam266a9TlBLddRbzuDl3ck0zy6Wztn/qZh6NhvuougkKmKRu8UAQCsIgEHoN7gCp39irJVxoJJWKCAJXWwt8Y/SS4YpBmnBgwo@1C/hpNeVxr4x5@G8JJwROu4@tLof6wVb@ZlrPUk7OsV/L4r9vA4Ph6nixR8UuP/74rxmhOXHZq/oI5PSApCQAnW6t/QI "Octave – 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.
[Answer]
# [C (clang)](http://clang.llvm.org/), 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) {
        fprintf(stderr, "Please, only positive numbers\n");
        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 contmask = 0;
                u8 carry = 0;
                for (int j = 0; j < lw; j++) {
                    contmask|= mainbuf[p+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 (contmask >= 4) {
                    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!](https://tio.run/##tVfPb@s2DL73r@DesMGOkyYu3qFAfpy2wy7DDgPeoQsKxZZjtbJkyHKTvL787RllyY7tOFl7mFHAKUnx@0hJJB1NIk7E9nT6mYmIlzGFRaFjJu/T1V1HxNmmL1NMbCtZTBMmKJSPUIqCbQWNIUqJurtjQkNGmPDMD6K20bhSjMzvt6e1D@93gA9LwDNa@GkJD7XQPEmOGDrxEJ8qNYYvf3FKCgq5km8MSTBcJ2KQgh@A@f@IL/68WauoLpWASWhFR4uEPFLKtqmGJRAtmVcxCdduoWHi9AsIP0BlbMFzWTDN3iiIMttQVfwHlYZLRvYZMplZxS5lnDYEVisIAmPgvx/ndsl0Ct9SoukbVff39/B3SjEZL2WBQREE1xI2FDZsC1TIcps2MAX7ThHGOQ4gnM0cpNHuWKzTS3WlLx9HoLN8UyZokBHOZeRV9i48ozc7bA0ia2Dgxpi@M@s/S84h16rAnafwnSo5hpgSjtgIzXTtawToqLJrnD0/b0rGNRPPucwjWQrtEuTDBMJxFZtMPFzs@@eYnBsW7216K4Uj@jRbozCcn7eBiZjumwycvagK77xBRrZVeORwW6yL1rZ1jkuVdV0myXNK8pwKvCrnBbUBt25qj2el89hz0PZvnkv/s3nHwGDkw2I@LBZ0r58drZ7aceJ4MYxBn0zjeIdLuTlE8ws9FohICp2R4vXSf21AlDoMaxOpoKokL5UeX0hlh@8gGGJjnhrux7LZ/Dx4Wc8HjRH9FbqGMIWvw8Ztq1@WH7EKlja4YUtTeV5M2mCJ@dtdC@hzTo9AOVaIj7gKJ3xn/b3OP4psC8PTtXwebwb68JlAP8TOekavs1tOrzl@GN1wffxEgE1W4Cpby3R1m2gr2ZPwRui30ozJCP8/jPquDiy6XJAHZrsvLXmwHKoVhn9TLFZ4wa6FYQx5dZauVKUGCPXDWNcjbNfDK0VtoBL/GETpInT/q4IwMKumG2AorVZjlPPuGkeroXhGPHa6DJ5HNLMt7tdOC0pMOUXx1I5dplH/JssNlnnseowDx@GI6gKnmx1VIBPQO4nDDIlSGneo4MI/UJuyAvBPUBrjCMgESigkJQ4/OA5k5JUCwYqcHy5CR4pDO@ea@NO5mQeBuVUZzaL84HXGkXF9lsfQnk@G0910ElZ1ErYwQ8Q8CJj/Xt8Ith5hVlopd7XU0HVjwaI1Zbj0fUupiZpok4dYCoybxDig8hgHH0UOBR4VorTp1omSGYSfp4VHuM6L5bFGYe9oWEX3sNe9wND8fe9GribTuKuv/QPk5tQ@q07cc4bN9z1RlHo1K4bDtHXQFrv9qES2QtbjYePftn18LZpLMHrAoRkO3QZf2eatwQkHwYMZBrtHG00WWGHfD/loiV7OGWoA7ZSDrz7g/nKiqOx7oPsuaA28d8D7PnAzVo6jwVUrDGvoHiiE3eeXFSVC@aEnv9ry1YBx7aTvfICyRDM18lQQ@tOHrrWZ2tLW4KRkEPVGgrzU5sPPQ7sAvj76Q7XKfCbQHccvyaLaohw/mzR@WpnvrpZR7Sqc@e0vu@pgOQaoOJ5OJ/yM@Rc "C (clang) – 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) {
        fprintf(stderr, "Please, only positive numbers\n");
        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 contmask = 0;
                u8 carry = 0;
                for (int j = 0; j < lw; j++) {
                    contmask|= mainbuf[p+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 (contmask >= 4) {
                    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 contmask = 0;
                u8 carry = 0;
                for (int j = 0; j < 2; j++) {
                    contmask|= mainbuf[p+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++) {
                    contmask|= mainbuf[p+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++) {
                    contmask|= mainbuf[p+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 (contmask >= 4) {
                    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);
}
```
```
[Answer]
# [Java (JDK)](http://jdk.java.net/), 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!](https://tio.run/##tVRNc4IwEL37K/bSDpSR@tGeqL12OlNPHh0OUSKGhmCTWPsx/na7gSAR61QPZRgCu@893m4@MvJOulnyutuxfFVIDRkGwrVmPJzMiRBURh03w4lIwzHRy6izWs84m8OcE6VgTJiA7w7YoNJE4/BesARyTHkTLZlIpzEQmSrfIMHqg7LjCATd1FFv8qk0zUMm/AixTGjIWZJwOmFfFKGWFAr6oZ@F9vYoxRL6QkWql4gawA14GPXBM55D9Sa11@j4cAsDHwIYGnYlMI3RZCrRd@XHhBrN2H2PHAZ6XT2dyULaopClMWBI6EU4PDjO8TsIqiY1yKxCZi1ktkcCWyASRqO6ZKcVZZ3X10bkVLpWKYufsniaxfjHpllRla0LPYnYAuWK/irWO6VhE9uOfWxNjzZLxil4Wq5puxd/du2Sth3bcUxfwZ01XXX3EXp@mxDULa143X4ZvW0x3f9DF/p/yQS/yWTnGMAYejhiXmoAY4Ers58dmBUFp0SYXT7jZjeaOYr@cYYOVtGha6dEB/UIQ7@xtyC4IttFGEaF8GEmKXmN6oX3L5vTnmfFWocrPAq147Y8u@y6x3u72/Xve3j9AA "Java (JDK) – 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.*
[Answer]
# [Rust](https://www.rust-lang.org/), 1903 bytes, score ~1000000 on TIO, ~3500000 on my computer
A simple port of [@ovs's Python answer](https://codegolf.stackexchange.com/a/242442/9288), 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](https://codegolf.stackexchange.com/a/242695/9288): 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!](https://tio.run/##pVVNc9owEL3zKxYmw8gB3GA3U0YpXHvsrReGYRRbBFEsHNlOnCb89ZKV5G9I25nqgqXdfe/taleoLElPp4gF6rBW2Z4nfUiYDGOx5/DaA1zk6kGJkIqQy3QMV4n4xavN/SGTIVMvlOexcmC@gFcbpNeepxBlKTxmLOUK5mCA3GB/kJw4d5Xf81aTWaNAT@K4TL6Qt2H@BjksahYHanC9LDACcqaagA0bz1MuQ9ICj1hcgS/AczCyFVo7r1E@BmwOas1ZsCVvGHSdw3AOfofOxCxvVjCaw2e4LuiX01XbDZFAgJAwdV1TSZjAtJNV6bfTfgJGf/Ot6QUSW88R7IyUQgYxME7LfBFHr1EZ1YQjO4PwD1HINflvLo3QKZ1ex17vg8xNtl5deJ2wwVs1OTCF1h5pLrCcF9PAX0yiQXmeuA5ri2okaa@kQ3/snX8ZOWZ3PN71sAQbWY2o7xFJIfO9MZi5hEz/ODBZwA8efEXDouiYchg1GI7iEw/6y5s7E4WqTKWslLKR5yCLuSjJ@kTbLNMYvNtbp6NmZhx8j8JwidSrD0TNPtCERkqfRbpdByxmgUhfCjx3zyVxnFpfOdaF@XyuWQLZrBrsywn4pfxYCZkakxE@6@quxQpTrrmtGnKYvZVVTzZ2gTl3XeteT2w11X/wKamIGMMOX1T8cNl9QhyTk1W@65x03iKjVIY810rTkNIgiimNhCww7XVj@zWMLLfGNpKpTJ8MXo@DsW0Mg7tquNk@PRaljBiyNCsmSw1cPlHK1APqrmJdmW7JtLHHe1qnW7xs/dDmbsxUwik1bYz3e/iJTVA7418OD1Iy@CaeOER4HyCz6J6rQaFO0xf3RYjUxdr4nuMmjwrfdPgEnvulLiFWw6/DbF9p6Y05s13htL0aPjMytHEtR9tcw7rx8Px4Op2mN2b9DjZ79pCcJt/fAQ "Rust – Try It Online")
---
# [Rust](https://www.rust-lang.org/), 2044 bytes, score ~750000 on TIO, ~6500000 on my computer
This one is a port of [@Polichinelle's C answer](https://codegolf.stackexchange.com/a/242612/9288), and borrow the trick in [@Aiden's answer](https://codegolf.stackexchange.com/a/242695/9288): 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!](https://tio.run/##pVVNc9owEL3nVyxMxiM12AV7MmFM4dJDj731kmEYYUQQxcK15YQ04a/XXUn@AGdonKkutqy3u@/trtZpnqmiiFmU7hdpvuNZDzImV4nYcXi5Alzk@iEVq1CsuFQDuM7Eb15vlvtcrlj6HPJDklKYzuDFGum14wriXMEyX8MUjBcv2u0lJ3RSg542OpI9FIqnhHpMPpNX5/AKB5g1ISg0nvXK1CoMYx6HYfbEEuLoSMbNAJwyKoY5MzHH98M5siF4rN8cCCjc2O1Ib3tB22q9T0GAkDDyPKMeXBi1yDTeRe1dWO9vcHrdgA@fShhujNs5zGbgX8LX2FEnnFvjJmfA4wVthkGlD5nVQul/K/0g849kw@2cjVOf3TLSrdqYq9J1lYvmS7fql/CO@SEmkbQJ0sXCfWvxwRy8q9@WYgrDtuO3b8bS7I7HyRWGXst65AQ@kSHkgT8AM2cg1w8K7gx@8OgLHsxKLtVw0c4w7iOPevfDCZT9a5RaKs2dl@XFroL1iB0XGjwA//aWttiMDSDwQ3DuMfT8AqnxBU54GIZPQm0WEUtYJNRz6c/bcUkobfh5/KC4XFXH5RyMcayVc5BlkI8p/aeAoKKfpEIqc2SIj9u8G7LCpGtqs4YxzN7SanoASyvsdBCtXtCY7TuYKhQRA9hSfUuEx5YZoUaTZb5tfWl1p9HTI/2XY38AdduVl2Y7P0HbJjuWeYiZwDSfyJVaqv5tcPkYhix9wKC1rSfVhoxO9pjkhdpgpV51CbyEpRkPQ9ODWJz9T6xgA8b/H48U6X8TjxxiTCbIPF7ytF@y0@HLZBMitdJ14FMv@5UqpPgZfO@u0Y@ygsbMNoWmfnJJbEnpOeoEMyaOtTsD2s5wmq7B78eiKO5uh7j@ROsde8gK93vhflWYHq7cKMmnkilU9Rc "Rust – Try It Online")
[Answer]
# [F# 6](https://fsharp.org/), 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
    }
type Adjustment =
    {
        Coordinates: Coordinates
        Value: int
    }
module SandPile =
    let maxHeight = 4
    let applyAdjustments (adjustments: Adjustment seq) =
        adjustments
        |> PSeq.groupBy (fun x -> x.Coordinates)
        |> PSeq.map (fun (_, changes) ->
            changes |> Seq.reduce (fun total next -> { total with Value = total.Value + next.Value })
        )
    
    let calculateDeltas (adjustments: Adjustment array) =
        adjustments
        |> PSeq.collect (fun originalAdjustment ->
            if originalAdjustment.Value >= maxHeight then
                let distribution = originalAdjustment.Value / maxHeight
                seq {
                    yield { originalAdjustment with Value = (originalAdjustment.Value % maxHeight) }
                    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
                [ originalAdjustment ]
        )
        |> PSeq.filter (fun adjustment -> adjustment.Value <> 0)
    
    let isStable (adjustments: Adjustment array) =
        adjustments
        |> PSeq.tryFind (fun x -> x.Value >= maxHeight)
        |> Option.isNone
    
    let printOut (adjustments: Adjustment array) =
        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
```
[Answer]
# [Rust](https://www.rust-lang.org/), 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.pad();
        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> {
    fn pad(&mut self) {
        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 {
            (size & !round_mask) + padding
        };
        let mut data = Vec::with_capacity_in(
            aligned_size * (aligned_size + padding),
            AlignedAlloc::<SIMD_ALIGN>,
        );
        data.extend(
            std::iter::repeat_with::<T, _>(Default::default)
                .take(aligned_size * (aligned_size + padding)),
        );
        let mut new = Mat {
            data,
            row_size: aligned_size + padding,
        };
        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.pad();
        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](https://play.rust-lang.org/?version=nightly&mode=debug&edition=2021&gist=cf91931df67694bdeb6163bca46e7db8)
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)
```
]