paper-unconditional-mertens-at-large-base.rs
4.0 kB · rust · 156 lines
1use figures::{ink, save, Board, Color};2use mrlyrs::core::error::Result;3use std::f64::consts::TAU;45const NAME: &str = "paper-unconditional-mertens-at-large-base";6const BASE: i64 = 10;7const LEVEL: u32 = 3;8const Z: i64 = 8;9const COUNTS: [usize; 4] = [732, 202, 26, 40];1011// DISSECTION1213fn fraction(a: i64, y: i64, cap: i64) -> (i64, i64) {14 let (mut p0, mut q0, mut p1, mut q1) = (0i64, 1i64, 1i64, 0i64);15 let (mut n, mut d) = (a, y);16 let mut best = (0, 1);17 while d > 0 {18 let c = n / d;19 let (p2, q2) = (c * p1 + p0, c * q1 + q0);20 if q2 > cap {21 break;22 }23 (p0, q0, p1, q1) = (p1, q1, p2, q2);24 best = (p2, q2);25 (n, d) = (d, n - c * d);26 }27 best28}2930fn gcd(a: i64, b: i64) -> i64 {31 if b == 0 {32 a.abs()33 } else {34 gcd(b, a % b)35 }36}3738fn smooth(mut d: i64) -> bool {39 loop {40 let g = gcd(d, BASE);41 if g == 1 {42 return d == 1;43 }44 while d % g == 0 {45 d /= g;46 }47 }48}4950fn region(a: i64, y: i64) -> usize {51 let cap = (y as f64).powf(0.6).floor() as i64;52 let low = (y as f64).powf(0.4);53 let (l, d) = fraction(a, y, cap);54 let h = (a * d - l * y).abs();55 assert!(d <= cap && h * cap <= y && gcd(l, d) == 1);56 if d as f64 >= low {57 058 } else if d < Z && h < Z {59 if smooth(d) {60 361 } else {62 263 }64 } else {65 166 }67}6869// DRAWING7071fn band(72 board: &mut Board,73 center: (f64, f64),74 radii: (f64, f64),75 angles: (f64, f64),76 color: Color,77) {78 let at = |r: f64, t: f64| (center.0 + r * t.sin(), center.1 - r * t.cos());79 let steps = (((angles.1 - angles.0) / 0.004).ceil() as usize).max(1);80 let mut pts = Vec::with_capacity(2 * steps + 2);81 for i in 0..=steps {82 pts.push(at(83 radii.1,84 angles.0 + (angles.1 - angles.0) * i as f64 / steps as f64,85 ));86 }87 for i in (0..=steps).rev() {88 pts.push(at(89 radii.0,90 angles.0 + (angles.1 - angles.0) * i as f64 / steps as f64,91 ));92 }93 board.polygon(&pts, color);94}9596fn runs(regions: &[usize], pass: usize) -> Vec<(usize, usize)> {97 let n = regions.len();98 let start = (0..n).find(|a| regions[*a] != pass).unwrap_or(0);99 let mut out = Vec::new();100 let mut open: Option<usize> = None;101 for step in 1..=n {102 let a = start + step;103 let inside = regions[a % n] == pass;104 match (open, inside) {105 (None, true) => open = Some(a),106 (Some(first), false) => {107 out.push((first, a));108 open = None;109 }110 _ => {}111 }112 }113 if let Some(first) = open {114 out.push((first, start + n + 1));115 }116 out117}118119fn main() -> Result<()> {120 let y = BASE.pow(LEVEL);121 let regions: Vec<usize> = (0..y).map(|a| region(a, y)).collect();122 let mut counts = [0usize; 4];123 for r in ®ions {124 counts[*r] += 1;125 }126 assert_eq!(counts, COUNTS);127 assert_eq!(counts.iter().sum::<usize>(), y as usize);128129 let mut board = Board::square();130 let frame = board.frame(0.08);131 let center = frame.center();132 let outer = frame.radius();133 let slot = TAU / y as f64;134 let width = 0.12 * outer;135 let tracks = [0.28, 0.48, 0.68, 0.88];136 let inks = [ink::blue(), ink::dim(), ink::orange(), ink::yellow()];137 for pass in 0..4 {138 let mid = tracks[pass] * outer + width / 2.0;139 board.ring(center.0, center.1, mid, 1.0, ink::line());140 let mut held = 0;141 for (first, end) in runs(®ions, pass) {142 held += end - first;143 let angles = ((first as f64 - 0.5) * slot, (end as f64 - 0.5) * slot);144 band(145 &mut board,146 center,147 (mid - width / 2.0, mid + width / 2.0),148 angles,149 inks[pass],150 );151 }152 assert_eq!(held, counts[pass]);153 }154 save(NAME, &board)?;155 Ok(())156}