Tribonacci Rauzy tile in C

/* Arithmos: deep-rauzy-tribonacci
 * Compile: cc -std=c11 -O2 deep-rauzy-tribonacci.c -lm -o deep-rauzy-tribonacci
 * Run:     ./deep-rauzy-tribonacci
 * Output:  deep-rauzy-tribonacci.svg (open this file in a browser)
 * Optional output path: ./deep-rauzy-tribonacci my-image.svg
 * Edit the constants in draw() to explore another case.
 */
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static FILE *nt_out;
#define NT_PI 3.14159265358979323846

/* The small SVG writer keeps this program free of graphics dependencies.
 * Coordinates are pixels on a 1000 x 700 drawing surface.
 * t runs from 0 to 1 through mint, blue, rose, and gold.
 */
static inline void nt_color(double t, char hex[8]) {
    const double stops[4][3] = {
        {91,227,201}, {128,146,240}, {218,138,220}, {244,200,127}
    };
    t = fmax(0.0, fmin(1.0,t)) * 3.0;
    int band = (int)fmin(2.0,floor(t));
    double blend = t - band;
    int r[3];
    for (int k=0;k<3;k++) r[k]=(int)lround(stops[band][k]*(1.0-blend)+stops[band+1][k]*blend);
    snprintf(hex,8,"#%02x%02x%02x",r[0],r[1],r[2]);
}
static inline void nt_line(double x,double y,double X,double Y,double t,double alpha,double width) {
    char color[8];nt_color(t,color);
    fprintf(nt_out,"<path d=\"M%.3f %.3f L%.3f %.3f\" fill=\"none\" stroke=\"%s\" stroke-opacity=\"%.3f\" stroke-width=\"%.3f\"/>\n",x,y,X,Y,color,alpha,width);
}
static inline void nt_dot(double x,double y,double r,double t,double alpha) {
    char color[8];nt_color(t,color);
    fprintf(nt_out,"<circle cx=\"%.3f\" cy=\"%.3f\" r=\"%.3f\" fill=\"%s\" fill-opacity=\"%.3f\"/>\n",x,y,r,color,alpha);
}
static inline void nt_circle(double x,double y,double r,double t,double alpha,double width) {
    char color[8];nt_color(t,color);
    fprintf(nt_out,"<circle cx=\"%.3f\" cy=\"%.3f\" r=\"%.3f\" fill=\"none\" stroke=\"%s\" stroke-opacity=\"%.3f\" stroke-width=\"%.3f\"/>\n",x,y,r,color,alpha,width);
}
static inline void nt_rect(double x,double y,double w,double h,double t,double alpha) {
    char color[8];nt_color(t,color);
    fprintf(nt_out,"<rect x=\"%.3f\" y=\"%.3f\" width=\"%.3f\" height=\"%.3f\" fill=\"%s\" fill-opacity=\"%.3f\"/>\n",x,y,w,h,color,alpha);
}
static inline void nt_text(double x,double y,const char *text) {
    fprintf(nt_out,"<text x=\"%.3f\" y=\"%.3f\" fill=\"#ededf3\" font-family=\"monospace\" font-size=\"16\">",x,y);
    for (;*text;text++) {
        if (*text=='&') fputs("&amp;",nt_out);
        else if (*text=='<') fputs("&lt;",nt_out);
        else if (*text=='>') fputs("&gt;",nt_out);
        else fputc(*text,nt_out);
    }
    fputs("</text>\n",nt_out);
}
static inline int nt_gcd(int a,int b) {a=abs(a);b=abs(b);while(b){int r=a%b;a=b;b=r;}return a;}
static inline int nt_prime(int n) {if(n<2)return 0;for(int d=2;d<=n/d;d++)if(n%d==0)return 0;return 1;}

/* Tribonacci Rauzy tile. Edit N for the number of nonempty prefixes.
   1->12, 2->13, 3->1. Project exact letter counts; color the FINAL letter. */
static void draw(void) {
    enum { N = 50000 };
    unsigned char *word = malloc(N), *next = malloc(N);
    double *px = malloc(N * sizeof *px), *py = malloc(N * sizeof *py);
    if (!word || !next || !px || !py) {
        free(word); free(next); free(px); free(py);
        nt_text(50, 80, "Not enough memory for this point count."); return;
    }
    int len = 1; word[0] = 1;
    while (len < N) {
        int m = 0;
        for (int i = 0; i < len && m < N; ++i) {
            next[m++] = 1;
            if (word[i] != 3 && m < N) next[m++] = word[i] + 1;
        }
        unsigned char *swap = word; word = next; next = swap; len = m;
    }
    double beta = 1.84;
    for (int j = 0; j < 12; ++j)
        beta -= (beta*beta*beta-beta*beta-beta-1)/(3*beta*beta-2*beta-1);
    double ar = (1-beta)/2, ai = -sqrt(1/beta-ar*ar), norm = ar*ar+ai*ai;
    double vx[3] = {1, ar-1, ar/norm}, vy[3] = {0, ai, -ai/norm};
    int counts[3] = {0,0,0};
    double xmin=1e30, xmax=-1e30, ymin=1e30, ymax=-1e30;
    for (int i=0; i<N; ++i) {
        ++counts[word[i]-1];
        px[i]=counts[0]*vx[0]+counts[1]*vx[1]+counts[2]*vx[2];
        py[i]=counts[0]*vy[0]+counts[1]*vy[1]+counts[2]*vy[2];
        xmin=fmin(xmin,px[i]); xmax=fmax(xmax,px[i]);
        ymin=fmin(ymin,py[i]); ymax=fmax(ymax,py[i]);
    }
    double scale=fmin(900/fmax(xmax-xmin,1e-12),580/fmax(ymax-ymin,1e-12));
    for (int i=0; i<N; ++i)
        nt_dot(500+(px[i]-(xmin+xmax)/2)*scale,
               340-(py[i]-(ymin+ymax)/2)*scale,1.3,(word[i]-1)/2.0,.72);
    nt_text(50,675,"Three colors: final prefix letters 1, 2, 3. Finite projected sample.");
    free(word); free(next); free(px); free(py);
}

int main(int argc, char **argv) {
    if (argc > 2) {
        fprintf(stderr, "Usage: %s [OUTPUT.svg]\n", argv[0]);
        return EXIT_FAILURE;
    }
    const char *filename = argc == 2 ? argv[1] : "deep-rauzy-tribonacci.svg";
    nt_out = fopen(filename, "wb");
    if (!nt_out) { perror(filename); return EXIT_FAILURE; }
    fputs("<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"1000\" height=\"700\" viewBox=\"0 0 1000 700\">\n"
          "<rect width=\"1000\" height=\"700\" fill=\"#171721\"/>\n", nt_out);
    draw();
    fputs("</svg>\n", nt_out);
    int failed = ferror(nt_out);
    if (fclose(nt_out) != 0) failed = 1;
    if (failed) { fputs("Could not finish writing the image.\n", stderr); return EXIT_FAILURE; }
    printf("Wrote %s\n", filename);
    return EXIT_SUCCESS;
}