Hilbert curve in 3D

hilbert3d.c

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <inttypes.h>
#include <assert.h>

enum dim : uint8_t {
    X, Y, Z, DIM_COUNT
};

enum pos : uint8_t {
    NEAR_BOTTOM_LEFT,
    NEAR_BOTTOM_RIGHT,
    NEAR_TOP_LEFT,
    NEAR_TOP_RIGHT,
    FAR_BOTTOM_LEFT,
    FAR_BOTTOM_RIGHT,
    FAR_TOP_LEFT,
    FAR_TOP_RIGHT,
    POS_COUNT
};

struct pattern {
    enum pos origin;
    enum dim dim;
};

static struct pattern next_pattern(struct pattern prev, enum pos pos, size_t* index) {
    struct table_elem {
        unsigned origin : 3;
        unsigned dim : 2;
        unsigned index : 3;
    };

    static const struct table_elem table[POS_COUNT][DIM_COUNT][POS_COUNT] = {
#define INV(x) INV_##x
#define INV_0 1
#define INV_1 0
#define CAT(x, y, z, w) x##y##z##w
#define POS(x, y, z) CAT(AT_, x, y, z)
#define AT_000 NEAR_BOTTOM_LEFT
#define AT_010 NEAR_BOTTOM_RIGHT
#define AT_001 NEAR_TOP_LEFT
#define AT_011 NEAR_TOP_RIGHT
#define AT_100 FAR_BOTTOM_LEFT
#define AT_110 FAR_BOTTOM_RIGHT
#define AT_101 FAR_TOP_LEFT
#define AT_111 FAR_TOP_RIGHT
#define ELEM(x, y, z) \
        [POS(x, y, z)] = { \
            [X] = { \
                [POS(     x,      y,      z)] = { POS(x, y, z), Y, 0 }, \
                [POS(     x, INV(y),      z)] = { POS(x, y, z), Z, 1 }, \
                [POS(     x, INV(y), INV(z))] = { POS(x, y, z), Z, 2 }, \
                [POS(     x,      y, INV(z))] = { POS(x, INV(y), INV(z)), X, 3 }, \
                [POS(INV(x),      y, INV(z))] = { POS(x, INV(y), INV(z)), X, 4 }, \
                [POS(INV(x), INV(y), INV(z))] = { POS(INV(x), y, INV(z)), Z, 5 }, \
                [POS(INV(x), INV(y),      z)] = { POS(INV(x), y, INV(z)), Z, 6 }, \
                [POS(INV(x),      y,      z)] = { POS(INV(x), INV(y), z), Y, 7 }, \
            }, \
            [Y] = { \
                [POS(     x,      y,      z)] = { POS(x, y, z), Z, 0 }, \
                [POS(     x,      y, INV(z))] = { POS(x, y, z), X, 1 }, \
                [POS(INV(x),      y, INV(z))] = { POS(x, y, z), X, 2 }, \
                [POS(INV(x),      y,      z)] = { POS(INV(x), y, INV(z)), Y, 3 }, \
                [POS(INV(x), INV(y),      z)] = { POS(INV(x), y, INV(z)), Y, 4 }, \
                [POS(INV(x), INV(y), INV(z))] = { POS(INV(x), INV(y), z), X, 5 }, \
                [POS(     x, INV(y), INV(z))] = { POS(INV(x), INV(y), z), X, 6 }, \
                [POS(     x, INV(y),      z)] = { POS(x, INV(y), INV(z)), Z, 7 }, \
            }, \
            [Z] = { \
                [POS(     x,      y,      z)] = { POS(x, y, z), X, 0 }, \
                [POS(INV(x),      y,      z)] = { POS(x, y, z), Y, 1 }, \
                [POS(INV(x), INV(y),      z)] = { POS(x, y, z), Y, 2 }, \
                [POS(     x, INV(y),      z)] = { POS(INV(x), INV(y), z), Z, 3 }, \
                [POS(     x, INV(y), INV(z))] = { POS(INV(x), INV(y), z), Z, 4 }, \
                [POS(INV(x), INV(y), INV(z))] = { POS(x, INV(y), INV(z)), Y, 5 }, \
                [POS(INV(x),      y, INV(z))] = { POS(x, INV(y), INV(z)), Y, 6 }, \
                [POS(     x,      y, INV(z))] = { POS(INV(x), y, INV(z)), X, 7 }, \
            } \
        }

        ELEM(0, 0, 0),
        ELEM(0, 1, 0),
        ELEM(1, 1, 0),
        ELEM(1, 0, 0),
        ELEM(0, 0, 1),
        ELEM(0, 1, 1),
        ELEM(1, 1, 1),
        ELEM(1, 0, 1)
    };

