main.rs
36.2 kB · rust · 1203 lines
1use mrlynum::spirograph::{frame, trace, track, Kind, Pencil};2use std::f64::consts::{PI, TAU};34// THE TRACK56fn rim(a: i64, b: i64, out: bool) -> i64 {7 if out {8 a + b9 } else {10 a - b11 }12}1314fn skew(a: i64, b: i64, out: bool) -> i64 {15 if out {16 a + 2 * b17 } else {18 a - 2 * b19 }20}2122fn gcd(mut x: i64, mut y: i64) -> i64 {23 while y != 0 {24 let t = x % y;25 x = y;26 y = t;27 }28 x29}3031// THE ROOTS3233fn bisect(f: &dyn Fn(f64) -> f64, mut lo: f64, mut hi: f64) -> f64 {34 let sign = f(lo).signum();35 for _ in 0..200 {36 let mid = 0.5 * (lo + hi);37 if f(mid).signum() == sign {38 lo = mid;39 } else {40 hi = mid;41 }42 }43 0.5 * (lo + hi)44}4546fn sweep(f: &dyn Fn(f64) -> f64, lo: f64, hi: f64, steps: usize) -> Vec<f64> {47 let mut out = Vec::new();48 let mut prev = f(lo);49 let mut left = lo;50 for k in 1..=steps {51 let x = lo + (hi - lo) * k as f64 / steps as f64;52 let cur = f(x);53 if prev * cur < 0.0 {54 out.push(bisect(f, left, x));55 }56 prev = cur;57 left = x;58 }59 out60}6162// THE CROSSING EQUATION6364fn arm_roots(a: i64, b: i64, out: bool, lam: f64, eps: f64, steps: usize) -> Vec<f64> {65 let (r, bf) = (rim(a, b, out) as f64, b as f64);66 let f = |d: f64| r * (bf * d).sin() - eps * lam * bf * (r * d).sin();67 sweep(&f, 1e-7, PI - 1e-7, steps)68}6970fn arms(a: i64, b: i64, out: bool, lam: f64, steps: usize) -> usize {71 arm_roots(a, b, out, lam, 1.0, steps).len() + arm_roots(a, b, out, lam, -1.0, steps).len()72}7374fn crossings(a: i64, b: i64, out: bool, lam: f64, steps: usize) -> usize {75 a as usize * arms(a, b, out, lam, steps) / 276}7778// THE TANGENCY7980fn tangent_angles(a: i64, b: i64, out: bool, steps: usize) -> Vec<f64> {81 let m = skew(a, b, out);82 let (af, mf) = (a as f64, m as f64);83 let g = |d: f64| af * (mf * d).sin() - mf * (af * d).sin();84 let lo = 0.4 / a.max(m.abs()) as f64;85 let mut out_angles = vec![0.0];86 for d in sweep(&g, lo, PI / 2.0, steps) {87 if PI / 2.0 - d > 1e-6 {88 out_angles.push(d);89 out_angles.push(PI - d);90 }91 }92 if g(PI / 2.0).abs() < 1e-9 {93 out_angles.push(PI / 2.0);94 }95 out_angles96}9798fn reach_at(a: i64, b: i64, out: bool, d: f64) -> f64 {99 let (r, bf) = (rim(a, b, out) as f64, b as f64);100 let (s, c) = ((r * d).sin(), (r * d).cos());101 if s.abs() > c.abs() {102 (r * (bf * d).sin() / (bf * s)).abs()103 } else {104 ((bf * d).cos() / c).abs()105 }106}107108fn steps_of(a: i64, b: i64, out: bool, steps: usize) -> Vec<f64> {109 let mut out_reach: Vec<f64> = tangent_angles(a, b, out, steps)110 .into_iter()111 .map(|d| reach_at(a, b, out, d))112 .collect();113 out_reach.sort_by(|x, y| x.partial_cmp(y).unwrap());114 out_reach115}116117fn law(a: i64, b: i64, out: bool, lam: f64, steps: &[f64]) -> i64 {118 let sign = skew(a, b, out).signum();119 let past = steps.iter().filter(|&&t| t < lam).count() as i64;120 a * (b - 1) + sign * a * past121}122123// THE SIGNS124125fn parity(e: i64) -> i64 {126 if e % 2 == 0 {127 1128 } else {129 -1130 }131}132133fn turns(b: i64, r: i64) -> i64 {134 let flip = (r - b).signum();135 let mut marks: Vec<(i64, i64, i64)> = Vec::new();136 for j in 1..r {137 marks.push((j, r, flip * parity(j + 1 + b * j / r)));138 }139 for k in 1..b {140 marks.push((k, b, flip * parity(k + r * k / b)));141 }142 marks.sort_by_key(|&(n, d, _)| (n * (b * r / d), d));143 let mut walk = vec![1_i64];144 walk.extend(marks.iter().map(|&(_, _, s)| s));145 1 + walk.windows(2).filter(|w| w[0] != w[1]).count() as i64146}147148// THE POLYLINE149150fn seat(p: (f64, f64)) -> Pencil {151 Pencil {152 x: p.0,153 y: p.1,154 seat: (0, 0),155 kind: Kind::Fill,156 }157}158159fn draw(160 a: i64,161 b: i64,162 out: bool,163 pens: &[(f64, f64)],164 samples: usize,165) -> (Vec<Vec<f64>>, [f64; 4]) {166 let kind = if out { "out" } else { "in" };167 let path = track(kind, a as usize, b as usize, 4, 1).unwrap();168 let pencils: Vec<Pencil> = pens.iter().map(|&p| seat(p)).collect();169 let flat = trace(&path, &pencils, samples).unwrap();170 let curves = pencils171 .iter()172 .enumerate()173 .map(|(k, _)| {174 flat[k * samples * 2..(k + 1) * samples * 2]175 .iter()176 .map(|&v| v as f64)177 .collect()178 })179 .collect();180 (curves, frame(&path, &pencils))181}182183fn hits(184 ax: f64,185 ay: f64,186 bx: f64,187 by: f64,188 cx: f64,189 cy: f64,190 dx: f64,191 dy: f64,192) -> Option<(f64, f64)> {193 let (ux, uy) = (bx - ax, by - ay);194 let (vx, vy) = (dx - cx, dy - cy);195 let det = ux * vy - uy * vx;196 if det == 0.0 {197 return None;198 }199 let (wx, wy) = (cx - ax, cy - ay);200 let t = (wx * vy - wy * vx) / det;201 let s = (wx * uy - wy * ux) / det;202 if (0.0..1.0).contains(&t) && (0.0..1.0).contains(&s) {203 Some((ax + t * ux, ay + t * uy))204 } else {205 None206 }207}208209fn nodes(210 curves: &[Vec<f64>],211 box_: [f64; 4],212 grid: usize,213 axes: i64,214 tilt: f64,215) -> (usize, usize, f64, usize, usize) {216 let mut segs: Vec<[f64; 4]> = Vec::new();217 let mut owner: Vec<u32> = Vec::new();218 let mut index: Vec<u32> = Vec::new();219 for (k, curve) in curves.iter().enumerate() {220 let n = curve.len() / 2 - 1;221 for i in 0..n {222 segs.push([223 curve[2 * i],224 curve[2 * i + 1],225 curve[2 * i + 2],226 curve[2 * i + 3],227 ]);228 owner.push(k as u32);229 index.push(i as u32);230 }231 }232 let last: Vec<u32> = curves.iter().map(|c| (c.len() / 2 - 2) as u32).collect();233 let (wx, wy) = (box_[2] - box_[0], box_[3] - box_[1]);234 let mut bins: Vec<Vec<u32>> = vec![Vec::new(); grid * grid];235 let put = |v: f64, lo: f64, w: f64| {236 (((v - lo) / w * grid as f64) as isize).clamp(0, grid as isize - 1) as usize237 };238 for (id, s) in segs.iter().enumerate() {239 let (x0, x1) = (240 put(s[0].min(s[2]), box_[0], wx),241 put(s[0].max(s[2]), box_[0], wx),242 );243 let (y0, y1) = (244 put(s[1].min(s[3]), box_[1], wy),245 put(s[1].max(s[3]), box_[1], wy),246 );247 for cx in x0..=x1 {248 for cy in y0..=y1 {249 bins[cx * grid + cy].push(id as u32);250 }251 }252 }253 let mut seen: std::collections::HashSet<u64> = std::collections::HashSet::new();254 let (mut own, mut between, mut off) = (0, 0, 0.0_f64);255 let (mut mine, mut yours): (Vec<(f64, f64)>, Vec<(f64, f64)>) = (Vec::new(), Vec::new());256 for bin in &bins {257 for u in 0..bin.len() {258 for v in u + 1..bin.len() {259 let (i, j) = (bin[u].min(bin[v]) as usize, bin[u].max(bin[v]) as usize);260 if owner[i] == owner[j] {261 let (p, q) = (index[i], index[j]);262 if q - p == 1 || (p == 0 && q == last[owner[i] as usize]) {263 continue;264 }265 }266 if !seen.insert((i as u64) << 32 | j as u64) {267 continue;268 }269 let (s, t) = (segs[i], segs[j]);270 if let Some(spot) = hits(s[0], s[1], s[2], s[3], t[0], t[1], t[2], t[3]) {271 if owner[i] == owner[j] {272 own += 1;273 off = off.max(stray(spot.0, spot.1, axes, tilt));274 mine.push(spot);275 } else {276 between += 1;277 yours.push(spot);278 }279 }280 }281 }282 }283 let tol = (wx.max(wy)) * 2e-3;284 (own, between, off, thin(&mine, tol), thin(&yours, tol))285}286287fn thin(spots: &[(f64, f64)], tol: f64) -> usize {288 let mut kept: Vec<(f64, f64)> = Vec::new();289 for &spot in spots {290 if !kept291 .iter()292 .any(|k: &(f64, f64)| (k.0 - spot.0).hypot(k.1 - spot.1) < tol)293 {294 kept.push(spot);295 }296 }297 kept.len()298}299300fn stray(x: f64, y: f64, axes: i64, tilt: f64) -> f64 {301 if axes <= 0 {302 return 0.0;303 }304 let (r, theta) = (x.hypot(y), y.atan2(x));305 (0..axes)306 .map(|k| (r * (theta - tilt - PI * k as f64 / axes as f64).sin()).abs())307 .fold(f64::MAX, f64::min)308}309310// THE SURDS311312#[derive(Clone, Copy)]313struct Surd {314 a: i128,315 b: i128,316 den: i128,317 root: i128,318}319320fn square_part(mut d: i128) -> (i128, i128) {321 let (mut f, mut k) = (1_i128, 2_i128);322 while k * k <= d {323 while d % (k * k) == 0 {324 d /= k * k;325 f *= k;326 }327 k += 1;328 }329 (f, d)330}331332impl Surd {333 fn tidy(self) -> Surd {334 let g = gcd(335 gcd(336 self.a.unsigned_abs().min(i64::MAX as u128) as i64,337 self.b.unsigned_abs().min(i64::MAX as u128) as i64,338 ),339 self.den.unsigned_abs().min(i64::MAX as u128) as i64,340 )341 .max(1) as i128;342 let sign = if self.den < 0 { -1 } else { 1 };343 Surd {344 a: sign * self.a / g,345 b: sign * self.b / g,346 den: sign * self.den / g,347 root: self.root,348 }349 }350351 fn plus(self, other: Surd) -> Option<Surd> {352 Some(353 Surd {354 a: self355 .a356 .checked_mul(other.den)?357 .checked_add(other.a.checked_mul(self.den)?)?,358 b: self359 .b360 .checked_mul(other.den)?361 .checked_add(other.b.checked_mul(self.den)?)?,362 den: self.den.checked_mul(other.den)?,363 root: self.root,364 }365 .tidy(),366 )367 }368369 fn times(self, other: Surd) -> Option<Surd> {370 Some(371 Surd {372 a: self373 .a374 .checked_mul(other.a)?375 .checked_add(self.b.checked_mul(other.b)?.checked_mul(self.root)?)?,376 b: self377 .a378 .checked_mul(other.b)?379 .checked_add(self.b.checked_mul(other.a)?)?,380 den: self.den.checked_mul(other.den)?,381 root: self.root,382 }383 .tidy(),384 )385 }386387 fn over(self, other: Surd) -> Option<Surd> {388 let base = other389 .a390 .checked_mul(other.a)?391 .checked_sub(other.b.checked_mul(other.b)?.checked_mul(other.root)?)?;392 if base == 0 {393 return None;394 }395 self.times(Surd {396 a: other.den.checked_mul(other.a)?,397 b: -other.den.checked_mul(other.b)?,398 den: base,399 root: self.root,400 })401 }402403 fn value(self) -> f64 {404 (self.a as f64 + self.b as f64 * (self.root as f64).sqrt()) / self.den as f64405 }406407 fn text(self) -> String {408 if self.b == 0 {409 return format!("{}/{}", self.a, self.den);410 }411 format!(412 "({} {} {} sqrt {})/{}",413 self.a,414 if self.b < 0 { "-" } else { "+" },415 self.b.abs(),416 self.root,417 self.den418 )419 }420}421422fn poly_surd(p: &[i128], u: Surd) -> Option<Surd> {423 let mut acc = Surd {424 a: 0,425 b: 0,426 den: 1,427 root: u.root,428 };429 for &c in p.iter().rev() {430 acc = acc.times(u)?.plus(Surd {431 a: c,432 b: 0,433 den: 1,434 root: u.root,435 })?;436 }437 Some(acc)438}439440fn quadratic_roots(q: &[i128]) -> Vec<Surd> {441 if q.len() != 3 {442 return Vec::new();443 }444 let (c, b, a) = (q[0], q[1], q[2]);445 let disc = b * b - 4 * a * c;446 if disc < 0 {447 return Vec::new();448 }449 let (f, d) = square_part(disc);450 [1_i128, -1]451 .iter()452 .map(|&s| {453 Surd {454 a: -b,455 b: s * f,456 den: 2 * a,457 root: d,458 }459 .tidy()460 })461 .collect()462}463464// THE PAIRS465466fn pair_arms(a: i64, b: i64, out: bool, p: (f64, f64), q: (f64, f64), steps: usize) -> usize {467 let (r, bf) = (rim(a, b, out) as f64, b as f64);468 let sw = if out { 1.0 } else { -1.0 };469 let mid = (0.5 * (p.0 + q.0), 0.5 * (p.1 + q.1));470 let gap = (0.5 * (p.0 - q.0), 0.5 * (p.1 - q.1));471 let big = mid.0 * mid.0 + mid.1 * mid.1;472 let small = gap.0 * gap.0 + gap.1 * gap.1;473 let tilt = mid.1 * gap.0 - mid.0 * gap.1;474 let f = |d: f64| {475 let (s, c) = ((r * d).sin(), (r * d).cos());476 let load = big * s * s + small * c * c - sw * tilt * (2.0 * r * d).sin();477 let w = (bf * d).sin();478 r * r * w * w - bf * bf * load479 };480 sweep(&f, 0.0, TAU, steps).len()481}482483fn pair_count(a: i64, b: i64, out: bool, p: (f64, f64), q: (f64, f64), steps: usize) -> usize {484 a as usize * pair_arms(a, b, out, p, q, steps) / 2485}486487// THE ALGEBRA488489fn poly_add(p: &[i128], q: &[i128], sp: i128, sq: i128) -> Vec<i128> {490 let n = p.len().max(q.len());491 (0..n)492 .map(|k| sp * p.get(k).copied().unwrap_or(0) + sq * q.get(k).copied().unwrap_or(0))493 .collect()494}495496fn poly_shift(p: &[i128]) -> Vec<i128> {497 let mut out = vec![0];498 out.extend_from_slice(p);499 out500}501502fn poly_trim(mut p: Vec<i128>) -> Vec<i128> {503 while p.len() > 1 && *p.last().unwrap() == 0 {504 p.pop();505 }506 p507}508509fn angle_polys(top: i64) -> Vec<(Vec<i128>, Vec<i128>)> {510 let mut out = vec![(vec![0_i128], vec![1_i128])];511 for k in 0..top {512 let (s, c) = (out[k as usize].0.clone(), out[k as usize].1.clone());513 let next = if k % 2 == 0 {514 let sin = poly_add(&c, &poly_add(&s, &poly_shift(&s), 1, -1), 1, 1);515 let cos = poly_add(&c, &poly_shift(&s), 1, -1);516 (poly_trim(sin), poly_trim(cos))517 } else {518 let sin = poly_add(&s, &c, 1, 1);519 let cos = poly_add(520 &poly_add(&c, &poly_shift(&c), 1, -1),521 &poly_shift(&s),522 1,523 -1,524 );525 (poly_trim(sin), poly_trim(cos))526 };527 out.push(next);528 }529 out530}531532fn poly_at(p: &[i128], u: f64) -> f64 {533 p.iter().rev().fold(0.0, |acc, &c| acc * u + c as f64)534}535536fn poly_text(p: &[i128]) -> String {537 let mut parts = Vec::new();538 for (k, &c) in p.iter().enumerate().rev() {539 if c == 0 {540 continue;541 }542 parts.push(match k {543 0 => format!("{c}"),544 1 => format!("{c} u"),545 _ => format!("{c} u^{k}"),546 });547 }548 if parts.is_empty() {549 "0".to_string()550 } else {551 parts.join(" + ").replace("+ -", "- ")552 }553}554555fn poly_mul(p: &[i128], q: &[i128]) -> Vec<i128> {556 let mut out = vec![0_i128; p.len() + q.len() - 1];557 for (i, &x) in p.iter().enumerate() {558 for (j, &y) in q.iter().enumerate() {559 out[i + j] += x * y;560 }561 }562 poly_trim(out)563}564565fn cos_square(tab: &[(Vec<i128>, Vec<i128>)], k: i64) -> Vec<i128> {566 let part = &tab[k as usize].1;567 let sq = poly_mul(part, part);568 if k % 2 == 0 {569 sq570 } else {571 poly_mul(&[1, -1], &sq)572 }573}574575fn reach_square(a: i64, b: i64, out: bool) -> (Vec<i128>, Vec<i128>) {576 let r = rim(a, b, out);577 let tab = angle_polys(b.max(r));578 (cos_square(&tab, b), cos_square(&tab, r))579}580581fn rational_root(p: &[i128], u: f64) -> Option<(i128, i128)> {582 for den in 1_i128..=120 {583 let num = (u * den as f64).round() as i128;584 if num < 0 || gcd(num as i64, den as i64) != 1 {585 continue;586 }587 if (num as f64 / den as f64 - u).abs() > 1e-9 {588 continue;589 }590 let deg = p.len() - 1;591 let mut acc = 0_i128;592 let mut ok = true;593 for (k, &c) in p.iter().enumerate() {594 match num595 .checked_pow(k as u32)596 .and_then(|x| {597 den.checked_pow((deg - k) as u32)598 .and_then(|y| x.checked_mul(y))599 })600 .and_then(|z| z.checked_mul(c))601 .and_then(|z| acc.checked_add(z))602 {603 Some(v) => acc = v,604 None => ok = false,605 }606 }607 if ok && acc == 0 {608 return Some((num, den));609 }610 }611 None612}613614fn poly_ratio(p: &[i128], num: i128, den: i128) -> Option<(i128, i128)> {615 let deg = p.len() - 1;616 let mut acc = 0_i128;617 for (k, &c) in p.iter().enumerate() {618 let term = num619 .checked_pow(k as u32)?620 .checked_mul(den.checked_pow((deg - k) as u32)?)?621 .checked_mul(c)?;622 acc = acc.checked_add(term)?;623 }624 Some((acc, den.checked_pow(deg as u32)?))625}626627fn tangency_poly(a: i64, m: i64) -> Vec<i128> {628 let m = m.abs();629 let tab = angle_polys(a.max(m));630 let (pa, pm) = (&tab[a as usize].0, &tab[m as usize].0);631 poly_trim(poly_add(pm, pa, a as i128, -(m as i128)))632}633634fn cases() -> Vec<(i64, i64, bool)> {635 let mut out = Vec::new();636 for b in 1..=6 {637 for a in b + 1..=11 {638 if gcd(a, b) != 1 {639 continue;640 }641 if skew(a, b, false) != 0 {642 out.push((a, b, false));643 }644 out.push((a, b, true));645 }646 }647 out648}649650fn side(out: bool) -> f64 {651 if out {652 1.0653 } else {654 -1.0655 }656}657658fn tag(out: bool) -> &'static str {659 if out {660 "out"661 } else {662 "in "663 }664}665666fn walk(a: i64, b: i64, out: bool, steps: &[f64]) -> String {667 let mut text = format!("{}", a * (b - 1));668 let mut k = 0;669 while k < steps.len() {670 let mut n = 1;671 while k + n < steps.len() && (steps[k + n] - steps[k]).abs() < 1e-9 {672 n += 1;673 }674 k += n;675 text += &format!(676 " |{:.6}| {}",677 steps[k - 1],678 law(a, b, out, steps[k - 1] + 1e-9, steps)679 );680 }681 text682}683684fn exact_reach(a: i64, b: i64, out: bool, num: i128, den: i128) -> Option<(i128, i128)> {685 let (top, bot) = reach_square(a, b, out);686 let (tn, td) = poly_ratio(&top, num, den)?;687 let (bn, bd) = poly_ratio(&bot, num, den)?;688 if bn == 0 {689 return None;690 }691 let (mut p, mut q) = (tn.checked_mul(bd)?, bn.checked_mul(td)?);692 if q < 0 {693 p = -p;694 q = -q;695 }696 let g = gcd(p.unsigned_abs().min(i64::MAX as u128) as i64, q as i64).max(1) as i128;697 Some((p / g, q / g))698}699700fn axes() {701 println!(702 "THE AXES every self crossing of one trochoid sits on one of the a mirror lines k pi / a"703 );704 for (a, b, out) in [705 (5, 1, false),706 (7, 2, false),707 (8, 3, false),708 (5, 2, true),709 (7, 3, true),710 ] {711 for lam in [1.4, 2.6, 3.9] {712 for alpha in [0.0_f64, 0.3] {713 let seat = (lam * alpha.cos(), lam * alpha.sin());714 let tilt = -side(out) * b as f64 * alpha / a as f64;715 let (curves, box_) = draw(a, b, out, &[seat], 24001);716 let (own, _, off, _, _) = nodes(&curves, box_, 260, a, tilt);717 let size = (box_[2] - box_[0]).max(box_[3] - box_[1]);718 println!(719 " {} {a}/{b} reach {lam} seat angle {alpha} nodes {own} worst stray {:.3e} of frame {size:.3} relative {:.2e}",720 tag(out),721 off,722 off / size723 );724 }725 }726 }727}728729fn staircase() {730 println!("THE STAIRCASE count on each step, tangency reaches between, angles counted with multiplicity");731 for &(a, b, out) in cases().iter().filter(|c| c.0 <= 8) {732 let steps = steps_of(a, b, out, 400_000);733 let m = skew(a, b, out);734 println!(735 " {} {a}/{b} m {m} angles {} of {} {}",736 tag(out),737 steps.len(),738 turns(b, rim(a, b, out)),739 walk(a, b, out, &steps)740 );741 }742 println!("THE ALGEBRA u = sin^2 of the tangency angle, Q(u) = 0 its integer equation, lam^2 rational in u");743 for (a, b, out) in [744 (5, 1, false),745 (7, 2, false),746 (8, 3, false),747 (11, 4, false),748 (3, 1, true),749 (5, 1, true),750 (3, 2, true),751 (5, 2, true),752 ] {753 let m = skew(a, b, out);754 let p = tangency_poly(a, m);755 let q: Vec<i128> = p[1..].to_vec();756 let stem = p[0];757 let (top, bot) = reach_square(a, b, out);758 println!(759 " {} {a}/{b} P(0) = {stem} Q(u) = {} lam^2 = ({}) / ({})",760 tag(out),761 poly_text(&q),762 poly_text(&top),763 poly_text(&bot)764 );765 for d in tangent_angles(a, b, out, 400_000) {766 if d == 0.0 || d > PI / 2.0 + 1e-9 {767 continue;768 }769 let u = d.sin() * d.sin();770 let lam = reach_at(a, b, out, d);771 let ring = a % 2 == 0 && (u - 1.0).abs() < 1e-12;772 let root = if ring {773 Some((1, 1))774 } else {775 rational_root(&q, u)776 };777 let exact = if ring {778 let r = rim(a, b, out) as i128;779 Some((r * r, (b * b) as i128))780 } else {781 root.and_then(|(n, dn)| exact_reach(a, b, out, n, dn))782 };783 let mut surd = String::new();784 if exact.is_none() {785 for pick in quadratic_roots(&q) {786 if (pick.value() - u).abs() > 1e-9 {787 continue;788 }789 if let Some(v) = poly_surd(&top, pick)790 .and_then(|t| poly_surd(&bot, pick).and_then(|w| t.over(w)))791 {792 surd = format!("= {} check {:.12}", v.text(), v.value());793 }794 }795 }796 println!(797 " u {u:.12} {} Q(u) {:.2e} lam {lam:.12} lam^2 {:.12} {}",798 match root {799 Some((n, dn)) => format!("= {n}/{dn}"),800 None => "irrational".to_string(),801 },802 if ring {803 0.0804 } else {805 poly_at(&q, u) / q.iter().map(|c| c.unsigned_abs() as f64).sum::<f64>()806 },807 lam * lam,808 match exact {809 Some((n, dn)) => format!("= {n}/{dn}"),810 None => surd,811 }812 );813 }814 }815}816817fn signs(top: i64) {818 println!(819 "THE SIGNS the threshold count from the exact critical values of the crossing integral"820 );821 let (mut pairs, mut wrong) = (0_usize, 0_usize);822 for r in 1..=top {823 for b in 1..=top {824 if b == r || gcd(b, r) != 1 {825 continue;826 }827 pairs += 1;828 if turns(b, r) != (r - b).abs() {829 wrong += 1;830 println!(831 " SIGNS b {b} rho {r} count {} want {}",832 turns(b, r),833 (r - b).abs()834 );835 }836 }837 }838 println!(" frequency pairs {pairs} to {top} failures {wrong}");839}840841fn census(top: i64) {842 println!("THE CENSUS thresholds counted against min(|m|, a), the closed form against the crossing equation");843 let (mut cases, mut wrong, mut ends) = (0_usize, 0_usize, 0_usize);844 let mut fractions = 0_usize;845 for a in 2..=top {846 for b in 1..a {847 if gcd(a, b) != 1 {848 continue;849 }850 fractions += 1;851 for out in [false, true] {852 let m = skew(a, b, out);853 if m == 0 {854 continue;855 }856 cases += 1;857 let steps = steps_of(a, b, out, 600_000);858 if steps.len() as i64 != turns(b, rim(a, b, out)) {859 wrong += 1;860 println!(861 " STEPS {} {a}/{b} m {m} found {} want {}",862 tag(out),863 steps.len(),864 turns(b, rim(a, b, out))865 );866 }867 if law(a, b, out, steps.last().unwrap() * 4.0 + 4.0, &steps)868 != a * (rim(a, b, out) - 1)869 {870 ends += 1;871 println!(" LIMIT {} {a}/{b} m {m}", tag(out));872 }873 let mut marks = vec![steps[0] * 0.5];874 for k in 1..steps.len() {875 if steps[k] - steps[k - 1] > 1e-6 {876 marks.push((steps[k] * steps[k - 1]).sqrt());877 }878 }879 marks.push(steps.last().unwrap() * 1.7 + 0.3);880 for lam in marks {881 let root = crossings(a, b, out, lam, 60_000) as i64;882 let want = law(a, b, out, lam, &steps);883 if root != want {884 wrong += 1;885 println!(886 " LAW {} {a}/{b} reach {lam:.6} equation {root} law {want}",887 tag(out)888 );889 }890 }891 }892 }893 }894 println!(" fractions {fractions} cases {cases} step-count or law failures {wrong} limit failures {ends}");895}896897fn spot(a: i64, b: i64, out: bool, lam: f64, phi: f64) -> (f64, f64) {898 let (r, bf, sw) = (rim(a, b, out) as f64, b as f64, side(out));899 let (u, v) = (bf * phi, sw * r * phi);900 (901 r * u.cos() + lam * bf * v.cos(),902 r * u.sin() + lam * bf * v.sin(),903 )904}905906fn mirror(top: i64) {907 println!("THE RECIPROCITY swapping the wheel and the rim frequency inverts every threshold");908 let (mut pairs, mut worst) = (0_usize, 0.0_f64);909 for b in 1..top {910 for r in 1..top {911 if b == r || b + r > top || gcd(b, r) != 1 {912 continue;913 }914 pairs += 1;915 let here = steps_of(b + r, b, false, 400_000);916 let there = steps_of(b + r, r, false, 400_000);917 for (x, y) in here.iter().zip(there.iter().rev()) {918 worst = worst.max((x * y - 1.0).abs());919 }920 }921 }922 println!(" wheel and rim swaps {pairs} to {top} worst deviation of the product from one {worst:.2e}");923}924925fn contact() {926 println!("THE CONTACT on a tangency reach the meetings are the transversal ones plus a/2 per tangency angle, each a tacnode");927 for (a, b, out) in [928 (5, 1, false),929 (7, 2, false),930 (11, 4, false),931 (3, 1, true),932 (5, 2, true),933 ] {934 let steps = steps_of(a, b, out, 400_000);935 let angles = tangent_angles(a, b, out, 400_000);936 let mut hits: Vec<f64> = Vec::new();937 for &t in steps.iter().filter(|&&t| t > 1.0 + 1e-9) {938 if !hits.iter().any(|h| (h - t).abs() < 1e-7) {939 hits.push(t);940 }941 }942 for hit in hits {943 let here: Vec<f64> = angles944 .iter()945 .copied()946 .filter(|&d| d > 1e-9 && (reach_at(a, b, out, d) - hit).abs() < 1e-7)947 .collect();948 let simple = arms(a, b, out, hit, 400_000);949 let under = law(a, b, out, hit * (1.0 - 1e-7), &steps);950 let (r, bf) = (rim(a, b, out) as f64, b as f64);951 let d0 = here[0];952 let eps = if (r * (bf * d0).sin() - hit * bf * (r * d0).sin()).abs() < 1e-6 {953 1.0954 } else {955 -1.0956 };957 let sigma = if eps * side(out) > 0.0 {958 PI / a as f64959 } else {960 0.0961 };962 let (p0, p1) = (963 spot(a, b, out, hit, sigma + d0),964 spot(a, b, out, hit, sigma - d0),965 );966 println!(967 " {} {a}/{b} reach {hit:.9} angles {} simple roots {simple} meetings {} under {under} over {} contact gap {:.2e} at radius {:.6}",968 tag(out),969 here.len(),970 a as usize * (simple + here.len()) / 2,971 law(a, b, out, hit * (1.0 + 1e-7), &steps),972 (p0.0 - p1.0).hypot(p0.1 - p1.1),973 p0.0.hypot(p0.1)974 );975 }976 }977 println!("THE BIRTHS a root is born at each end of (0, pi) as the reach passes 1, so the jump of a is two births of a/2");978 for (a, b, out) in [(3, 1, false), (5, 1, false), (7, 2, false), (3, 1, true)] {979 for lam in [0.98, 1.02] {980 let mut ends = Vec::new();981 for eps in [1.0_f64, -1.0] {982 let roots = arm_roots(a, b, out, lam, eps, 400_000);983 ends.push(format!(984 "e{}:{}",985 if eps > 0.0 { "+" } else { "-" },986 roots987 .iter()988 .map(|d| format!("{d:.4}"))989 .collect::<Vec<String>>()990 .join(",")991 ));992 }993 println!(994 " {} {a}/{b} reach {lam} roots {} count {}",995 tag(out),996 ends.join(" "),997 crossings(a, b, out, lam, 400_000)998 );999 }1000 }1001}10021003fn centre() {1004 println!("THE CENTRE at reach rho/b the curve runs through the origin with a branches and the point count drops by C(a,2) - 1");1005 for (a, b, out) in [1006 (3, 1, false),1007 (5, 1, false),1008 (5, 2, false),1009 (5, 3, false),1010 (5, 4, false),1011 (7, 2, false),1012 (7, 4, false),1013 (3, 1, true),1014 (3, 2, true),1015 (5, 2, true),1016 (4, 1, false),1017 ] {1018 let steps = steps_of(a, b, out, 400_000);1019 let hub = rim(a, b, out) as f64 / b as f64;1020 let mut row = Vec::new();1021 for (name, lam) in [("under", hub * 0.97), ("hub", hub), ("over", hub * 1.03)] {1022 let (curves, box_) = draw(a, b, out, &[(lam, 0.0)], 24001);1023 let (own, _, _, spots, _) = nodes(&curves, box_, 260, 0, 0.0);1024 row.push(format!(1025 "{name} pairs {own} points {spots} law {}",1026 law(a, b, out, lam, &steps)1027 ));1028 }1029 println!(1030 " {} {a}/{b} rho/b {hub:.6} {} hub law less C(a,2)-1 {}{}",1031 tag(out),1032 row.join(" "),1033 law(a, b, out, hub, &steps) - a * (a - 1) / 2 + 1,1034 if a % 2 == 0 {1035 " even a, the hub is a threshold too"1036 } else {1037 ""1038 }1039 );1040 }1041}10421043fn hunt() {1044 println!("THE SWEEP reach 0.5 to 4, trace against the closed form, two sample counts");1045 let (mut tests, mut wrong, mut worst_gap) = (0_usize, 0_usize, 0.0_f64);1046 let mut blur = 0.0_f64;1047 for (a, b, out) in cases() {1048 let steps = steps_of(a, b, out, 400_000);1049 let inband: Vec<f64> = steps1050 .iter()1051 .copied()1052 .filter(|&t| (0.4987..=4.0014).contains(&t))1053 .collect();1054 let mut seen: Vec<(f64, f64, i64, i64)> = Vec::new();1055 let mut prev: Option<(f64, i64)> = None;1056 for k in 0..=202 {1057 let lam = 0.4987 + 0.01734 * k as f64;1058 let want = law(a, b, out, lam, &steps);1059 let root = crossings(a, b, out, lam, 40_000) as i64;1060 if root != want {1061 wrong += 1;1062 println!(1063 " ROOTS {} {a}/{b} reach {lam:.4} equation {root} law {want}",1064 tag(out)1065 );1066 }1067 let mut got = 0;1068 for samples in [4001, 12001] {1069 let (curves, box_) = draw(a, b, out, &[(lam, 0.0)], samples);1070 got = nodes(&curves, box_, 220, 0, 0.0).0 as i64;1071 tests += 1;1072 if got != want {1073 wrong += 1;1074 let gap = inband1075 .iter()1076 .map(|t| (t - lam).abs())1077 .fold(f64::MAX, f64::min);1078 blur = blur.max(gap);1079 println!(" MISS {} {a}/{b} reach {lam:.4} samples {samples} trace {got} law {want} gap {gap:.4}", tag(out));1080 }1081 }1082 if let Some((last, count)) = prev {1083 if got != count {1084 seen.push((last, lam, count, got));1085 worst_gap = worst_gap.max(lam - last);1086 }1087 }1088 prev = Some((lam, got));1089 }1090 let jumps: Vec<String> = seen1091 .iter()1092 .map(|(x, y, c0, c1)| {1093 let inside = inband.iter().filter(|t| **t > *x && **t <= *y).count();1094 format!("({x:.3},{y:.3}) {c0}->{c1} holds {inside}")1095 })1096 .collect();1097 let caught: usize = seen1098 .iter()1099 .map(|(x, y, _, _)| inband.iter().filter(|t| **t > *x && **t <= *y).count())1100 .sum();1101 println!(1102 " {} {a}/{b} predicted {:?} jumps {} caught {caught} of {}",1103 tag(out),1104 inband1105 .iter()1106 .map(|t| (t * 1e4).round() / 1e4)1107 .collect::<Vec<f64>>(),1108 jumps.join(" "),1109 inband.len()1110 );1111 }1112 println!(" tests {tests} disagreements {wrong} widest jump bracket {worst_gap:.4} widest blur {blur:.4}");1113}11141115fn pairs() {1116 println!("THE PAIRS two pencils on one wheel, crossings between the two curves");1117 for (a, b, out) in [(5, 1, false), (7, 2, false), (8, 3, false), (5, 2, true)] {1118 for lam in [0.6, 1.3, 2.4, 3.7] {1119 let p = (lam, 0.0);1120 for (name, q) in [1121 ("even", (lam * 0.6, lam * 0.8)),1122 ("short", (lam * 0.45, lam * 0.35)),1123 ("centre", (0.0, 0.0)),1124 ] {1125 let (curves, box_) = draw(a, b, out, &[p, q], 12001);1126 let (_, between, _, _, spots) = nodes(&curves, box_, 220, 0, 0.0);1127 println!(1128 " {} {a}/{b} reach {lam} seat {name} trace pairs {between} points {spots} law {} 2ab {}",1129 tag(out),1130 pair_count(a, b, out, p, q, 400_000),1131 2 * a * b1132 );1133 }1134 }1135 }1136 println!("THE PAIR POINTS the pair law counts parameter pairs, and a point count needs no crossing of one curve to sit on the other");1137 for (a, b, out, lam, name) in [1138 (3, 1, false, 5.0_f64.sqrt() - 1.0, "witness"),1139 (3, 1, false, (5.0_f64.sqrt() - 1.0) * 0.98, "under"),1140 (3, 1, false, (5.0_f64.sqrt() - 1.0) * 1.02, "over"),1141 (5, 1, false, 2.6, "generic"),1142 ] {1143 let (p, q) = ((lam, 0.0), (0.0, 0.0));1144 let (curves, box_) = draw(a, b, out, &[p, q], 24001);1145 let (_, between, _, _, spots) = nodes(&curves, box_, 260, 0, 0.0);1146 println!(1147 " {} {a}/{b} reach {lam:.9} {name} law pairs {} trace pairs {between} trace points {spots} 2ab {}",1148 tag(out),1149 pair_count(a, b, out, p, q, 400_000),1150 2 * a * b1151 );1152 }1153 println!(1154 "THE PAIR THRESHOLD two seats at one reach, half angle nu apart, the reach where 2ab dies"1155 );1156 for (a, b, out) in [(5, 1, false), (7, 2, false), (8, 3, false), (5, 2, true)] {1157 let mut row = Vec::new();1158 for nu in [1.0_f64, 0.5, 0.2, 0.05, 0.01, 0.001, 0.0001] {1159 let alive = |lam: f64| {1160 let p = (lam * nu.cos(), lam * nu.sin());1161 let q = (lam * nu.cos(), -lam * nu.sin());1162 pair_count(a, b, out, p, q, 20_000) as i64 == 2 * a * b1163 };1164 let (mut lo, mut hi) = (1.0, f64::NAN);1165 for k in 0..=600 {1166 let lam = 1.0 + 1e-6 * 1e7_f64.powf(k as f64 / 600.0);1167 if alive(lam) {1168 lo = lam;1169 } else {1170 hi = lam;1171 break;1172 }1173 }1174 for _ in 0..40 {1175 let mid = 0.5 * (lo + hi);1176 if alive(mid) {1177 lo = mid;1178 } else {1179 hi = mid;1180 }1181 }1182 row.push(format!("{nu}:{hi:.6}"));1183 }1184 println!(1185 " {} {a}/{b} 2ab {} first departure by nu {}",1186 tag(out),1187 2 * a * b,1188 row.join(" ")1189 );1190 }1191}11921193fn main() {1194 axes();1195 staircase();1196 signs(220);1197 census(24);1198 mirror(26);1199 contact();1200 centre();1201 pairs();1202 hunt();1203}