life2.c

#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>

#ifdef __SIZEOF_INT128__
typedef __uint128_t block_t;
#else
#warning "128-bit integers not supported, using 64-bit blocks!"
typedef uint64_t block_t;
#endif

#define WIDTH (1920)
#define HEIGHT (1080)

#define NUM_STEPS (1 << 20)

#define OUTPUT_PBM

#define BLOCK_BITS (8ULL * sizeof(block_t))
#define MAX_STANDARD_BITS (8ULL * sizeof(uint64_t))
#define ROT_SHIFT (BLOCK_BITS - 1ULL)

#define BLOCKS ((WIDTH * HEIGHT) / BLOCK_BITS)
#define W_SHIFT (WIDTH / BLOCK_BITS)

// Xorshiro256++ with an arbitrary seed

static uint64_t s[4] = {7537, 1675226715, 35763, 17253};

static inline uint64_t b_rotl(const uint64_t x, int k) { return (x << k) | (x >> (64 - k)); }

uint64_t xrand(void) {
    const uint64_t result = b_rotl(s[0] + s[3], 23) + s[0];

    const uint64_t t = s[1] << 17;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = b_rotl(s[3], 45);

    return result;
}

#define AS_U64(w, x, y) (((uint64_t *)(&w[(x)]))[(y)])

// Return the value of the specified cell, if 'val' is 1 or 0 it also sets the value of the cell to 'val'
static inline int get_set_helper(block_t world[BLOCKS], int i, int j, int val) {
    uint64_t idx = j * WIDTH + i;
    uint64_t block, sub_block, bit;
    uint64_t temp;

    block = idx / BLOCK_BITS;
    bit = (idx) % BLOCK_BITS;

    sub_block = bit / MAX_STANDARD_BITS;
    bit %= MAX_STANDARD_BITS;

    temp = AS_U64(world, block, sub_block);

    if (val == 0 || val == 1) {
        AS_U64(world, block, sub_block) = temp & ~(1ULL << bit) | (((uint64_t)val) << bit);
    }

    return (temp >> bit) & 1;
}

// Wrapper functions for convenience
static inline int get_cell(block_t world[BLOCKS], int i, int j) { return get_set_helper(world, i, j, -1); }
static inline void set_cell(block_t world[BLOCKS], int i, int j, int val) { get_set_helper(world, i, j, val); }

void rotl(block_t dst[BLOCKS], block_t src[BLOCKS]) {
    for (int i = BLOCKS - 1; i > 0; i--)
        dst[i] = (src[i] << 1) | (src[i - 1] >> ROT_SHIFT);

    dst[0] = (src[0] << 1) | (src[BLOCKS - 1] >> ROT_SHIFT);
}

void rotr(block_t dst[BLOCKS], block_t src[BLOCKS]) {
    for (int i = 0; i < BLOCKS - 1; i++)
        dst[i] = (src[i] >> 1) | ((src[i + 1] & 1ULL) << ROT_SHIFT);

    dst[BLOCKS - 1] = (src[BLOCKS - 1] >> 1) | ((src[0] & 1ULL) << ROT_SHIFT);
}

int main(void) {
    block_t world[4][BLOCKS];

#define FLIP(w) ((w) ^ 2) // equiv. ((w) == 1 ? 3 : 1)
    int w = 1;

    for (int i = 0; i < WIDTH; i++) {
        for (int j = 0; j < HEIGHT; j++) {
            set_cell(world[w], i, j, xrand() & 1);
        }
    }

    for (int step = 0; step < NUM_STEPS; step++) {

        rotl(world[0], world[w]);
        rotr(world[2], world[w]);

        for (int i = 0; i < BLOCKS; i++) {

#define W(a, b) (world[(a) ? (a) + 1 : w][(i + BLOCKS + ((b)*W_SHIFT)) % BLOCKS])
#define Q(a, b, c, d)                                                                                                  \
    block_t a = W(-1, (d));                                                                                            \
    block_t b = W(0, (d));                                                                                             \
    block_t c = W(1, (d));

            Q(A, B, C, -1);
            Q(D, Y, E, 0);
            Q(F, G, H, 1);

            block_t P, Q, R, S;

            P = A ^ B ^ C ^ D;
            Q = E ^ F ^ G ^ H;

            R = P | Q;
            S = P ^ Q;

            R ^= ((A & B) | (C & D)) & ((E & F) | (G & H));
            R ^= ((A & B) ^ (C & D)) | ((E & F) ^ (G & H));
            R ^= ((A ^ B) & (C ^ D)) | ((E ^ F) & (G ^ H));
            S |= ((A | B) & (C | D)) & ((E | F) & (G | H));

            world[FLIP(w)][i] = (Y | S) & (R ^ S);
        }

        w = FLIP(w);

#ifdef OUTPUT_PBM
        printf("P1\n%d %d\n", WIDTH, HEIGHT);
        for (int j = 0; j < HEIGHT; j++)
            for (int i = 0; i < WIDTH; i++)
                putchar('0' + get_cell(world[w], i, j));
#endif
    }

#ifndef OUTPUT_PBM
    int population_count = 0;

    for (int j = 0; j < HEIGHT; j++)
        for (int i = 0; i < WIDTH; i++)
            population_count += get_cell(world[w], i, j);

    printf("After %d generations, population count is %d.", NUM_STEPS, population_count);
#endif

    return 0;
}