research-apollonian.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, save, Board, Color, Ramp};3use mrlynum::factor::gcd;4use mrlynum::lattice::totients;56const TOP: i64 = 2048;7const DEEP: usize = 32;8const CIRCLES: usize = 2448;9const RESTING: usize = 323;10const THICK: f64 = 1.7;11const HAIR: f64 = 1.2;12const RULE: f64 = 2.0;13const OVER: f64 = 0.035;1415#[derive(Clone, Copy)]16struct Circle {17    k: i64,18    x: i64,19    y: i64,20}2122type Quad = [Circle; 4];2324fn root() -> Quad {25    [26        Circle { k: 0, x: 0, y: -1 },27        Circle { k: 0, x: 0, y: 1 },28        Circle { k: 2, x: 0, y: 1 },29        Circle { k: 2, x: 2, y: 1 },30    ]31}3233fn descartes(q: &Quad) -> bool {34    let sum: i64 = q.iter().map(|c| c.k).sum();35    let squares: i64 = q.iter().map(|c| c.k * c.k).sum();36    sum * sum == 2 * squares37}3839fn reflect(q: &Quad, i: usize) -> Circle {40    let (mut k, mut x, mut y) = (0i64, 0i64, 0i64);41    for (j, c) in q.iter().enumerate() {42        if j != i {43            k += c.k;44            x += c.x;45            y += c.y;46        }47    }48    Circle {49        k: 2 * k - q[i].k,50        x: 2 * x - q[i].x,51        y: 2 * y - q[i].y,52    }53}5455fn swap(q: &Quad, i: usize) -> Quad {56    let mut out = *q;57    out[i] = reflect(q, i);58    out59}6061fn grow(cap: i64) -> Vec<Circle> {62    let seed = root();63    let mut out = Vec::new();64    let mut stack: Vec<(Quad, usize)> = Vec::new();65    assert!(descartes(&seed));66    for i in 0..2 {67        if reflect(&seed, i).k <= cap {68            stack.push((swap(&seed, i), i));69        }70    }71    while let Some((q, last)) = stack.pop() {72        assert!(descartes(&q));73        out.push(q[last]);74        for j in 0..4 {75            if j == last {76                continue;77            }78            let next = reflect(&q, j);79            if next.k > cap || next.k <= q[j].k {80                continue;81            }82            stack.push((swap(&q, j), j));83        }84    }85    out86}8788fn fraction(c: Circle) -> (i64, i64) {89    let b = ((c.k / 2) as f64).sqrt().round() as i64;90    if b <= 0 {91        return (0, 0);92    }93    (c.x / (2 * b), b)94}9596fn mark(board: &mut Board, at: (f64, f64), r: f64, thick: f64, span: (f64, f64), color: Color) {97    let half = thick / 2.0;98    let reach = r + half + 1.0;99    let x0 = (at.0 - reach).max(span.0 - 1.0).max(0.0).floor() as usize;100    let x1 = (at.0 + reach).min(span.1 + 1.0).max(0.0).ceil() as usize;101    let y0 = (at.1 - reach).max(0.0).floor() as usize;102    let y1 = (at.1 + reach).max(0.0).ceil() as usize;103    for py in y0..y1.min(board.height) {104        for px in x0..x1.min(board.width) {105            let (fx, fy) = (px as f64 + 0.5, py as f64 + 0.5);106            let ring = (((fx - at.0).powi(2) + (fy - at.1).powi(2)).sqrt() - r).abs() - half;107            let cut = (span.0 - fx).max(fx - span.1);108            board.blend(px, py, color, 0.5 - ring.max(cut));109        }110    }111}112113fn main() -> Result<()> {114    let seed = root();115    let grown = grow(TOP);116    let phi = totients(DEEP);117    let nodes = phi[1..=DEEP].iter().sum::<u64>() as usize - 1;118    let resting: Vec<Circle> = grown.iter().copied().filter(|c| c.y == 1).collect();119120    assert_eq!(grown.len(), CIRCLES);121    assert_eq!(resting.len(), RESTING);122    assert_eq!(resting.len(), nodes);123    assert!(grown.iter().all(|c| c.k > 0 && c.x > 0 && c.x < c.k));124    let mut deepest = 1;125    for c in &resting {126        let (a, b) = fraction(*c);127        assert_eq!((c.k, c.x, c.y), (2 * b * b, 2 * a * b, 1));128        assert_eq!(gcd(a as usize, b as usize), 1);129        assert!(0 < a && a < b && b <= DEEP as i64);130        deepest = deepest.max(b);131    }132    assert_eq!(deepest, DEEP as i64);133    assert_eq!(fraction(seed[2]), (0, 1));134    assert_eq!(fraction(seed[3]), (1, 1));135136    let mut board = Board::square();137    let frame = board.frame(0.08);138    let span = (frame.x, frame.x + frame.w);139    let over = frame.w * OVER;140    let place = |c: &Circle| {141        (142            (143                frame.x + frame.w * c.x as f64 / c.k as f64,144                frame.y + frame.h * (1.0 - c.y as f64 / c.k as f64),145            ),146            frame.w / c.k as f64,147        )148    };149    for foot in [frame.y, frame.y + frame.h] {150        board.segment(151            (span.0 - over, foot),152            (span.1 + over, foot),153            RULE,154            ink::dim(),155        );156    }157    let rest = ink::fade(ink::dim(), 0.75);158    for c in grown.iter().filter(|c| c.y != 1) {159        let (at, r) = place(c);160        mark(&mut board, at, r, HAIR, span, rest);161    }162    let ramp = Ramp::tone(ink::blue(), ink::fg());163    let reach = (DEEP as f64).ln();164    for c in resting.iter().chain([seed[2], seed[3]].iter()) {165        let (at, r) = place(c);166        let (_, b) = fraction(*c);167        mark(168            &mut board,169            at,170            r,171            THICK,172            span,173            ramp.at((b as f64).ln() / reach),174        );175    }176    save("research-apollonian", &board)?;177    Ok(())178}