CA.c

// Compile with:
// $ gcc -Ofast -march=native -fopenmp ./CA.c -o CA
// Output is stream of PPM frames that can be played with mpv:
// $ ./CA | mpv - --scale=nearest

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <time.h>

// Output parameters

// Downscaling factor
#define S 2

// Output dimensions
#define W (1080 / S)
#define H (1920 / S)

// Space reserved for PPM header
#define MAX_HEADER 32

// Skips every other frame, useful if CA has flickering behavior
#define NO_FLICKER 1

// Numbers that are EPS apart are considered equal
#define EPS 1e-15

// How many neighbors each cell has
#define NUM_NEIGH 9

// Relative positions of neighbors of a cell
const int shift[NUM_NEIGH][2] = {
    {0, 0},
    {-1, 1}, {-1, -1}, {1, -1}, {1, 1},
    {-1, 0}, { 0, -1}, {0,  1}, {1, 0},
};

// Cell values lie within [0, MAX_VAL]
const double MAX_VAL = 1;

// Constants defining behavior of the CA
const double A = 0.11;
const double B = 0.80;

// Simple absolute value function
static inline double fabs(const double x){
    if (x > 0) return x;
    return -x;
}

// Compute new cell value based on its neighbors and itself
static inline double rule(const uint16_t i, const uint16_t j, const double world[W][H]){
    double cen = world[i][j];
    double sum = 0;

    for(int n = 0; n < NUM_NEIGH; n ++){
        // Visit all neighbors, taking wrapping into account
        sum += world[(W + i + shift[n][0]) % W][(H + j + shift[n][1]) % H];
    }

    double hash = A * MAX_VAL * (NUM_NEIGH + 1) + B * ((2.0/NUM_NEIGH) * sum  - 2 * cen) - sum;
    double val = 0;

    double min = 1e15;
    double temp;

    // Find the neighbor with value closest to 'hash' calculated earlier
    for(int n = 0; n < 9; n ++){
        temp = world[(W + i + shift[n][0]) % W][(H + j + shift[n][1]) % H];
        double dist = fabs(temp - hash);
        // Choose the smaller neighbor value in case of tie
        if (dist < min || (dist - min <= EPS && temp < val)){
            min = dist;
            val = temp;
        }
    }

    // Cell 'steals' value from its neighbor
    return val;
}

int main(void){
    static double world[2][W][H];
    uint16_t i, j, flag = 0;

    // Set up output buffer with three colors per pixel + PPM header
    static char out_buff[W * H * 3 + MAX_HEADER];
    char *image;

    int header_len = snprintf(out_buff, MAX_HEADER, "P6\n%d %d\n255\n", H, W);
    int out_len = 3 * W * H + header_len;

    // Pointer to the actual pixel values
    image = out_buff + header_len;

    srand48(time(NULL));

    // Randomize the world
    for(i = 0; i < W; i ++)
        for(j = 0; j < H; j ++)
            world[0][i][j] = drand48();


    while(1) {
        // Parallelize computation with OpenMP, not required but improves
        // performance
        #pragma omp parallel for collapse(2) schedule(guided)
        for(i = 0; i < W; i ++){
            for(j = 0; j < H; j ++){
                uint64_t color;
                world[!flag][i][j] = rule(i, j, world[flag]);

                // Transform a double into a pseudorandom int value
                color = 18446744073709551615.0 * world[!flag][i][j];
                color ^= color << 31;
                color ^= color >> 45;

                color &= 0xffffff;

                // Update a 3-byte value within internal representation of output
                // with the color generated above.
                *((uint32_t *) &image[3 * (i * H + j)]) &= ~0xffffff;
                *((uint32_t *) &image[3 * (i * H + j)]) |= color;
            }
        }

        flag ^= 1;

        // Output the image
        if (!NO_FLICKER || flag)
            fwrite(out_buff, 1, out_len, stdout);
    }

    return 0;
}