Initial commit.

.gitignore

@@ -0,0 +1,10 @@
+sandpile
+sandpile-video
+sandpile_1d
+sandpile_add
+sandpile_fifo
+sandpile_pointers
+sandpile_quadran
+sandpile_quarter
+sandpile_quarter_fifo
+sandpile_struct

README.md

@@ -0,0 +1,80 @@
+# Searching for a fast Abelian Sandpile implementation
+
+For context about the Abelian Sandpile Model see: https://en.wikipedia.org/wiki/Abelian_sandpile_model
+
+I'm focusing only on the "flattening" step which I like to call "apply_gravity" in my code.
+
+I'm focusing only in flattening a single huge pile of sand placed in the middle.
+
+Here the performance I'm getting on an intel `i9-9980HK`:
+
+## sandpile_quarter.c
+
+This is my fastest implemtation, it relies on the 8-fold symetry of
+flattening a centered pile. I in fact only rely on a 4-fold symmetry
+for simplicity, so it could be enchanced again.
+
+> Mean +- std dev: 208 ms +- 5 ms
+
+## sandpile_1d.c
+
+This one does not rely on pointers of pointers to represent a 2d
+array. So no `[x][y]` access.  Instead it represents the surface as an
+1d line, in which, if the current cell is `i`:
+
+- the previous cell is at `surface[i - 1]`
+- the cell in the next row is at `surface[i + width]`
+- the cell in the previous row is at `surface[i - width]`
+- the next cell is at `surface[i + 1]`
+
+> Mean +- std dev: 731 ms +- 9 ms
+
+
+## sandpile.c
+
+This is the "dumb" implementation.
+
+> Mean +- std dev: 750 ms +- 10 ms
+
+
+## sandpile_struct.c
+
+This one is just to measure the cost of using a struct to store together the width and the values.
+
+> Mean +- std dev: 936 ms +- 16 ms
+
+
+## sandpile_quarter_fifo.c
+
+This is a merge of the `fifo` implementation and the `quarter` one.
+
+> Mean +- std dev: 831 ms +- 14 ms
+
+
+## sandpile_pointers.c
+
+This one stores, for each cells, 4 pointers to the west, south, north,
+and east cells so incrementing them is easy.
+
+> Mean +- std dev: 933 ms +- 27 ms
+
+
+## sandpile_add.c
+
+This one starts with a small (1024 typically) power of two stack of
+sand grains, applies gravity on it, and multiply the resulting grid by
+two. This is done iteratively until the needed number of grains is
+reached.
+
+The resulting grid is the same, but the performance are not better.
+
+> Mean +- std dev: 1.18 sec +- 0.03 sec
+
+
+## sandpile_fifo.c
+
+This one is the "smart" one, instead of doing a full scan of the grid,
+a FIFO of points to fire is kept up to date each time a cell is
+firered.
+
+> Mean +- std dev: 2.81 sec +- 0.05 sec

check.sh