    struct table_elem elem = table[prev.origin][prev.dim][pos];
    *index = elem.index;
    return (struct pattern) { .origin = elem.origin, .dim = elem.dim };
}

static inline enum pos make_pos(bool is_right, bool is_top, bool is_far) {
    return (is_far ? 4 : 0) | (is_top ? 2 : 0) | (is_right ? 1 : 0);
}

static inline size_t hilbert(uint32_t x, uint32_t y, uint32_t z, int k, enum pos start_pos, enum dim start_dim) {
    struct pattern pattern = { start_pos, start_dim };
    size_t order = 0;
    for (int i = k - 1; i >= 0; --i) {
        const enum pos pos = make_pos((x >> i) & 1, (y >> i) & 1, (z >> i) & 1);
        size_t index;
        pattern = next_pattern(pattern, pos, &index);
        order += index * (((size_t)1) << (3 * i));
    }
    return order;
}

static inline bool find_neighbor(size_t w, size_t h, size_t d, size_t order[d][h][w], uint32_t* x, uint32_t* y, uint32_t* z) {
    const uint32_t min_k = *z < 1 ? *z : *z - 1;
    const uint32_t max_k = *z + 1 >= d ? *z : *z + 1;
    const uint32_t min_j = *y < 1 ? *y : *y - 1;
    const uint32_t max_j = *y + 1 >= h ? *y : *y + 1;
    const uint32_t min_i = *x < 1 ? *x : *x - 1;
    const uint32_t max_i = *x + 1 >= w ? *x : *x + 1;

    for (uint32_t k = min_k; k <= max_k; k++) {
        for (uint32_t j = min_j; j <= max_j; j++) {
            for (uint32_t i = min_i; i <= max_i; i++) {
                if (i == *x && j == *y && k == *z)
                    continue;
                if (order[k][j][i] == order[*z][*y][*x] + 1) {
                    *z = k;
                    *y = j;
                    *x = i;
                    return true;
                }
            }
        }
    }
    return false;
}

int main() {
    static const int k = 3;
    static const uint32_t w = UINT32_C(1) << k;
    static const uint32_t h = UINT32_C(1) << k;
    static const uint32_t d = UINT32_C(1) << k;
    static const enum pos start_pos = NEAR_TOP_RIGHT;
    static const enum dim start_dim = X;

    size_t order[d][h][w];
    for (uint32_t z = 0; z < d; ++z) {
        for (uint32_t y = 0; y < h; ++y) {
            for (uint32_t x = 0; x < w; ++x)
                order[z][y][x] = hilbert(x, y, z, k, start_pos, start_dim);
        }
    }

    bool seen[w * h * d] = {};
    for (uint32_t z = 0; z < d; ++z) {
        for (uint32_t y = 0; y < h; ++y) {
            for (uint32_t x = 0; x < w; ++x) {
                size_t o = order[z][y][x];
                if (seen[o]) {
                    printf("element %zu already present\n", o);
                    return 1;
                }
                seen[o] = true;

                uint32_t i = x, j = y, k = z;
                if (o != w * h * d - 1 && !find_neighbor(w, h, d, order, &i, &j, &k)) {
                    printf("element at %"PRIu32", %"PRIu32", %"PRIu32" is invalid\n", x, y, z);
                    return 1;
                }
            }
        }
    }

    for (uint32_t z = 0; z < d; ++z) {
        for (uint32_t y = 0; y < h; ++y) {
            for (uint32_t x = 0; x < w; ++x)
                printf("%"PRIu32", %"PRIu32", %"PRIu32": %zu\n", x, y, z, order[z][y][x]);
        }
    }
    return 0;
}