wiki-horocycle-flow.rs

1use std::collections::BTreeSet;23use figures::{ink, save, Board, Color, Ramp};4use mrlyrs::core::error::Result;5use mrlyrs::num::factor::gcd;67const HEIGHT: f64 = 0.01;8const ROOF: f64 = 1.6;9const SAMPLES: usize = 20000;10const THICK: f64 = 1.5;11const EDGE: f64 = 2.6;12const STEP: f64 = 1.5;1314struct Canvas {15    left: f64,16    bottom: f64,17    scale: f64,18    floor: f64,19    top: f64,20}2122impl Canvas {23    fn at(&self, x: f64, y: f64) -> (f64, f64) {24        (25            self.left + x * self.scale,26            self.bottom - (y - self.floor) * self.scale,27        )28    }2930    fn inside(&self, x: f64, y: f64) -> bool {31        x.abs() <= 0.5 && x * x + y * y >= 1.0 && y <= self.top32    }33}3435fn main() -> Result<()> {36    let mut board = Board::square();37    let frame = board.frame(0.08);38    let floor = 3f64.sqrt() / 2.0;39    let scale = frame.h / ROOF;40    let canvas = Canvas {41        left: frame.x + frame.w / 2.0,42        bottom: frame.y + frame.h,43        scale,44        floor,45        top: floor + ROOF,46    };47    let cusps = cusps();48    let deepest = cusps.iter().map(|&(_, c)| c).max().unwrap_or(1);49    let ramp = Ramp::tone(ink::blue(), ink::pink());50    let mut drawn = 0usize;51    for &(a, c) in &cusps {52        let tone = ramp.at((c - 1) as f64 / (deepest - 1).max(1) as f64);53        drawn += trace(&mut board, &canvas, a, c, tone);54    }55    outline(&mut board, &canvas);56    let hits = fold(&cusps);57    assert_eq!((cusps.len(), drawn, hits), (243, 110, 239));58    save("wiki-horocycle-flow", &board)?;59    Ok(())60}6162fn cusps() -> Vec<(i64, i64)> {63    let floor = 3f64.sqrt() / 2.0;64    let mut out = Vec::new();65    let mut c = 1i64;66    while 1.0 / ((c * c) as f64 * HEIGHT) >= floor {67        let radius = 1.0 / (2.0 * (c * c) as f64 * HEIGHT);68        let reach = ((radius + 0.5) * c as f64).floor() as i64;69        for a in -reach..=reach {70            if gcd(a.unsigned_abs() as u128, c as u128) == 1 {71                out.push((a, c));72            }73        }74        c += 1;75    }76    out77}7879fn trace(board: &mut Board, canvas: &Canvas, a: i64, c: i64, tone: Color) -> usize {80    let centre = a as f64 / c as f64;81    let radius = 1.0 / (2.0 * (c * c) as f64 * HEIGHT);82    let dt = STEP / (radius * canvas.scale);83    let steps = (std::f64::consts::TAU / dt).ceil() as usize;84    let mut run: Vec<(f64, f64)> = Vec::new();85    let mut runs = 0usize;86    for k in 0..=steps {87        let t = k as f64 * std::f64::consts::TAU / steps as f64 - std::f64::consts::FRAC_PI_2;88        let x = centre + radius * t.cos();89        let y = radius + radius * t.sin();90        if canvas.inside(x, y) {91            run.push(canvas.at(x, y));92        } else if !run.is_empty() {93            board.polyline(&run, THICK, tone);94            run.clear();95            runs += 1;96        }97    }98    if !run.is_empty() {99        board.polyline(&run, THICK, tone);100        runs += 1;101    }102    runs103}104105fn outline(board: &mut Board, canvas: &Canvas) {106    let tone = ink::fg();107    let floor = canvas.floor;108    let top = canvas.top;109    board.segment(canvas.at(-0.5, floor), canvas.at(-0.5, top), EDGE, tone);110    board.segment(canvas.at(0.5, floor), canvas.at(0.5, top), EDGE, tone);111    let pts: Vec<(f64, f64)> = (0..=240)112        .map(|k| {113            let t = std::f64::consts::PI * (1.0 / 3.0 + k as f64 / 720.0);114            canvas.at(t.cos(), t.sin())115        })116        .collect();117    board.polyline(&pts, EDGE, tone);118}119120fn fold(cusps: &[(i64, i64)]) -> usize {121    let known: BTreeSet<(i64, i64)> = cusps.iter().copied().collect();122    let mut hit = BTreeSet::new();123    for k in 0..SAMPLES {124        let (mut x, mut y) = ((k as f64 + 0.5) / SAMPLES as f64, HEIGHT);125        let mut m = [[1i64, 0], [0, 1]];126        loop {127            let n = x.round();128            x -= n;129            m = [130                [m[0][0] - n as i64 * m[1][0], m[0][1] - n as i64 * m[1][1]],131                m[1],132            ];133            let norm = x * x + y * y;134            if norm >= 1.0 {135                break;136            }137            (x, y) = (-x / norm, y / norm);138            m = [[-m[1][0], -m[1][1]], m[0]];139        }140        let (mut a, mut c) = (m[0][0], m[1][0]);141        if c < 0 {142            (a, c) = (-a, -c);143        }144        let centre = a as f64 / c as f64;145        let radius = 1.0 / (2.0 * (c * c) as f64 * HEIGHT);146        let off = ((x - centre).powi(2) + (y - radius).powi(2)).sqrt() - radius;147        assert!(off.abs() < 1e-6 * radius);148        assert!(known.contains(&(a, c)));149        hit.insert((a, c));150    }151    hit.len()152}