Integral Apollonian gasket in C

/* Arithmos: geo-integralApollonian
 * Compile: cc -std=c11 -O2 geo-integralApollonian.c -lm -o geo-integralApollonian
 * Run:     ./geo-integralApollonian
 * Output:  geo-integralApollonian.svg (open this file in a browser)
 * Optional output path: ./geo-integralApollonian 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;}

/* Edit MAX_BEND (30..1200), DEPTH (2..18). b is signed reciprocal radius.
   Store integer triples (b,u,v)=(b,b*x,b*y), postponing division until drawing. */
typedef struct {int64_t b,u,v;} Bend;
typedef struct {Bend q[4];int depth,last;} PackingState;
static void draw(void) {
    const int MAX_BEND=360,DEPTH=12,CAPACITY=20000;
    static Bend circles[20000];static PackingState queue[20000];
    Bend seed[4]={{-1,0,0},{2,-1,0},{2,1,0},{3,0,2}};
    int count=4,head=0,tail=1;
    memcpy(circles,seed,sizeof seed);memcpy(queue[0].q,seed,sizeof seed);
    queue[0].depth=0;queue[0].last=-1;
    while(head<tail) {
        PackingState state=queue[head++];if(state.depth>=DEPTH)continue;
        for(int i=0;i<4;i++) {
            if(i==state.last)continue;
            Bend fresh={-state.q[i].b,-state.q[i].u,-state.q[i].v};
            for(int j=0;j<4;j++)if(j!=i){fresh.b+=2*state.q[j].b;fresh.u+=2*state.q[j].u;fresh.v+=2*state.q[j].v;}
            /* Descartes reflection: b'_i=2(sum of the other bends)-b_i. */
            if(fresh.b<=0||fresh.b>MAX_BEND)continue;
            int seen=0;
            for(int j=0;j<count;j++)if(circles[j].b==fresh.b&&circles[j].u==fresh.u&&circles[j].v==fresh.v){seen=1;break;}
            if(seen)continue;
            if(count>=CAPACITY||tail>=CAPACITY){fprintf(stderr,"Packing capacity exceeded\n");return;}
            /* Exact tangency test for this new circle and its three neighbors. */
            for(int j=0;j<4;j++)if(j!=i){
                Bend z=state.q[j];int64_t dx=fresh.u*z.b-z.u*fresh.b,dy=fresh.v*z.b-z.v*fresh.b,s=fresh.b+z.b;
                if(dx*dx+dy*dy!=s*s){fprintf(stderr,"Tangency invariant failed\n");return;}
            }
            circles[count++]=fresh;queue[tail]=state;queue[tail].q[i]=fresh;
            queue[tail].depth=state.depth+1;queue[tail++].last=i;
        }
    }
    for(int j=0;j<count;j++) {
        Bend z=circles[j];double b=(double)z.b,x=500+280*z.u/b,y=350-280*z.v/b,r=280/fabs(b),t=b<0?0:log(b)/log(MAX_BEND);
        nt_circle(x,y,r,t,.9,b<0?1.6:1.05);
        if(z.b>0&&z.b<=20){char label[30];snprintf(label,sizeof label,"%lld",(long long)z.b);nt_text(x-4,y+4,label);}
    }
    nt_text(30,35,"Integral Apollonian packing: every circle carries an exact integer bend");
    nt_text(30,670,"Negative bend -1 is the enclosing circle. Smaller circles have larger positive bends.");
}

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] : "geo-integralApollonian.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;
}