arithmos
MAKE THE MATHEMATICS YOURSELF

Conway’s quadratic river in C

A complete program. Standard C, a compiler, and the rule behind the shape.

Output generated by this C program: Conway’s quadratic river
Generated by the C program below Download image ↗
FROM SOURCE TO SHAPE

Run it in three steps.

  1. Save the source.

    Download geo-fareyQuadraticRiver.c into a folder on your computer.

  2. Compile it.

    In that folder, run this with GCC or Clang:

    cc -std=c11 -O2 geo-fareyQuadraticRiver.c -lm -o geo-fareyQuadraticRiver
  3. Make the image.
    ./geo-fareyQuadraticRiver

    Open geo-fareyQuadraticRiver.svg in a browser to see the result.

On Windows with GCC, name the executable geo-fareyQuadraticRiver.exe and run it from the same folder.

THE RULE IN THE PROGRAM

How the picture is built

Generate Farey triangles using sums and differences of primitive rational vectors, map their vertices into the Poincare disk, and highlight dual-tree edges separating opposite signs of p^2-Dq^2.

Make it your own

DEPTH = 6 (1..8); RADICAND = 2 (2..30). A square radicand advances by one.

The highlighted river uses straight connections between chosen triangle centers. It is distinct from the separately drawn smooth geodesic with quadratic-irrational endpoints. The tessellation is truncated. Each run produces one image; use Graphic mode for the interactive animation.

The complete source

geo-fareyQuadraticRiver.c · 139 lines · no graphics libraries
Download .c ↓
/* 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;
}
All 40 explorations, ready to compile.Download all C examples ↓