MAKE THE MATHEMATICS YOURSELF

Yellowstone geysers in C

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

Output generated by this C program: Yellowstone geysers
Generated by the C program below Download image ↗
FROM SOURCE TO SHAPE

Run it in three steps.

  1. Save the source.

    Download yellowstone.c into a folder on your computer.

  2. Compile it.

    In that folder, run this with GCC or Clang:

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

    Open yellowstone.svg in a browser to see the result.

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

THE RULE IN THE PROGRAM

How the picture is built

Start 1,2,3. Each later term is the least unused positive integer sharing a factor greater than one with the term two places earlier, while remaining coprime to the previous term. Join index-value points.

Make it your own

N=300 terms; SEARCH_LIMIT=100000 candidate values. Increase the search bound if the program reports that the prefix cannot be completed.

The spikes come from the exact greedy gcd rule. Apparent bands in a finite graph are not a proof of an asymptotic formula. N must be at least three, SEARCH_LIMIT at least three; the bounded search can stop explicitly. Each run produces one image; use Graphic mode for the interactive animation.

The complete source

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

/* Greedy Yellowstone permutation: begin 1,2,3; share a factor with n-2,
   be coprime to n-1, and never reuse a value. Edit N (search is bounded). */
static void draw(void) {
    enum { N=300, SEARCH_LIMIT=100000 };
    int a[N];unsigned char *used=calloc(SEARCH_LIMIT+1,1);
    if(!used){nt_text(50,80,"Not enough memory.");return;}
    a[0]=1;a[1]=2;a[2]=3;used[1]=used[2]=used[3]=1;
    int maximum=3;
    for(int i=3;i<N;++i) {
        int chosen=0;
        for(int k=2;k<=SEARCH_LIMIT;++k)
            if(!used[k]&&nt_gcd(k,a[i-2])>1&&nt_gcd(k,a[i-1])==1){chosen=k;break;}
        if(!chosen){free(used);nt_text(50,80,"Increase SEARCH_LIMIT for this prefix.");return;}
        a[i]=chosen;used[chosen]=1;if(chosen>maximum)maximum=chosen;
    }
    nt_line(65,630,940,630,.2,.4,1);nt_line(65,60,65,630,.2,.4,1);
    for(int i=0;i<N;++i) {
        double x=65+875*i/(double)(N-1),y=630-560.0*a[i]/maximum;
        if(i)nt_line(65+875*(i-1)/(double)(N-1),630-560.0*a[i-1]/maximum,x,y,i/(double)N,.72,.85);
        nt_dot(x,y,1.4,i/(double)N,.9);
    }
    char label[64];snprintf(label,sizeof(label),"Term index (1 to %d)",N);nt_text(65,665,label);nt_text(70,45,"Value: arithmetic constraints create the spikes");
    free(used);
}

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] : "yellowstone.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 ↓