/* Arithmos: geo-fareyQuadraticRiver
 * Compile: cc -std=c11 -O2 geo-fareyQuadraticRiver.c -lm -o geo-fareyQuadraticRiver
 * Run:     ./geo-fareyQuadraticRiver
 * Output:  geo-fareyQuadraticRiver.svg (open this file in a browser)
 * Optional output path: ./geo-fareyQuadraticRiver 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 DEPTH (1..8), RADICAND (2..30); square inputs advance to the next integer.
   The bright river joins triangle centers; the smooth axis is a separate object. */
typedef struct {int p,q;} Rational;
typedef struct {Rational v[3];int depth,entry;double x,y;} FareyTriangle;
typedef struct {Rational a,b;int depth,left,right;} FareyEdge;
static Rational rational(int p,int q) {
    int g=nt_gcd(abs(p),abs(q));p/=g;q/=g;if(q<0||(!q&&p<0)){p=-p;q=-q;}
    Rational v={p,q};return v;
}
static int equal_rational(Rational a,Rational b){return a.p*b.q==b.p*a.q;}
static void ideal(Rational a,double *x,double *y){double d=(double)a.p*a.p+(double)a.q*a.q;*x=((double)a.p*a.p-(double)a.q*a.q)/d;*y=-2.0*a.p*a.q/d;}
static int quadratic(Rational a,int D){return a.p*a.p-D*a.q*a.q;}
static void geodesic(double ax,double ay,double bx,double by,int steps,double color,double alpha,double width) {
    double det=ax*by-ay*bx;
    if(fabs(det)<1e-12){nt_line(500+280*ax,350-280*ay,500+280*bx,350-280*by,color,alpha,width);return;}
    double cx=(by-ay)/det,cy=(ax-bx)/det,r=sqrt(fmax(0,cx*cx+cy*cy-1));
    double start=atan2(ay-cy,ax-cx),delta=atan2(by-cy,bx-cx)-start;
    while(delta>NT_PI)delta-=2*NT_PI;while(delta< -NT_PI)delta+=2*NT_PI;
    double px=ax,py=ay;
    for(int j=1;j<=steps;j++){double a=start+delta*j/steps,x=cx+r*cos(a),y=cy+r*sin(a);if(j==steps){x=bx;y=by;}nt_line(500+280*px,350-280*py,500+280*x,350-280*y,color,alpha,width);px=x;py=y;}
}
static void draw(void) {
    const int DEPTH=6,RADICAND=2;
    int D=RADICAND;if((int)sqrt(D)*(int)sqrt(D)==D)D++;
    static FareyTriangle tr[800];static FareyEdge ed[1600];int nt=1,ne=0;
    tr[0]=(FareyTriangle){{{0,1},{1,1},{1,0}},0,-1,0,0};
    for(int h=0;h<nt;h++) {
        double kx=0,ky=0;
        for(int k=0;k<3;k++){double x,y;ideal(tr[h].v[k],&x,&y);kx+=x/3;ky+=y/3;}
        double den=1+sqrt(fmax(0,1-kx*kx-ky*ky));tr[h].x=kx/den;tr[h].y=ky/den;
        for(int i=0;i<3;i++) {
            Rational a=tr[h].v[i],b=tr[h].v[(i+1)%3],op=tr[h].v[(i+2)%3];int e=-1;
            for(int j=0;j<ne;j++)if((equal_rational(ed[j].a,a)&&equal_rational(ed[j].b,b))||(equal_rational(ed[j].a,b)&&equal_rational(ed[j].b,a))){e=j;break;}
            if(e<0){if(ne>=1600){fprintf(stderr,"Farey edge capacity exceeded\n");return;}e=ne++;ed[e]=(FareyEdge){a,b,tr[h].depth,h,-1};}
            else ed[e].right=h;
            if(tr[h].depth>=DEPTH||e==tr[h].entry)continue;
            Rational sum=rational(a.p+b.p,a.q+b.q),dif=rational(a.p-b.p,a.q-b.q),next=equal_rational(sum,op)?dif:sum;
            if(nt>=800){fprintf(stderr,"Farey triangle capacity exceeded\n");return;}
            tr[nt++]=(FareyTriangle){{a,b,next},tr[h].depth+1,e,0,0};
        }
    }
    for(int e=0;e<ne;e++) {
        double ax,ay,bx,by;ideal(ed[e].a,&ax,&ay);ideal(ed[e].b,&bx,&by);
        int steps=(int)round(40/(1+ed[e].depth*.4));if(steps<5)steps=5;
        geodesic(ax,ay,bx,by,steps,.08+.55*ed[e].depth/DEPTH,.46,.85);
        nt_dot(500+280*ax,350-280*ay,1.35,quadratic(ed[e].a,D)>0?.2:.55,.8);
        nt_dot(500+280*bx,350-280*by,1.35,quadratic(ed[e].b,D)>0?.2:.55,.8);
    }
    for(int e=0;e<ne;e++)if(ed[e].right>=0&&(int64_t)quadratic(ed[e].a,D)*quadratic(ed[e].b,D)<0) {
        FareyTriangle a=tr[ed[e].left],b=tr[ed[e].right];
        nt_line(500+280*a.x,350-280*a.y,500+280*b.x,350-280*b.y,.95,1,3.2);
        nt_dot(500+280*a.x,350-280*a.y,2.4,.95,1);nt_dot(500+280*b.x,350-280*b.y,2.4,.95,1);
    }
    double x=(D-1.0)/(D+1),y=2*sqrt(D)/(D+1);
    geodesic(x,-y,x,y,160,.73,.8,1.5);nt_dot(500+280*x,350+280*y,4,.73,1);nt_dot(500+280*x,350-280*y,4,.73,1);
    nt_circle(500,350,280,.05,.65,1.2);
    nt_text(30,35,"Farey quadratic river: Q(p,q) = p^2 - Dq^2 separates positive and negative regions");
    nt_text(30,670,"Bright dual-tree links trace the river; the smooth axis ends at +sqrt(D), -sqrt(D).");
}

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-fareyQuadraticRiver.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;
}