@@ -0,0 +1,12 @@
+for file in sandpile.c sandpile_*.c
+do
+    echo "## $file"
+    dest=${file%.c}
+    cc -DNDEBUG -Ofast -W -Wall -pedantic --std=c99 "$file" -lm -lpng -o "$dest"
+    rm -f "$dest-100000.png"
+    python -m pyperf command --fast --quiet "./$dest" 100000
+    ./"$dest" 16 --stdout
+    ./"$dest" 100000 "$dest-100000.png"
+    echo
+done
+sha1sum ./*-100000.png

common.c

@@ -0,0 +1,106 @@
+#include <stdbool.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <png.h>
+
+typedef struct args {
+    int height;
+    char *png_filename;
+    bool print_to_stdout;
+} args;
+
+
+int parse_args(int argc, char **argv, args *args)
+{
+    if (argc < 2) {
+        fprintf(stderr, "Usage:\n");
+        fprintf(stderr, "    %s --stdout  # To print as an ASCII table.\n", argv[0]);
+        fprintf(stderr, "    %s file.png  # To write a PNG file.\n", argv[0]);
+        return 1;
+    }
+
+    if (sscanf(argv[1], "%i", &(args->height)) == EOF) {
+        fprintf(stderr, "Can't read height as an integer.\n");
+        return 1;
+    }
+
+    if (argc > 2) {
+        if (strcmp(argv[2], "--stdout") == 0)
+            args->print_to_stdout = true;
+        else
+            args->png_filename = argv[2];
+    }
+    return 0;
+}
+
+static void print_line(int width)
+{
+    printf("+");
+    for (int x = 0; x < width; x++)
+        printf("---+");
+    printf("\n");
+}
+
+
+void terrain_to_stdout(int width, int **terrain)
+{
+    for (int x = 0; x < width; x++) {
+        print_line(width);
+        printf("|");
+        for (int y = 0; y < width; y++) {
+            if (terrain[x][y] == 0)
+                printf("   |");
+            else if (terrain[x][y] < 10)
+                printf(" %i |", terrain[x][y]);
+            else
+                printf("%3i|", terrain[x][y]);
+        }
+        printf("\n");
+    }
+    print_line(width);
+}
+
+#define ONE_GRAIN "\x3E\x92\xCC"
+#define TWO_GRAINS "\x2A\x62\x8F"
+#define THREE_GRAINS "\x13\x29\x3D"
+#define CMAP "\xFF\xFF\xFF" ONE_GRAIN TWO_GRAINS THREE_GRAINS
+
+int terrain_to_png(int width, int **terrain, char *filename)
+{
+    png_image image;
+    png_byte *bytes;
+
+    bytes = malloc(sizeof(*bytes) * width * width);
+    memset(&image, 0, (sizeof image));
+    image.version = PNG_IMAGE_VERSION;
+    image.format = PNG_FORMAT_FLAG_COLORMAP|PNG_FORMAT_FLAG_COLOR;
+    image.width = width;
+    image.height = width;
+    image.colormap_entries = 4;
+    for (int x = 0; x < width; x++)
+        for (int y = 0; y < width; y++)
+            bytes[x * width + y] = terrain[x][y];
+    if (png_image_write_to_file(&image, filename, 0, bytes, width, CMAP))
+        return 0;
+
+    fprintf(stderr, "pngtopng: error: %s\n", image.message);
+    return 1;
+}
+
+void save_terrain(int width, int **terrain, args *args)
+{
+    if (args->png_filename)
+        terrain_to_png(width, terrain, args->png_filename);
+    if (args->print_to_stdout)
+        terrain_to_stdout(width, terrain);
+}
+
+int sandpile_width(int ngrains) {
+    int width;
+
+    width = pow(ngrains / 2 / 3.1415926535, 0.5) * 2 + 1;
+    width += 1 - width % 2;  // Ensure width is odd (terrain has a center)
+    return width;
+}

sandpile-video.c

@@ -0,0 +1,78 @@
+// Compile using:
+//  cc -W -Wall --pedantic --std=c99 -Ofast sandpile-video.c -lpng -lm -o sandpile-video
+// The run it to produce png files:
+//  sandpile-video 2000
+// Then use ffmpeg to convert to a video
+//  ffmpeg -pattern_type glob -i '*.png' out.mp4
+
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+
+void apply_gravity(int width, int **terrain)
+{
+    bool did_something;
+    int div;
+
+    while (1) {
+        did_something = false;
+
+        for (int x = 0; x < width; x++) {
+            for (int y = 0; y < width; y++) {
+                if (terrain[x][y] >= 4) {
+                    did_something = true;
+                    div = terrain[x][y] / 4;
+                    terrain[x][y] = terrain[x][y] % 4;
+                    terrain[x - 1][y] += div;
+                    terrain[x + 1][y] += div;
+                    terrain[x][y + 1] += div;
+                    terrain[x][y - 1] += div;
+                }
+            }
+        }
+        if (!did_something)
+            return;
+    }
+}
+
+
+int **sandpile_new(int width)
+{
+    int **terrain = calloc(width, sizeof(int*));
+    int *data = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+    return terrain;
+}
+
+#define STEP 100
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    int **terrain = sandpile_new(width);
+
+    char imgname[255] = {0};
+    for (int i = 0; i < args.height; i += STEP) {
+        terrain[width / 2][width / 2] += STEP;
+        apply_gravity(width, terrain);
+        sprintf(imgname, "sandpile-%i-%010i.png", args.height, i);
+        if (save_terrain(width, terrain, imgname) == 0)
+            printf("%i/%i saved\n", i, args.height);
+        else
+            printf("Can't save terrain\n");
+    }
+    return EXIT_SUCCESS;
+}

sandpile.c

@@ -0,0 +1,61 @@
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+
+void apply_gravity(int width, int **terrain)
+{
+    bool did_something;
+    int div;
+
+    while (1) {
+        did_something = false;
+        for (int x = 0; x < width; x++) {
+            for (int y = 0; y < width; y++) {
+                if (terrain[x][y] >= 4) {
+                    did_something = true;
+                    div = terrain[x][y] / 4;
+                    terrain[x][y] %= 4;
+                    terrain[x - 1][y] += div;
+                    terrain[x + 1][y] += div;
+                    terrain[x][y + 1] += div;
+                    terrain[x][y - 1] += div;
+                }
+            }
+        }
+        if (!did_something)
+            return;
+    }
+}
+
+
+int **sandpile_new(int width)
+{
+    int **terrain = calloc(width, sizeof(int*));
+    int *data = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+    return terrain;
+}
+
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    int **terrain = sandpile_new(width);
+    terrain[width / 2][width / 2] = args.height;
+    apply_gravity(width, terrain);
+    save_terrain(width, terrain, &args);
+    return EXIT_SUCCESS;
+}

sandpile_1d.c

@@ -0,0 +1,64 @@
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+
+void apply_gravity(int width, int *terrain)
+{
+    bool did_something;
+    int div;
+    int grid_size = width * width;
+
+    while (1) {
+        did_something = false;
+        for (int x = 0; x < grid_size; x++) {
+            if (terrain[x] >= 4) {
+                did_something = true;
+                div = terrain[x] / 4;
+                terrain[x] %= 4;
+                terrain[x - 1] += div;
+                terrain[x + 1] += div;
+                terrain[x - width] += div;
+                terrain[x + width] += div;
+            }
+        }
+        if (!did_something) {
+#ifdef DEBUG
+            fprintf(stderr, "looped %i times:\n", loops);
+            for (int i = 0; i < 9; i++) {
+                fprintf(stderr, "- %"PRIu64" times to break a pile of %i grains\n", by_height[i], i);
+            }
+            fprintf(stderr, "- %"PRIu64" times to break a pile of more than 13 grains\n", by_height[9]);
+            fprintf(stderr, "Total %i collapses\n", collapses);
+#endif
+            return;
+        }
+    }
+}
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    int *terrain = calloc(width * width, sizeof(int*));
+    terrain[width * (width / 2) + (width / 2)] = args.height;
+    apply_gravity(width, terrain);
+
+    // "convert" to 2D for printing:
+    int **terrain_2d = calloc(width, sizeof(int*));
+    for (int i = 0; i < width; i++) {
+        terrain_2d[i] = terrain + i * width;
+    }
+
+    save_terrain(width, terrain_2d, &args);
+    return EXIT_SUCCESS;
+}

sandpile_add.c

@@ -0,0 +1,83 @@
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+
+void sandpile_multiply(int width, int **terrain, int mul)
+{
+    for (int x = 0; x < width; x++)
+        for (int y = 0; y < width; y++)
+            terrain[x][y] *= mul;
+}
+
+void apply_gravity(int width, int **terrain)
+{
+    bool did_something;
+    int div;
+
+    while (1) {
+        did_something = false;
+        for (int x = 0; x < width; x++) {
+            for (int y = 0; y < width; y++) {
+                if (terrain[x][y] >= 4) {
+                    did_something = true;
+                    div = terrain[x][y] / 4;
+                    terrain[x][y] %= 4;
+                    terrain[x - 1][y] += div;
+                    terrain[x + 1][y] += div;
+                    terrain[x][y + 1] += div;
+                    terrain[x][y - 1] += div;
+                }
+            }
+        }
+        if (!did_something)
+            return;
+    }
+}
+
+
+int **sandpile_new(int width)
+{
+    int **terrain = calloc(width, sizeof(int*));
+    int *data = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+    return terrain;
+}
+
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    int **terrain = sandpile_new(width);
+    terrain[width / 2][width / 2] = args.height;
+    int target_power_of_two = (int)log2(args.height);
+    int current_height = args.height > 1024 ? 1024 : args.height;
+    terrain[width / 2][width / 2] = current_height;
+
+    apply_gravity(width, terrain);
+    for (int power_of_two = 10; power_of_two < target_power_of_two; power_of_two++) {
+        current_height *= 2;
+        sandpile_multiply(width, terrain, 2);
+        apply_gravity(width, terrain);
+    }
+
+    if (current_height != args.height) {
+        terrain[width / 2][width / 2] += args.height - current_height;
+        apply_gravity(width, terrain);
+    }
+
+    save_terrain(width, terrain, &args);
+    return EXIT_SUCCESS;
+}

sandpile_fifo.c

@@ -0,0 +1,126 @@
+#include <assert.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+
+typedef struct point {
+    int x;
+    int y;
+} point_t;
+
+
+typedef struct todo {
+    int width;
+    int height;
+    int size;
+    int idx;
+    int qty;
+    point_t *points;
+    bool *dedup;  // Used to avoid adding twice the same point, stored as x*width+y
+} todo_t;
+
+
+// The todo is deduplicated list of points.
+todo_t *todo_new(size_t width, size_t height)
+{
+    todo_t *todo = malloc(sizeof(*todo));
+
+    todo->qty = 0;
+    todo->width = width;
+    todo->height = height;
+    todo->size = width * 4;
+    todo->dedup = calloc(todo->width * todo->height, sizeof(bool));
+    todo->points = malloc(todo->size * sizeof(point_t));
+    return todo;
+}
+
+static void todo_write(todo_t *todo, point_t point)
+{
+    if (todo->dedup[point.x * todo->width + point.y])
+        return;
+    if (todo->qty == todo->size - 1) {
+        todo->points = realloc(todo->points, sizeof(point_t) * todo->size * 2);
+        todo->size *= 4;
+    }
+    todo->dedup[point.x * todo->width + point.y] = true;
+    todo->points[todo->qty] = point;
+    todo->qty += 1;
+}
+
+static point_t todo_read(todo_t *todo)
+{
+    assert(todo->qty > 0);
+    todo->qty -= 1;
+    todo->dedup[todo->points[todo->qty].x * todo->width + todo->points[todo->qty].y] = false;
+    return todo->points[todo->qty];
+}
+
+void unstable_to_todo(int width, int **terrain, todo_t *todo)
+{
+    for (int x = 0; x < width; x++)
+        for (int y = 0; y < width; y++)
+            if (terrain[x][y] >= 4)
+                todo_write(todo, (point_t){x, y});
+}
+
+void apply_gravity(int width, int **terrain)
+{
+    int div;
+    todo_t *todo;
+    point_t point;
+    unsigned int collapses = 0;
+
+    todo = todo_new(width, width);
+    unstable_to_todo(width, terrain, todo);
+    while (todo->qty) {
+        point = todo_read(todo);
+        collapses += 1;
+        div = terrain[point.x][point.y] / 4;
+        terrain[point.x][point.y] %= 4;
+        terrain[point.x - 1][point.y] += div;
+        terrain[point.x + 1][point.y] += div;
+        terrain[point.x][point.y + 1] += div;
+        terrain[point.x][point.y - 1] += div;
+        if (terrain[point.x - 1][point.y] >= 4)
+            todo_write(todo, (point_t){point.x-1, point.y});
+        if (terrain[point.x + 1][point.y] >= 4)
+            todo_write(todo, (point_t){point.x+1, point.y});
+        if (terrain[point.x][point.y + 1] >= 4)
+            todo_write(todo, (point_t){point.x, point.y+1});
+        if (terrain[point.x][point.y - 1] >= 4)
+            todo_write(todo, (point_t){point.x, point.y-1});
+    }
+}
+
+
+int **sandpile_new(int width)
+{
+    int **terrain = calloc(width, sizeof(int*));
+    int *data = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+    return terrain;
+}
+
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    int **terrain = sandpile_new(width);
+    terrain[width / 2][width / 2] = args.height;
+    apply_gravity(width, terrain);
+    save_terrain(width, terrain, &args);
+    return EXIT_SUCCESS;
+}

sandpile_pointers.c

@@ -0,0 +1,109 @@
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+typedef struct cell cell;
+
+struct cell {
+    cell *west;
+    cell *north;
+    cell *south;
+    cell *east;
+    int height;
+};
+
+
+static inline void fire(cell *cell) {
+    int div, rest;
+
+    div = cell->height / 4;
+    rest = cell->height % 4;
+    cell->height = rest;
+    cell->west->height += div;
+    cell->south->height += div;
+    cell->north->height += div;
+    cell->east->height += div;
+}
+
+
+void apply_gravity(int width, cell *cells)
+{
+    bool did_something;
+
+    while (1) {
+        did_something = false;
+        for (int x = 0; x < width*width; x++) {
+            if (cells[x].height >= 4) {
+                did_something = true;
+                fire(cells+x);
+            }
+        }
+        if (!did_something)
+            return;
+    }
+}
+
+cell *sandpile_new(int width, int height)
+{
+    cell **terrain = calloc(width, sizeof(cell*));
+    cell *data = calloc(width * width, sizeof(cell));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+    terrain[width / 2][width / 2].height = height;
+    for (int x = 0; x < width; x++) {
+        for (int y = 0; y < width; y++) {
+            terrain[x][y].west = &terrain[x-1][y];
+            terrain[x][y].south = &terrain[x][y-1];
+            terrain[x][y].north = &terrain[x][y+1];
+            terrain[x][y].east = &terrain[x+1][y];
+        }
+    }
+    return data;
+}
+
+int **cells_to_array(int width, cell *cells) {
+    // Create terrain:
+    int **terrain = calloc(width, sizeof(int*));
+    int *data = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+
+    // Organize cells in a 2d structure:
+    cell **cell_array = calloc(width, sizeof(cell*));
+    for (int i = 0; i < width; i++) {
+        cell_array[i] = cells + i * width;
+    }
+
+    // Copy cells to terrain:
+    for (int x = 0; x < width; x++)
+        for (int y = 0; y < width; y++)
+            terrain[x][y] = cell_array[x][y].height;
+    return terrain;
+}
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    cell *cells = sandpile_new(width, args.height);
+    apply_gravity(width, cells);
+
+    if (args.png_filename != NULL || args.print_to_stdout) {
+        // Convert to int** for easy display:
+        int **terrain = cells_to_array(width, cells);
+        save_terrain(width, terrain, &args);
+    }
+
+    return EXIT_SUCCESS;
+}

sandpile_quarter.c

@@ -0,0 +1,106 @@
+/*
+Some performance measures:
+
++--------+--------+
+| height | time   |
++--------+--------+
+| 2**16  | 0.1s   |
+| 2**17  | 0.4s   |
+| 2**18  | 1.5s   |
+| 2**19  | 6.5s   |
+| 2**20  | 23s    |
+| 2**21  | 1m36s  |
+| 2**22  | 6m18s  |
+| 2**23  | 25m4s  |
+| 2**24  | 96m54s |
++--------+--------+
+
+*/
+
+
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+typedef struct terrain {
+    int width;
+    int **cells;
+} terrain;
+
+
+void propagate_to_quartiles(terrain *quartile, terrain *full)
+{
+    for (int x = 0; x <= full->width / 2; x++)
+        for (int y = 0; y <= full->width / 2; y++) {
+            full->cells[x][y] = quartile->cells[x][y];
+            full->cells[full->width-x-1][y] = quartile->cells[x][y];
+            full->cells[x][full->width-y-1] = quartile->cells[x][y];
+            full->cells[full->width-x-1][full->width-y-1] = quartile->cells[x][y];
+        }
+}
+
+void apply_gravity(terrain *land)
+{
+    bool did_something;
+    int div;
+
+    while (1) {
+        did_something = false;
+        for (int x = 0; x <= land->width - 2; x++) {
+            for (int y = 0; y <= land->width - 2; y++) {
+                if (land->cells[x][y] >= 4) {
+                    did_something = true;
+                    div = land->cells[x][y] / 4;
+                    land->cells[x][y] = land->cells[x][y] % 4;
+                    land->cells[x - 1][y] += div;
+                    land->cells[x + 1][y] += div;
+                    land->cells[x][y + 1] += div;
+                    land->cells[x][y - 1] += div;
+                    if (x == land->width - 3)
+                        land->cells[x + 1][y] += div;
+                    if (y == land->width - 3)
+                        land->cells[x][y + 1] += div;
+                }
+            }
+        }
+        if (!did_something)
+            return;
+    }
+}
+
+terrain *terrain_new(int width)
+{
+    terrain *land = malloc(sizeof(terrain));
+    int **lines = calloc(width, sizeof(int*));
+    int *cells = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        lines[i] = cells + i * width;
+    }
+    land->width = width;
+    land->cells = lines;
+    return land;
+}
+
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+
+    terrain *land = terrain_new(width);
+    terrain *small_land = terrain_new(width / 2 + 2);
+    small_land->cells[width / 2][width / 2] = args.height;
+    apply_gravity(small_land);
+    propagate_to_quartiles(small_land, land);
+    save_terrain(land->width, land->cells, &args);
+    return EXIT_SUCCESS;
+}

