main.rs
33.3 kB · rust · 1014 lines
1use mrlynum::spirograph::{distinct, pencils, point, trace, track, Kind, Pencil};2use std::f64::consts::{PI, TAU};34// THE CURVE56#[derive(Clone, Copy)]7struct Wheel {8 a: usize,9 b: usize,10 side: f64,11 amp: f64,12}1314fn wheel(a: usize, b: usize, inside: bool) -> Wheel {15 let side = if inside { -1.0 } else { 1.0 };16 Wheel {17 a,18 b,19 side,20 amp: (a as f64 / b as f64 + side).abs(),21 }22}2324fn lean(w: &Wheel, rho: f64, x: f64) -> f64 {25 let k = (rho - w.amp) / (rho + w.amp);26 if x >= PI - 1e-13 {27 return PI / 2.0 + if k > 0.0 { PI / 2.0 } else { -PI / 2.0 };28 }29 x / 2.0 + (k * (x / 2.0).tan()).atan()30}3132fn sweep(w: &Wheel, rho: f64, x: f64) -> f64 {33 w.b as f64 * w.side * x / w.a as f64 + lean(w, rho, x)34}3536fn ring(w: &Wheel, rho: f64, x: f64) -> f64 {37 (w.amp * w.amp + rho * rho + 2.0 * w.amp * rho * x.cos()).sqrt()38}3940fn mark_at(w: &Wheel, rho: f64, r: f64) -> f64 {41 let c = ((r * r - w.amp * w.amp - rho * rho) / (2.0 * w.amp * rho)).clamp(-1.0, 1.0);42 4.0 * w.a as f64 * sweep(w, rho, c.acos()) / PI43}4445// THE SEATS4647#[derive(Clone, Copy)]48struct Seat {49 half: (f64, f64),50 corner: bool,51}5253fn carpet() -> Vec<Seat> {54 let mut out = Vec::new();55 for i in 0..3i64 {56 for j in 0..3i64 {57 if i == 1 && j == 1 {58 continue;59 }60 let (hx, hy) = ((j - 1) as f64, (1 - i) as f64);61 out.push(Seat {62 half: (hx, hy),63 corner: hx != 0.0 && hy != 0.0,64 });65 }66 }67 out68}6970fn reach_unit() -> f64 {71 1.0 / (1.5_f64).hypot(1.5)72}7374fn seat_arm(seat: &Seat, t: f64) -> f64 {75 let u = t * reach_unit();76 (seat.half.0 * u).hypot(seat.half.1 * u)77}7879fn seat_oct(w: &Wheel, seat: &Seat) -> i64 {80 let dial = match (seat.half.0 as i64, seat.half.1 as i64) {81 (1, 0) => 0,82 (1, 1) => 1,83 (0, 1) => 2,84 (-1, 1) => 3,85 (-1, 0) => 4,86 (-1, -1) => 5,87 (0, -1) => 6,88 _ => 7,89 };90 (-(w.b as i64) * (w.side as i64) * dial).rem_euclid(8)91}9293// THE CURVES9495fn band(w: &Wheel, rho: f64) -> (f64, f64) {96 ((w.amp - rho).abs(), w.amp + rho)97}9899#[derive(Clone, Copy)]100struct Curve {101 arm: f64,102 oct: i64,103 corner: bool,104}105106fn curves(w: &Wheel, seats: &[Seat], t: f64) -> (Vec<Curve>, Vec<usize>) {107 let mut out: Vec<Curve> = Vec::new();108 let mut owner = Vec::new();109 for seat in seats {110 let arm = seat_arm(seat, t);111 let oct = seat_oct(w, seat);112 let at = out113 .iter()114 .position(|c: &Curve| (c.arm - arm).abs() < 1e-12 && c.oct == oct);115 match at {116 Some(k) => owner.push(k),117 None => {118 owner.push(out.len());119 out.push(Curve {120 arm,121 oct,122 corner: seat.corner,123 });124 }125 }126 }127 (out, owner)128}129130// THE SOLVER131132fn meeting(cs: &[Curve], m: i64, n: i64, flat: bool) -> (usize, usize, i64, Vec<usize>) {133 let mut best = (0usize, 0usize, 0i64, Vec::new());134 for slot in 0..8i64 {135 let mut ends = 0usize;136 let mut seen: Vec<usize> = Vec::new();137 for (k, c) in cs.iter().enumerate() {138 let mark = if c.corner { m } else { n };139 let signs: &[i64] = if flat && !c.corner { &[1] } else { &[1, -1] };140 for &sign in signs {141 if (c.oct + sign * mark).rem_euclid(8) == slot {142 ends += 1;143 if !seen.contains(&k) {144 seen.push(k);145 }146 }147 }148 }149 if ends > best.0 {150 seen.sort_unstable();151 best = (ends, seen.len(), slot, seen);152 }153 }154 best155}156157#[derive(Clone, Copy)]158struct Align {159 reach: f64,160 turn: f64,161 ring: f64,162 mark_c: i64,163 mark_e: i64,164 ends: usize,165 kinds: usize,166 slot: i64,167 flat: bool,168}169170fn turn_for(w: &Wheel, arm: f64, mark: f64) -> f64 {171 let (mut lo, mut hi) = (0.0_f64, PI);172 for _ in 0..64 {173 let mid = 0.5 * (lo + hi);174 let got = 4.0 * w.a as f64 * sweep(w, arm, mid) / PI;175 if (got - mark) * w.side < 0.0 {176 lo = mid;177 } else {178 hi = mid;179 }180 }181 0.5 * (lo + hi)182}183184fn tame(w: &Wheel) -> f64 {185 1.5 * 1.0_f64.min(w.amp)186}187188fn arms(t: f64) -> (f64, f64) {189 (2.0 * t / 3.0, 2.0_f64.sqrt() * t / 3.0)190}191192fn edge_mark(w: &Wheel, t: f64, x: f64) -> Option<f64> {193 let (arm_c, arm_e) = arms(t);194 let r = ring(w, arm_c, x);195 let (lo, hi) = band(w, arm_e);196 if r <= lo || r >= hi {197 return None;198 }199 Some(mark_at(w, arm_e, r))200}201202fn alignments(w: &Wheel, tlo: f64, thi: f64, grid: usize) -> Vec<Align> {203 let sample = curves(w, &carpet(), 0.5 * (tlo + thi)).0;204 let reach = |i: usize| tlo + (thi - tlo) * i as f64 / grid as f64;205 let full = 4.0 * w.b as f64 * w.side;206 let marks: Vec<i64> = if w.side < 0.0 {207 ((full.ceil() as i64 + 1)..0).collect()208 } else {209 (1..(full.floor() as i64)).collect()210 };211 let hold = |t: f64, m: i64| -> f64 {212 let (arm_c, arm_e) = arms(t);213 mark_at(w, arm_e, ring(w, arm_c, turn_for(w, arm_c, m as f64)))214 };215 let mut out: Vec<Align> = Vec::new();216 for &m in &marks {217 let mut back = hold(reach(0), m);218 for cell in 1..=grid {219 let front = hold(reach(cell), m);220 {221 let (d0, d1) = (back, front);222 let (dlo, dhi) = (d0.min(d1), d0.max(d1));223 for n in (dlo.floor() as i64 + 1)..=(dhi.ceil() as i64) {224 if (n as f64) <= dlo || (n as f64) >= dhi {225 continue;226 }227 let (ends, kinds, slot, _) = meeting(&sample, m, n, false);228 if ends < 3 {229 continue;230 }231 let (mut u, mut v) = (reach(cell - 1), reach(cell));232 for _ in 0..90 {233 let mid = 0.5 * (u + v);234 if (hold(mid, m) - n as f64) * (d0 - n as f64) > 0.0 {235 u = mid;236 } else {237 v = mid;238 }239 }240 let t = 0.5 * (u + v);241 let (arm_c, _) = arms(t);242 let x = turn_for(w, arm_c, m as f64);243 if edge_mark(w, t, x).is_none() {244 continue;245 }246 out.push(Align {247 reach: t,248 turn: x,249 ring: ring(w, arm_c, x),250 mark_c: m,251 mark_e: n,252 ends,253 kinds,254 slot,255 flat: false,256 });257 }258 }259 back = front;260 }261 }262 for &m in &marks {263 for (n, outer) in [(0i64, true), (full.round() as i64, false)] {264 let (ends, kinds, slot, _) = meeting(&sample, m, n, true);265 if ends < 3 {266 continue;267 }268 let gap = |t: f64| -> f64 {269 let (arm_c, arm_e) = arms(t);270 let r = ring(w, arm_c, turn_for(w, arm_c, m as f64));271 let (lo, hi) = band(w, arm_e);272 r - if outer { hi } else { lo }273 };274 let mut back = gap(reach(0));275 for cell in 1..=grid {276 let front = gap(reach(cell));277 if back * front < 0.0 {278 let (mut u, mut v) = (reach(cell - 1), reach(cell));279 for _ in 0..90 {280 let mid = 0.5 * (u + v);281 if gap(mid) * back > 0.0 {282 u = mid;283 } else {284 v = mid;285 }286 }287 let t = 0.5 * (u + v);288 let (arm_c, _) = arms(t);289 let x = turn_for(w, arm_c, m as f64);290 out.push(Align {291 reach: t,292 turn: x,293 ring: ring(w, arm_c, x),294 mark_c: m,295 mark_e: n,296 ends,297 kinds,298 slot,299 flat: true,300 });301 }302 back = front;303 }304 }305 }306 out.sort_by(|p, q| p.reach.partial_cmp(&q.reach).unwrap());307 out308}309310// THE CENSUS311312struct Node {313 ring: f64,314 angle: f64,315 curves: Vec<usize>,316 branches: usize,317}318319fn arc(w: &Wheel, cs: &[Curve], k: usize, sign: f64, r: f64) -> Option<f64> {320 let (lo, hi) = band(w, cs[k].arm);321 if r < lo || r > hi {322 return None;323 }324 Some(cs[k].oct as f64 / 8.0 + sign * mark_at(w, cs[k].arm, r) / 8.0)325}326327fn census(w: &Wheel, cs: &[Curve], grid: usize) -> Vec<Node> {328 let ports: Vec<(usize, f64)> = (0..cs.len()).flat_map(|k| [(k, 1.0), (k, -1.0)]).collect();329 let mut raw: Vec<(f64, f64, usize, usize)> = Vec::new();330 for i in 0..ports.len() {331 for j in (i + 1)..ports.len() {332 if ports[i].0 == ports[j].0 && ports[i].1 == ports[j].1 {333 continue;334 }335 let (bi, bj) = (band(w, cs[ports[i].0].arm), band(w, cs[ports[j].0].arm));336 let (lo, hi) = (bi.0.max(bj.0), bi.1.min(bj.1));337 if hi <= lo {338 continue;339 }340 let (mid, half) = (0.5 * (lo + hi), 0.5 * (hi - lo));341 let ray = |s: f64| mid - half * (PI * s).cos();342 let gap = |s: f64| {343 let r = ray(s).clamp(lo, hi);344 let p = arc(w, cs, ports[i].0, ports[i].1, r);345 let q = arc(w, cs, ports[j].0, ports[j].1, r);346 match (p, q) {347 (Some(p), Some(q)) => Some(p - q),348 _ => None,349 }350 };351 let mut back = match gap(0.0) {352 Some(d) => d,353 None => continue,354 };355 for cell in 1..=grid {356 let s1 = cell as f64 / grid as f64;357 let s0 = (cell - 1) as f64 / grid as f64;358 let front = match gap(s1) {359 Some(d) => d,360 None => continue,361 };362 let (dlo, dhi) = (back.min(front), back.max(front));363 for goal in (dlo.floor() as i64 + 1)..=(dhi.ceil() as i64) {364 let goal = goal as f64;365 if goal <= dlo || goal >= dhi {366 continue;367 }368 let (mut u, mut v) = (s0, s1);369 for _ in 0..90 {370 let mid = 0.5 * (u + v);371 match gap(mid) {372 Some(d) => {373 if (d - goal) * (back - goal) > 0.0 {374 u = mid;375 } else {376 v = mid;377 }378 }379 None => break,380 }381 }382 let seat = 0.5 * (u + v);383 if seat < 1e-11 || seat > 1.0 - 1e-11 {384 continue;385 }386 let r = ray(seat).clamp(lo, hi);387 if let Some(p) = arc(w, cs, ports[i].0, ports[i].1, r) {388 raw.push((r, p - p.floor(), i, j));389 }390 }391 back = front;392 }393 }394 }395 let reach = cs.iter().map(|c| band(w, c.arm).1).fold(0.0_f64, f64::max);396 raw.sort_by(|p, q| p.0.partial_cmp(&q.0).unwrap());397 let mut out: Vec<Node> = Vec::new();398 let mut pool: Vec<Vec<(f64, f64, usize, usize)>> = Vec::new();399 for row in raw {400 let hit = pool.iter().position(|g: &Vec<(f64, f64, usize, usize)>| {401 (g[0].0 - row.0).abs() < 1e-9 * reach && {402 let d = (g[0].1 - row.1).abs();403 d.min(1.0 - d) < 1e-9404 }405 });406 match hit {407 Some(k) => pool[k].push(row),408 None => pool.push(vec![row]),409 }410 }411 for group in pool {412 let mut seen: Vec<usize> = Vec::new();413 let mut ends: Vec<usize> = Vec::new();414 for row in &group {415 for port in [row.2, row.3] {416 if !ends.contains(&port) {417 ends.push(port);418 }419 if !seen.contains(&ports[port].0) {420 seen.push(ports[port].0);421 }422 }423 }424 seen.sort_unstable();425 out.push(Node {426 ring: group[0].0,427 angle: group[0].1,428 branches: ends.len(),429 curves: seen,430 });431 }432 out433}434435// THE TRACES436437const CARPET: [u8; 9] = [1, 1, 1, 1, 0, 1, 1, 1, 1];438439fn near(px: f64, py: f64, x0: f64, y0: f64, x1: f64, y1: f64) -> f64 {440 let (dx, dy) = (x1 - x0, y1 - y0);441 let len = dx * dx + dy * dy;442 let s = if len > 0.0 {443 (((px - x0) * dx + (py - y0) * dy) / len).clamp(0.0, 1.0)444 } else {445 0.0446 };447 (px - x0 - s * dx).hypot(py - y0 - s * dy)448}449450fn probe(w: &Wheel, t: f64, r: f64, angle: f64, samples: usize) -> Vec<f64> {451 let kind = if w.side < 0.0 { "in" } else { "out" };452 let path = track(kind, w.a, w.b, 4, 1).unwrap();453 let pens = pencils(&CARPET, 3, 3, "fill", t, 0.0, 1).unwrap();454 let pts = trace(&path, &pens, samples).unwrap();455 let hit = w.b as f64 * r;456 let th = angle * TAU / w.a as f64;457 let (px, py) = (hit * th.cos(), hit * th.sin());458 (0..pens.len())459 .map(|k| {460 let head = k * samples * 2;461 let mut best = f64::MAX;462 for i in 0..samples - 1 {463 best = best.min(near(464 px,465 py,466 pts[head + 2 * i] as f64,467 pts[head + 2 * i + 1] as f64,468 pts[head + 2 * i + 2] as f64,469 pts[head + 2 * i + 3] as f64,470 ));471 }472 best473 })474 .collect()475}476477fn probe_wide(w: &Wheel, t: f64, r: f64, angle: f64, samples: usize) -> Vec<f64> {478 let kind = if w.side < 0.0 { "in" } else { "out" };479 let path = track(kind, w.a, w.b, 4, 1).unwrap();480 let pens = pencils(&CARPET, 3, 3, "fill", t, 0.0, 1).unwrap();481 let hit = w.b as f64 * r;482 let th = angle * TAU / w.a as f64;483 let (px, py) = (hit * th.cos(), hit * th.sin());484 pens.iter()485 .map(|pen| {486 let mut best = f64::MAX;487 let mut back = point(&path, pen, 0.0);488 for i in 1..samples {489 let front = point(&path, pen, path.total * i as f64 / (samples - 1) as f64);490 best = best.min(near(px, py, back.0, back.1, front.0, front.1));491 back = front;492 }493 best494 })495 .collect()496}497498fn f32_floor(v: f64) -> f64 {499 let a = v as f32;500 0.5 * (f32::from_bits(a.to_bits() + 1) - a) as f64501}502503fn residue(w: &Wheel, arm: f64, x: f64, mark: i64) -> f64 {504 let om = (x.cos(), x.sin());505 let mut left = (1.0, 0.0);506 let mut right = (1.0, 0.0);507 let pow = w.a as i64 + 2 * w.b as i64 * w.side as i64;508 let base = (w.amp + arm * om.0, arm * om.1);509 let flip = (w.amp * om.0 + arm, w.amp * om.1);510 for _ in 0..w.a {511 left = (512 left.0 * base.0 - left.1 * base.1,513 left.0 * base.1 + left.1 * base.0,514 );515 right = (516 right.0 * flip.0 - right.1 * flip.1,517 right.0 * flip.1 + right.1 * flip.0,518 );519 }520 let spin = pow as f64 * x;521 left = (522 left.0 * spin.cos() - left.1 * spin.sin(),523 left.0 * spin.sin() + left.1 * spin.cos(),524 );525 let quarter = mark.rem_euclid(4);526 let unit = [(1.0, 0.0), (0.0, 1.0), (-1.0, 0.0), (0.0, -1.0)][quarter as usize];527 let want = (528 right.0 * unit.0 - right.1 * unit.1,529 right.0 * unit.1 + right.1 * unit.0,530 );531 ((left.0 - want.0).hypot(left.1 - want.1)) / (left.0.hypot(left.1)).max(1e-12)532}533534// THE RUN535536fn defect(w: &Wheel, t: f64, m: f64, n: f64) -> f64 {537 let (arm_c, arm_e) = (2.0 * t / 3.0, 2.0_f64.sqrt() * t / 3.0);538 let r = ring(w, arm_c, turn_for(w, arm_c, m));539 mark_at(w, arm_e, r) - n540}541542fn deck() -> Vec<(usize, usize, bool)> {543 let mut out = Vec::new();544 for b in 1..=8usize {545 for a in (b + 1)..=13usize {546 let (mut p, mut q) = (a, b);547 while q > 0 {548 let r = p % q;549 p = q;550 q = r;551 }552 if p != 1 {553 continue;554 }555 for inside in [true, false] {556 out.push((a, b, inside));557 }558 }559 }560 out561}562563fn reduction() {564 println!("THE REDUCTION");565 let seats = carpet();566 let (mut worst, mut counts) = (0.0_f64, 0usize);567 for &(a, b, inside) in &deck() {568 let w = wheel(a, b, inside);569 let kind = if inside { "in" } else { "out" };570 let path = track(kind, a, b, 4, 1).unwrap();571 let pens = pencils(&CARPET, 3, 3, "fill", 0.83, 0.0, 1).unwrap();572 for pen in &pens {573 for k in 0..29 {574 let psi = TAU * k as f64 / 29.0;575 let (x, y) = point(&path, pen, path.total * psi / TAU);576 let spin = w.b as f64 * psi;577 let seat = (w.b as f64 + w.side * w.a as f64) * psi;578 let zx = b as f64 * (w.amp * spin.cos() + pen.x * seat.cos() - pen.y * seat.sin());579 let zy = b as f64 * (w.amp * spin.sin() + pen.x * seat.sin() + pen.y * seat.cos());580 worst = worst.max((x - zx).hypot(y - zy));581 }582 }583 let (cs, _) = curves(&w, &seats, 0.83);584 if cs.len() == distinct(&path, &pens, true) {585 counts += 1;586 }587 }588 println!(" z(psi) = r e^(i b psi) (A + p e^(i eps a psi)), A = abs(a/b + eps), at reach 0.83 and 29 phases a seat, worst gap to spirograph::point {worst:.3e}");589 println!(590 " offset classes against spirograph::distinct: {counts} of {} tracks agree",591 deck().len()592 );593 let w = wheel(7, 3, true);594 let path = track("in", 7, 3, 4, 1).unwrap();595 let (alpha, arm) = (0.3_f64, 0.5_f64);596 let seat = |turn: f64| Pencil {597 x: arm * turn.cos(),598 y: arm * turn.sin(),599 seat: (0, 0),600 kind: Kind::Fill,601 };602 let (plain, turned) = (seat(0.0), seat(alpha));603 let spin = w.b as f64 * alpha / w.a as f64;604 let mut gaps = [0.0_f64; 2];605 for (k, lean) in [-1.0_f64, 1.0].iter().enumerate() {606 let shift = lean * alpha / (w.side * w.a as f64);607 let angle = lean * w.side * spin;608 for j in 0..3 {609 let psi = TAU * (0.17 + 0.31 * j as f64);610 let (x0, y0) = point(&path, &plain, path.total * psi / TAU);611 let (x1, y1) = point(612 &path,613 &turned,614 path.total * (psi + shift).rem_euclid(TAU) / TAU,615 );616 let (c, sn) = (angle.cos(), angle.sin());617 gaps[k] = gaps[k].max((x1 - (x0 * c - y0 * sn)).hypot(y1 - (x0 * sn + y0 * c)));618 }619 }620 println!(621 " turning the seat by {alpha} on 7/3 inside turns the curve by minus b eps alpha over a: worst gap {:.3e}, against {:.3e} for the plus sign",622 gaps[0], gaps[1]623 );624}625626fn tally(nodes: &[Node], count: usize) -> (usize, usize, usize, usize, usize) {627 let mut own = vec![0usize; count];628 let mut duo = vec![0usize; count * count];629 let mut branch = 0;630 for node in nodes {631 branch = branch.max(node.branches);632 if node.curves.len() == 1 {633 own[node.curves[0]] += 1;634 } else if node.curves.len() == 2 {635 duo[node.curves[0] * count + node.curves[1]] += 1;636 }637 }638 let pairs: Vec<usize> = (0..count)639 .flat_map(|i| ((i + 1)..count).map(move |j| (i, j)))640 .map(|(i, j)| duo[i * count + j])641 .collect();642 (643 own.iter().cloned().min().unwrap_or(0),644 own.iter().cloned().max().unwrap_or(0),645 pairs.iter().cloned().min().unwrap_or(0),646 pairs.iter().cloned().max().unwrap_or(0),647 branch,648 )649}650651fn mark_law() {652 println!("THE MARK");653 println!(" m(x) = 4a(b eps x / a + arg(A + p e^(ix)))/pi runs 0 to 4 b eps while abs(p) < A, and is one to one for abs(p) < min(1, A)");654 println!(655 " a/b side reach abs(p) A tame self/curve pair/pair branch a(b-1) 2ab"656 );657 for &(a, b, inside, t) in &[658 (7usize, 3usize, true, 0.83),659 (7, 3, false, 0.83),660 (5, 2, true, 0.83),661 (11, 4, true, 0.83),662 (5, 1, false, 0.83),663 (13, 7, true, 0.83),664 (7, 5, true, 0.585),665 (7, 5, true, 0.615),666 (7, 5, true, 0.83),667 (4, 3, true, 0.49),668 (4, 3, true, 0.51),669 (9, 8, true, 0.83),670 ] {671 let w = wheel(a, b, inside);672 let (cs, _) = curves(&w, &carpet(), t);673 let nodes = census(&w, &cs, 1500);674 let (s0, s1, p0, p1, branch) = tally(&nodes, cs.len());675 let arm = cs.iter().map(|c| c.arm).fold(0.0_f64, f64::max);676 let tame = if arm < 1.0_f64.min(w.amp) {677 "yes"678 } else {679 "no "680 };681 let side = if inside { "in " } else { "out" };682 println!(683 " {a}/{b} {side} {t:6.3} {arm:6.4} {:6.4} {tame} {s0:2}..{s1:2} {p0:2}..{p1:2} {branch} {:5} {:4}",684 w.amp,685 a * (b - 1),686 2 * a * b687 );688 }689}690691fn carpet_table() {692 println!("THE ALIGNMENT");693 let w = wheel(7, 3, true);694 let (tlo, thi) = (0.05_f64, tame(&w).min(1.4999));695 let found = alignments(&w, tlo, thi, 6000);696 println!(" 7/3 inside, carpet fills, the window abs(p) < min(1, A) being the whole of abs(p) < 1 here: {} alignment reaches", found.len());697 let flats = found.iter().filter(|h| h.flat).count();698 let quads = found.iter().filter(|h| !h.flat && h.kinds == 4).count();699 println!(" the scan runs from reach {tlo} to {thi:.4}, which on the corner seats is abs(p) from {:.6} to {:.6}", 2.0 * tlo / 3.0, 2.0 * thi / 3.0);700 println!(" {} of them transversal and {flats} tangential, the tangential ones meeting an edge curve at its own apex", found.len() - flats);701 println!(" every transversal meeting carries four branches; {quads} of the {} carry four distinct curves and {} carry three, one curve bringing two branches as its own self crossing lands on the meeting", found.len() - flats, found.len() - flats - quads);702 println!(703 " reach abs(p) corner ring m_c m_e branches curves kind law residual"704 );705 for hit in &found {706 let arm_c = 2.0 * hit.reach / 3.0;707 let arm_e = 2.0_f64.sqrt() * hit.reach / 3.0;708 let turn_e = if hit.flat {709 if hit.mark_e == 0 {710 0.0711 } else {712 PI713 }714 } else {715 ((hit.ring * hit.ring - w.amp * w.amp - arm_e * arm_e) / (2.0 * w.amp * arm_e))716 .clamp(-1.0, 1.0)717 .acos()718 };719 let worst =720 residue(&w, arm_c, hit.turn, hit.mark_c).max(residue(&w, arm_e, turn_e, hit.mark_e));721 println!(722 " {:.12} {:.12} {:.6} {:3} {:3} {:6} {:5} {:5} {:.2e}",723 hit.reach,724 arm_c,725 hit.ring,726 hit.mark_c,727 hit.mark_e,728 hit.ends,729 hit.kinds,730 if hit.flat { "apex" } else { "cross" },731 worst732 );733 }734}735736fn node_of(cs: &[Curve], hit: &Align) -> Node {737 let (ends, _, _, seen) = meeting(cs, hit.mark_c, hit.mark_e, hit.flat);738 Node {739 ring: hit.ring,740 angle: hit.slot as f64 / 8.0,741 branches: ends,742 curves: seen,743 }744}745746fn pinch(w: &Wheel, t: f64) -> (Vec<Node>, usize, usize) {747 let (cs, _) = curves(w, &carpet(), t);748 let nodes = census(w, &cs, 8000);749 let deep = nodes.iter().filter(|n| n.branches > 2).count();750 let branch = nodes.iter().map(|n| n.branches).max().unwrap_or(0);751 (nodes, deep, branch)752}753754fn seat_hits(w: &Wheel, t: f64, node: &Node, samples: usize) -> (f64, f64, usize) {755 seat_gaps(w, t, node, probe(w, t, node.ring, node.angle, samples))756}757758fn seat_wide(w: &Wheel, t: f64, node: &Node, samples: usize) -> (f64, f64, usize) {759 seat_gaps(w, t, node, probe_wide(w, t, node.ring, node.angle, samples))760}761762fn seat_gaps(w: &Wheel, t: f64, node: &Node, gaps: Vec<f64>) -> (f64, f64, usize) {763 let (cs, owner) = curves(w, &carpet(), t);764 let _ = cs;765 let want: Vec<usize> = (0..owner.len())766 .filter(|&k| node.curves.contains(&owner[k]))767 .collect();768 let near = want.iter().map(|&k| gaps[k]).fold(0.0_f64, f64::max);769 let far = (0..gaps.len())770 .filter(|k| !want.contains(k))771 .map(|k| gaps[k])772 .fold(f64::INFINITY, f64::min);773 (near, far, want.len())774}775776fn witness() {777 println!("THE WITNESS");778 let w = wheel(7, 3, true);779 let (mut lo, mut hi) = (0.78_f64, 0.80_f64);780 let sign = defect(&w, lo, -6.0, -9.0);781 for _ in 0..200 {782 let mid = 0.5 * (lo + hi);783 if defect(&w, mid, -6.0, -9.0) * sign > 0.0 {784 lo = mid;785 } else {786 hi = mid;787 }788 }789 let t = 0.5 * (lo + hi);790 println!(791 " 7/3 inside, marks (-6, -9): reach {t:.12}, bracket width {:.1e}",792 hi - lo793 );794 println!(795 " abs(p) corner {:.12}, ring {:.12} wheel radii",796 2.0 * t / 3.0,797 ring(&w, 2.0 * t / 3.0, turn_for(&w, 2.0 * t / 3.0, -6.0))798 );799 let (nodes, deep, branch) = pinch(&w, t);800 println!(801 " census at the reach: {} classes, {deep} of them with {branch} branches, the rest 2",802 nodes.len()803 );804 let node = nodes.iter().find(|n| n.branches > 2).unwrap();805 let (cs, _) = curves(&w, &carpet(), t);806 let corners = node.curves.iter().filter(|&&k| cs[k].corner).count();807 println!(808 " the class at ring {:.9}, angle {:.9} of a turn over a: {corners} corner curves and {} edge curves",809 node.ring,810 node.angle,811 node.curves.len() - corners812 );813 let floor = f32_floor(w.b as f64 * node.ring);814 println!(" half an f32 step at that radius is {floor:.3e}, the floor any read of the f32 trace can reach");815 println!(" samples worst gap in f64 over the meeting seats the same read off the f32 trace least gap over the others");816 for samples in [2000usize, 8000, 32000, 128000] {817 let (wide, far, count) = seat_wide(&w, t, node, samples);818 let thin = seat_hits(&w, t, node, samples).0;819 println!(" {samples:7} {wide:.3e} over {count} seats {thin:.3e} {far:.3e}");820 }821 for shift in [-0.01_f64, 0.01] {822 let control = t + shift;823 let (plain, deep, branch) = pinch(&w, control);824 let (near, far, count) = seat_hits(&w, control, node, 8000);825 println!(826 " control reach {control:.6}: {} classes, {deep} past 2 branches, deepest {branch}; the same point now {near:.3e} from {count} seats, {far:.3e} from the rest",827 plain.len()828 );829 }830 let all = alignments(&w, 0.05, tame(&w).min(1.4999), 6000);831 let sample = curves(&w, &carpet(), 0.8).0;832 println!(" the tangential family, a corner pair meeting an edge curve at its own apex:");833 for hit in all.iter().filter(|h| h.flat) {834 let node = node_of(&sample, hit);835 let (wide, far, count) = seat_wide(&w, hit.reach, &node, 8000);836 println!(837 " reach {:.12}, marks ({}, {}), ring {:.9}: {count} seats within {wide:.3e} of the point in f64 and {:.3e} off the f32 trace, floor {:.3e}, the rest no nearer than {far:.3e}",838 hit.reach,839 hit.mark_c,840 hit.mark_e,841 hit.ring,842 seat_hits(&w, hit.reach, &node, 8000).0,843 f32_floor(w.b as f64 * hit.ring)844 );845 }846 let plain = pinch(&w, t - 0.01).0.len();847 println!(848 " the drop: {plain} classes and so {} nodes at a generic reach, {} classes and {} nodes at the alignment, five double points swallowed by each meeting",849 plain * w.a,850 nodes.len(),851 nodes.len() * w.a852 );853}854855fn gcd4(b: usize) -> usize {856 let (mut p, mut q) = (b, 4usize);857 while q > 0 {858 let r = p % q;859 p = q;860 q = r;861 }862 p863}864865#[derive(Clone)]866struct Band {867 tracks: usize,868 aligns: usize,869 bent: usize,870 held: usize,871 deep: usize,872 wide: f64,873 thin: f64,874 floor: f64,875 far: f64,876 breach: usize,877}878879fn survey() {880 println!("THE SWEEP");881 let mut rows: Vec<Band> = (0..8)882 .map(|_| Band {883 tracks: 0,884 aligns: 0,885 bent: 0,886 held: 0,887 deep: 0,888 wide: 0.0,889 thin: 0.0,890 floor: 0.0,891 far: f64::INFINITY,892 breach: 0,893 })894 .collect();895 let mut book: Vec<(usize, usize, bool, usize, usize)> = Vec::new();896 for &(a, b, inside) in &deck() {897 let w = wheel(a, b, inside);898 let cap = tame(&w).min(1.4999);899 let found = alignments(&w, 0.05, cap, 6000);900 book.push((b, a, inside, found.len(), 4 / gcd4(b)));901 let row = &mut rows[b - 1];902 row.tracks += 1;903 row.aligns += found.len();904 let shape = curves(&w, &carpet(), 0.5 * (0.05 + cap)).0;905 for hit in found.iter() {906 let node = if hit.flat {907 row.bent += 1;908 node_of(&shape, hit)909 } else {910 let (nodes, deep, branch) = pinch(&w, hit.reach);911 if branch < 3 {912 continue;913 }914 row.deep += deep;915 row.held += 1;916 let seen = nodes.iter().find(|n| n.branches > 2).unwrap();917 Node {918 ring: seen.ring,919 angle: seen.angle,920 branches: seen.branches,921 curves: seen.curves.clone(),922 }923 };924 let (wide, far, _) = seat_wide(&w, hit.reach, &node, 8000);925 row.wide = row.wide.max(wide);926 row.thin = row.thin.max(seat_hits(&w, hit.reach, &node, 8000).0);927 row.floor = row.floor.max(f32_floor(w.b as f64 * node.ring));928 if far.is_finite() {929 row.far = row.far.min(far);930 }931 }932 let control = match found.len() {933 0 => 0.7,934 1 => found[0].reach + 0.013,935 _ => 0.5 * (found[0].reach + found[1].reach),936 };937 if control > 0.05 && control < 1.4999 {938 let (_, _, branch) = pinch(&w, control);939 if branch > 2 {940 row.breach += 1;941 }942 }943 }944 println!(" b tracks alignments tangential held in the census classes past two branches worst gap f64 worst gap f32 f32 floor least other gap control breaches");945 for (k, row) in rows.iter().enumerate() {946 let other = if row.far.is_finite() {947 format!("{:.3e}", row.far)948 } else {949 "no other seat".to_string()950 };951 println!(952 " {} {:5} {:8} {:8} {:14} {:14} {:.3e} {:.3e} {:.3e} {other:13} {}",953 k + 1,954 row.tracks,955 row.aligns,956 row.bent,957 row.held,958 row.deep,959 row.wide,960 row.thin,961 row.floor,962 row.breach963 );964 }965 println!(" alignment reaches per track over the window abs(p) < min(1, A), a:count");966 for b in 1..=8usize {967 for inside in [true, false] {968 let line: Vec<String> = book969 .iter()970 .filter(|r| r.0 == b && r.2 == inside)971 .map(|r| format!("{}:{}", r.1, r.3))972 .collect();973 if line.is_empty() {974 continue;975 }976 println!(977 " b {b} {} classes {} {}",978 if inside { "in " } else { "out" },979 4 / gcd4(b),980 line.join(" ")981 );982 }983 }984 let total: usize = rows.iter().map(|r| r.aligns).sum();985 let breach: usize = rows.iter().map(|r| r.breach).sum();986 let held: usize = rows.iter().map(|r| r.held).sum();987 let bent: usize = rows.iter().map(|r| r.bent).sum();988 println!(" {total} alignment reaches over {} circle tracks, {bent} of them tangential; every one is read against mrlynum::spirograph::point and mrlynum::spirograph::trace and holds, {held} of the {} transversal ones also show a class past two branches in the census, and {breach} control reaches show one", deck().len(), total - bent);989 for (a, b, inside) in [990 (7usize, 3usize, true),991 (12, 7, true),992 (13, 7, true),993 (9, 5, true),994 (13, 8, true),995 (13, 7, false),996 ] {997 let probe = wheel(a, b, inside);998 let cap = tame(&probe).min(1.4999);999 let side = if inside { "inside" } else { "outside" };1000 println!(1001 " grid check on {a}/{b} {side}, reach up to {cap:.4}: {} at 6000, {} at 24000",1002 alignments(&probe, 0.05, cap, 6000).len(),1003 alignments(&probe, 0.05, cap, 24000).len()1004 );1005 }1006}10071008fn main() {1009 reduction();1010 mark_law();1011 carpet_table();1012 witness();1013 survey();1014}