main.rs
13.2 kB · rust · 433 lines
1use mrlynum::spirograph::{2 cover, disc, pencils, point, representatives, signed_area, track, Kind, Pencil, Track,3};4use std::f64::consts::{PI, TAU};56const SIDES: [usize; 4] = [256, 512, 1024, 2048];7const REACH: f64 = 0.9;8const TRACE: usize = 200_001;910struct Case {11 name: String,12 path: Track,13 pens: Vec<Pencil>,14}1516struct Read {17 covered: f64,18 hole: f64,19 wall: f64,20 winding: f64,21 areas: f64,22}2324struct Limit {25 covered: f64,26 bar: f64,27 hole: f64,28 ratio: f64,29}3031fn seat(d: f64, quarter: usize) -> Pencil {32 let mut cell = (2i64, 0i64);33 for _ in 0..quarter % 4 {34 cell = (-cell.1, cell.0);35 }36 let angle = TAU * (quarter % 4) as f64 / 4.0;37 Pencil {38 x: d * angle.cos(),39 y: d * angle.sin(),40 seat: cell,41 kind: Kind::Fill,42 }43}4445fn copies(d: f64, count: usize) -> Vec<Pencil> {46 (0..count).map(|k| seat(d, k * (4 / count))).collect()47}4849fn carpet() -> Vec<u8> {50 vec![1, 1, 1, 1, 0, 1, 1, 1, 1]51}5253fn one(name: String, kind: &str, a: usize, b: usize, pens: Vec<Pencil>) -> Case {54 Case {55 name,56 path: track(kind, a, b, 4, 1).unwrap(),57 pens,58 }59}6061fn cases() -> Vec<Case> {62 let mut out = Vec::new();63 for kind in ["in", "out"] {64 for (a, b) in ratios() {65 out.push(one(66 format!("one pencil {kind} {a}/{b}"),67 kind,68 a,69 b,70 vec![seat(REACH, 0)],71 ));72 }73 }74 for (a, b) in [(7usize, 3usize), (5, 2)] {75 for count in [2usize, 4] {76 out.push(one(77 format!("{count} copies in {a}/{b}"),78 "in",79 a,80 b,81 copies(REACH, count),82 ));83 }84 }85 for (pens, label) in [("fill", "carpet fills"), ("corners", "carpet corners")] {86 for (kind, a, b) in [("in", 7usize, 3usize), ("in", 5, 2), ("out", 5, 8)] {87 out.push(one(88 format!("{label} {kind} {a}/{b}"),89 kind,90 a,91 b,92 pencils(&carpet(), 3, 3, pens, REACH, 0.0, 1).unwrap(),93 ));94 }95 }96 out97}9899fn ratios() -> Vec<(usize, usize)> {100 let mut out = Vec::new();101 for b in 1..=4usize {102 for a in b + 1..=9usize {103 if gcd(a, b) == 1 {104 out.push((a, b));105 }106 }107 }108 out109}110111fn gcd(a: usize, b: usize) -> usize {112 if b == 0 {113 a114 } else {115 gcd(b, a % b)116 }117}118119fn curves(case: &Case) -> Vec<usize> {120 representatives(&case.path, &case.pens, true)121}122123fn ladder(case: &Case) -> Vec<Read> {124 SIDES125 .iter()126 .map(|&side| {127 let out = cover(&case.path, &case.pens, true, 2, side).unwrap();128 Read {129 covered: out.covered,130 hole: out.hole,131 wall: out.wall,132 winding: out.winding,133 areas: out.areas,134 }135 })136 .collect()137}138139fn richer(coarse: f64, fine: f64) -> f64 {140 2.0 * fine - coarse141}142143fn limit(values: &[f64]) -> (f64, f64) {144 let last = richer(values[values.len() - 2], values[values.len() - 1]);145 let before = richer(values[values.len() - 3], values[values.len() - 2]);146 (last, (last - before).abs())147}148149fn limits(reads: &[Vec<Read>]) -> Vec<Limit> {150 reads151 .iter()152 .map(|rungs| {153 let covered: Vec<f64> = rungs.iter().map(|r| r.covered).collect();154 let holes: Vec<f64> = rungs.iter().map(|r| r.hole).collect();155 let steps: Vec<f64> = (1..covered.len())156 .map(|k| covered[k] - covered[k - 1])157 .collect();158 let ratio = steps159 .windows(2)160 .map(|pair| (pair[1] / pair[0]).abs())161 .fold(0.0_f64, f64::max);162 let (value, bar) = limit(&covered);163 Limit {164 covered: value,165 bar,166 hole: limit(&holes).0,167 ratio,168 }169 })170 .collect()171}172173fn far(case: &Case) -> f64 {174 case.path.wheel175 * case176 .pens177 .iter()178 .map(|p| p.x.hypot(p.y))179 .fold(0.0_f64, f64::max)180}181182fn rho_of(case: &Case) -> f64 {183 disc(&case.path, &case.pens).unwrap().radius - far(case)184}185186fn bound(case: &Case, side: usize) -> f64 {187 let bounds = disc(&case.path, &case.pens).unwrap();188 let rho = rho_of(case);189 let arcs: f64 = curves(case)190 .iter()191 .map(|&k| {192 TAU * case.path.ratio.1 as f64 * rho * (1.0 + case.pens[k].x.hypot(case.pens[k].y))193 })194 .sum();195 2.0 * arcs * (2.0 * bounds.radius / side as f64) / (PI * bounds.radius * bounds.radius)196}197198fn simple(case: &Case) -> f64 {199 let (rho, d) = (rho_of(case), far(case));200 rho * (rho + case.path.side * d * d / case.path.wheel) / (rho + d).powi(2)201}202203fn annulus(case: &Case) -> f64 {204 let (rho, d) = (rho_of(case), far(case));205 4.0 * rho * d / (rho + d).powi(2)206}207208fn inscribed(case: &Case) -> f64 {209 let bounds = disc(&case.path, &case.pens).unwrap();210 (bounds.hole / bounds.radius).powi(2)211}212213fn lone(case: &Case) -> bool {214 curves(case).len() == 1 && case.path.ratio.1 == 1215}216217fn the_disc(all: &[Case]) -> (f64, usize) {218 println!("THE DISC: the two radii against a trace of {TRACE} points per curve");219 println!(220 "{:<22} {:>7} {:>10} {:>12} {:>10} {:>12} {:>10}",221 "case", "curves", "radius", "traced max", "hole", "traced min", "gap"222 );223 let mut worst = 0.0_f64;224 for case in all {225 let bounds = disc(&case.path, &case.pens).unwrap();226 let (mut high, mut low) = (0.0_f64, f64::MAX);227 for &k in curves(case).iter() {228 for step in 0..TRACE {229 let s = case.path.total * step as f64 / (TRACE - 1) as f64;230 let (x, y) = point(&case.path, &case.pens[k], s);231 let radius = (x - bounds.x).hypot(y - bounds.y);232 high = high.max(radius);233 low = low.min(radius);234 }235 }236 let gap = (high - bounds.radius).abs().max((low - bounds.hole).abs());237 worst = worst.max(gap);238 println!(239 "{:<22} {:>7} {:>10.6} {:>12.6} {:>10.6} {:>12.6} {:>10.2e}",240 case.name,241 curves(case).len(),242 bounds.radius,243 high,244 bounds.hole,245 low,246 gap247 );248 }249 println!();250 (worst, all.len())251}252253fn the_winding(all: &[Case], reads: &[Vec<Read>]) -> (f64, usize) {254 println!("THE WINDING SELF-CHECK: the scanline mean against the Green closed form, every side");255 println!(256 "{:<22} {:>7} {:>12} {:>10} {:>10} {:>10} {:>10} {:>11}",257 "case", "curves", "closed form", "gap 256", "gap 512", "gap 1024", "gap 2048", "bound 2048"258 );259 let (mut worst, mut loose) = (0.0_f64, 0usize);260 for (case, rungs) in all.iter().zip(reads) {261 let gaps: Vec<f64> = rungs.iter().map(|r| (r.winding - r.areas).abs()).collect();262 for (rung, gap) in gaps.iter().enumerate() {263 worst = worst.max(*gap);264 if *gap > bound(case, SIDES[rung]) {265 loose += 1;266 }267 }268 println!(269 "{:<22} {:>7} {:>12.6} {:>10.2e} {:>10.2e} {:>10.2e} {:>10.2e} {:>11.2e}",270 case.name,271 curves(case).len(),272 rungs[3].areas,273 gaps[0],274 gaps[1],275 gaps[2],276 gaps[3],277 bound(case, SIDES[3])278 );279 }280 println!();281 (worst, loose)282}283284fn the_ladder(all: &[Case], reads: &[Vec<Read>], marks: &[Limit]) {285 println!("THE COVER: the shape's share of the disc, raster by raster, the wall counted in");286 println!(287 "{:<22} {:>7} {:>9} {:>9} {:>9} {:>9} {:>10} {:>9} {:>9} {:>9}",288 "case", "curves", "256", "512", "1024", "2048", "covered*", "bar", "wall 2048", "hole*"289 );290 for ((case, rungs), mark) in all.iter().zip(reads).zip(marks) {291 println!(292 "{:<22} {:>7} {:>9.6} {:>9.6} {:>9.6} {:>9.6} {:>10.6} {:>9.6} {:>9.6} {:>9.6}",293 case.name,294 curves(case).len(),295 rungs[0].covered,296 rungs[1].covered,297 rungs[2].covered,298 rungs[3].covered,299 mark.covered,300 mark.bar,301 rungs[3].wall,302 mark.hole303 );304 }305 println!();306}307308fn the_hole(all: &[Case], marks: &[Limit]) -> (usize, (f64, String), (f64, String), f64) {309 println!(310 "THE HOLE: the centre flood against its inscribed disc and, at b = 1, the closed form"311 );312 println!(313 "{:<22} {:>10} {:>11} {:>10} {:>11} {:>10}",314 "case", "hole*", "inscribed", "slack", "b = 1 form", "gap"315 );316 let mut leaks = 0;317 let mut low = (f64::MAX, String::new());318 let mut high = (0.0_f64, String::new());319 let mut worst = 0.0_f64;320 for (case, mark) in all.iter().zip(marks) {321 let (floor, form) = (inscribed(case), simple(case));322 let slack = mark.hole - floor;323 if slack < 0.0 {324 leaks += 1;325 }326 if slack < low.0 {327 low = (slack, case.name.clone());328 }329 if slack > high.0 {330 high = (slack, case.name.clone());331 }332 if lone(case) {333 worst = worst.max((mark.hole - form).abs());334 }335 println!(336 "{:<22} {:>10.6} {:>11.6} {:>10.6} {:>11} {:>10}",337 case.name,338 mark.hole,339 floor,340 slack,341 if lone(case) {342 format!("{form:.6}")343 } else {344 "-".to_string()345 },346 if lone(case) {347 format!("{:.2e}", (mark.hole - form).abs())348 } else {349 "-".to_string()350 }351 );352 }353 println!();354 (leaks, low, high, worst)355}356357fn the_candidates(all: &[Case], marks: &[Limit]) -> (usize, usize) {358 println!("THE CANDIDATES: the extrapolated cover against three closed forms");359 println!(360 "{:<22} {:>10} {:>9} {:>10} {:>10} {:>10} {:>10}",361 "case", "covered*", "bar", "simple", "annulus", "areas", "verdict"362 );363 let (mut held, mut tried) = (0usize, 0usize);364 for (case, mark) in all.iter().zip(marks) {365 let bar = mark.bar.max(1e-6);366 let forms = [simple(case), annulus(case), areas(case)];367 let names = ["simple", "annulus", "areas"];368 let kept: Vec<&str> = names369 .iter()370 .zip(&forms)371 .filter(|(_, form)| (**form - mark.covered).abs() < 10.0 * bar)372 .map(|(name, _)| *name)373 .collect();374 held += kept.len();375 tried += forms.len();376 println!(377 "{:<22} {:>10.6} {:>9.6} {:>10.6} {:>10.6} {:>10.6} {:>10}",378 case.name,379 mark.covered,380 mark.bar,381 forms[0],382 forms[1],383 forms[2],384 if kept.is_empty() {385 "none".to_string()386 } else {387 kept.join(",")388 }389 );390 }391 println!();392 (held, tried)393}394395fn areas(case: &Case) -> f64 {396 let bounds = disc(&case.path, &case.pens).unwrap();397 curves(case)398 .iter()399 .filter_map(|&k| signed_area(&case.path, &case.pens[k]))400 .sum::<f64>()401 / (PI * bounds.radius * bounds.radius)402}403404fn main() {405 let all = cases();406 let reads: Vec<Vec<Read>> = all.iter().map(ladder).collect();407 let marks = limits(&reads);408 println!("ROULETTE COVER: the shape between the walls of a circle roulette");409 println!("reach {REACH} wheel radii, seats exact, sides {SIDES:?}");410 println!("ratios: every a/b in lowest terms with 1 <= b <= 4 and b < a <= 9, inside and out");411 println!();412 let (radii, seen) = the_disc(&all);413 let (winding, loose) = the_winding(&all, &reads);414 the_ladder(&all, &reads, &marks);415 let (leaks, low, high, form) = the_hole(&all, &marks);416 let (held, tried) = the_candidates(&all, &marks);417 let flat = all418 .iter()419 .zip(&marks)420 .filter(|(case, _)| lone(case))421 .map(|(_, mark)| mark.covered.abs())422 .fold(0.0_f64, f64::max);423 let ones = all.iter().filter(|case| lone(case)).count();424 let ratio = marks.iter().map(|mark| mark.ratio).fold(0.0_f64, f64::max);425 println!("THE REPORT");426 println!("- the two disc radii meet the trace on all {seen} cases, worst gap {radii:.2e}.");427 println!("- the scanline winding meets the closed form on all {seen} cases at all four sides, worst gap {winding:.2e}, {loose} readings outside the perimeter bound.");428 println!("- successive ladder differences fall by a factor of at most {ratio:.3} on every case and both steps, the pixel law the extrapolation rests on.");429 println!("- b = 1 and one curve, {ones} cases: the curve is simple, so the cover is the wall alone, at most {flat:.6} extrapolated.");430 println!("- b = 1 and one curve: the centre flood is the whole inside, the signed closed form holding to {form:.2e}.");431 println!("- the centre flood holds its inscribed disc on every case, {leaks} leaks, slack from {:.6} at {} to {:.6} at {}.", low.0, low.1, high.0, high.1);432 println!("- {held} of {tried} candidate readings survive ten bars: no closed form shows itself for the cover.");433}