sandpile_quarter_fifo.c

@@ -0,0 +1,143 @@
+#include <assert.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+typedef struct terrain {
+    int width;
+    int **cells;
+} terrain;
+
+
+void propagate_to_quartiles(terrain *quartile, terrain *full)
+{
+    for (int x = 0; x <= full->width / 2; x++)
+        for (int y = 0; y <= full->width / 2; y++) {
+            full->cells[x][y] = quartile->cells[x][y];
+            full->cells[full->width-x-1][y] = quartile->cells[x][y];
+            full->cells[x][full->width-y-1] = quartile->cells[x][y];
+            full->cells[full->width-x-1][full->width-y-1] = quartile->cells[x][y];
+        }
+}
+
+typedef struct point {
+    int x;
+    int y;
+} point_t;
+
+
+typedef struct todo {
+    int width;
+    int height;
+    int size;
+    int idx;
+    int qty;
+    point_t *points;
+    bool *dedup;  // Used to avoid adding twice the same point, stored as x*width+y
+} todo_t;
+
+
+// The todo is deduplicated list of points.
+todo_t *todo_new(size_t width, size_t height)
+{
+    todo_t *todo = malloc(sizeof(*todo));
+
+    todo->qty = 0;
+    todo->width = width;
+    todo->height = height;
+    todo->size = width * 4;
+    todo->dedup = calloc(todo->width * todo->height, sizeof(bool));
+    todo->points = malloc(todo->size * sizeof(point_t));
+    return todo;
+}
+
+static void todo_write(todo_t *todo, point_t point)
+{
+    if (todo->dedup[point.x * todo->width + point.y])
+        return;
+    if (todo->qty == todo->size - 1) {
+        todo->points = realloc(todo->points, sizeof(point_t) * todo->size * 2);
+        todo->size *= 4;
+    }
+    todo->dedup[point.x * todo->width + point.y] = true;
+    todo->points[todo->qty] = point;
+    todo->qty += 1;
+}
+
+static point_t todo_read(todo_t *todo)
+{
+    assert(todo->qty > 0);
+    todo->qty -= 1;
+    todo->dedup[todo->points[todo->qty].x * todo->width + todo->points[todo->qty].y] = false;
+    return todo->points[todo->qty];
+}
+
+
+void apply_gravity(terrain *land, int middle_column, todo_t *todo)
+{
+    int div;
+    point_t point;
+
+    while (todo->qty) {
+        point = todo_read(todo);
+        div = land->cells[point.x][point.y] / 4;
+        land->cells[point.x][point.y] = land->cells[point.x][point.y] % 4;
+        land->cells[point.x - 1][point.y] += div;
+        if (point.x < middle_column)
+            land->cells[point.x + 1][point.y] += div;
+        if (point.y < middle_column)
+            land->cells[point.x][point.y + 1] += div;
+        land->cells[point.x][point.y - 1] += div;
+        if (point.x == middle_column - 1)
+            land->cells[point.x + 1][point.y] += div;
+        if (point.y == middle_column - 1)
+            land->cells[point.x][point.y + 1] += div;
+        if (land->cells[point.x - 1][point.y] > 3)
+            todo_write(todo, (point_t){point.x-1, point.y});
+        if (point.x < middle_column && land->cells[point.x + 1][point.y] > 3)
+            todo_write(todo, (point_t){point.x+1, point.y});
+        if (point.y < middle_column && land->cells[point.x][point.y + 1] > 3)
+            todo_write(todo, (point_t){point.x, point.y+1});
+        if (land->cells[point.x][point.y - 1] > 3)
+            todo_write(todo, (point_t){point.x-1, point.y-1});
+    }
+}
+
+terrain *terrain_new(int width)
+{
+    terrain *land = malloc(sizeof(terrain));
+    int **lines = calloc(width, sizeof(int*));
+    int *cells = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        lines[i] = cells + i * width;
+    }
+    land->width = width;
+    land->cells = lines;
+    return land;
+}
+
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+
+    terrain *land = terrain_new(width);
+    terrain *small_land = terrain_new(width / 2 + 1);
+    small_land->cells[width / 2][width / 2] = args.height;
+    todo_t *todo = todo_new(width, args.height);
+    todo_write(todo, (point_t){width/2, width/2});
+    apply_gravity(small_land, width/2, todo);
+    propagate_to_quartiles(small_land, land);
+    save_terrain(land->width, land->cells, &args);
+    return EXIT_SUCCESS;
+}

