pi fractal

Makefile

CC     = c99
CFLAGS = -Wall -Wextra -Ofast
LDLIBS = -lm
pi: pi.c

pi.c

/* https://redd.it/70vzlx */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

#include <getopt.h>

#define RULE(n, p) (0x1f & (rules[n] >> ((p) * 5)))
static const unsigned long long rules[] = {
    0x0b7c1861e72b,0x1def3bbcf35a,0x1cef75ab779c,0x15b771dee3b6,0x11d77bae63bc,
    0x021844240800,0x0cb419d920f7,0x011844110442,0x1dd86d2c59b3,0x02287875dd00,
    0x0aae66fdebb6,0x1d2b55508630,0x1decbbdbf7b7,0x18661f6d777b,0x13bf5dce778e,
    0x1d3f7073b0f0,0x0a6f6b5ed516,0x1005db815ee4,0x01c9448d8625,0x15e4fbd9bf39, 
    0x0200c46080a1,0x1cb7776b537b,0x0210c4320421,0x18d771cee37d,0x1de4f7cbef35,
    0x1b38f4daf39c,0x16761f67d2c5,0x0010ce110861,0x07b17a7b29e3,0x030844108c21
};

static int
color(long long x, long long y, long long div)
{
    int c = 0;
    for (int i = 1; div > 0; i++) {
        c = RULE(c, y / div * 3 + x / div);
        x %= div;
        y %= div;
        div /= 3;
    }
    return c;
}

int main(int argc, char **argv)
{
    int n = 4;
    int scale = 10;
    long long xmin = 0;
    long long ymin = 0;
    long long xsize = LLONG_MAX;
    long long ysize = LLONG_MAX;

    int r, option;
    while ((option = getopt(argc, argv, "n:r:s:")) != -1) {
        switch (option) {
            case 'n':
                n = atoi(optarg);
                break;
            case 'r':
                r = sscanf(optarg,"%lldx%lld+%lld+%lld",
                           &xsize, &ysize, &xmin, &ymin);
                switch (r) {
                    case 0:
                        xsize = LLONG_MAX;
                    case 1:
                        ysize = LLONG_MAX;
                    case 2:
                        xmin = 0;
                    case 3:
                        ymin = 0;
                }
                break;
            case 's':
                scale = atoi(optarg);
                break;
            default:
                exit(EXIT_FAILURE);
        }
    }

    long long div = 1;
    for (int i = 1; i < n; i++)
        div *= 3;

    if (scale > 0) {
        /* Scale up */
        long long size = 3 * div;
        if (xmin + xsize > size)
            xsize = size - xmin;
        if (ymin + ysize > size)
            ysize = size - ymin;
        printf("P2\n%lld %lld\n29\n", xsize * scale, ysize * scale);
        for (long long y = ymin; y < ymin + ysize; y++)
            for (int sy = 0; sy < scale; sy++)
                for (long long x = xmin; x < xmin + xsize; x++)
                    for (int sx = 0; sx < scale; sx++)
                        printf("%d\n", color(x, y, div));

    } else {
        /* Scale down */
        scale = -scale;
        long long orig = 3 * div;
        long long size = orig / scale;

        int ksize = 2 * scale + 1;
        double *kernel = malloc(sizeof(*kernel) * ksize * ksize);
        double pi = atan(1) * 4;
        double stdev = 1.0;
        for (int y = 0; y < ksize; y++) {
            double fy = (y - scale) / (stdev * scale);
            for (int x = 0; x < ksize; x++) {
                double fx = (x - scale) / (stdev * scale);
                int i = y * ksize + x;
                double e = exp((fx * fx + fy * fy) / (-2 * stdev * stdev));
                kernel[i] = e / (2 * pi * stdev);
            }
        }

        printf("P2\n%lld %lld\n255\n", size, size);
        for (long long y = 0; y < orig; y += scale) {
            for (long long x = 0; x < orig; x += scale) {
                double acc = 0;
                double sum = 0;
                for (long ky = 0; ky < ksize; ky++) {
                    int dy = ky - scale;
                    long long sy = y + dy;
                    for (long kx = 0; kx < ksize; kx++) {
                        int dx = kx - scale;
                        long long sx = x + dx;
                        double k = kernel[ky * ksize + kx];
                        if (sx >= 0 && sx < orig && sy >= 0 && sy < orig) {
                            double c = color(x + dx, y + dy, div) / 29.0;
                            sum += k * c;
                            acc += k;
                        }
                    }
                }
                printf("%d\n", (int)(sum / acc * 255));
            }
        }

        free(kernel);
    }
}