sandpile_struct.c

@@ -0,0 +1,70 @@
+#include <stdlib.h>
+#include <stdbool.h>
+#include <string.h>
+#include <inttypes.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "common.c"
+
+
+typedef struct terrain {
+    int width;
+    int **cells;
+} terrain;
+
+
+void apply_gravity(terrain *land)
+{
+    bool did_something;
+    int div;
+    int width = land->width;
+    int **terrain = land->cells;
+
+    while (1) {
+        did_something = false;
+        for (int x = 0; x < width; x++) {
+            for (int y = 0; y < width; y++) {
+                if (terrain[x][y] >= 4) {
+                    did_something = true;
+                    div = terrain[x][y] / 4;
+                    terrain[x][y] %= 4;
+                    terrain[x - 1][y] += div;
+                    terrain[x + 1][y] += div;
+                    terrain[x][y + 1] += div;
+                    terrain[x][y - 1] += div;
+                }
+            }
+        }
+        if (!did_something)
+            return;
+    }
+}
+
+terrain *sandpile_new(int width)
+{
+    terrain *land = malloc(sizeof(terrain));
+    int **terrain = calloc(width, sizeof(int*));
+    int *data = calloc(width * width, sizeof(int));
+    for (int i = 0; i < width; i++) {
+        terrain[i] = data + i * width;
+    }
+    land->width = width;
+    land->cells = terrain;
+    return land;
+}
+
+int main(int argc, char **argv)
+{
+    args args = {0};
+
+    if (parse_args(argc, argv, &args))
+        return EXIT_FAILURE;
+
+    int width = sandpile_width(args.height);
+    terrain *land = sandpile_new(width);
+    land->cells[width / 2][width / 2] = args.height;
+    apply_gravity(land);
+    save_terrain(land->width, land->cells, &args);
+    return EXIT_SUCCESS;
+}