main.rs
11.4 kB · rust · 428 lines
1use std::collections::HashSet;2use std::env;34// THE OBJECTS56#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)]7struct Circle {8 k: i64,9 x: i64,10 y: i64,11}1213type Quad = [Circle; 4];1415fn form(u: [i64; 4], v: [i64; 4]) -> i128 {16 let su: i128 = u.iter().map(|&a| a as i128).sum();17 let sv: i128 = v.iter().map(|&a| a as i128).sum();18 let dot: i128 = (0..4).map(|i| u[i] as i128 * v[i] as i128).sum();19 su * sv - 2 * dot20}2122fn curvatures(q: &Quad) -> [i64; 4] {23 [q[0].k, q[1].k, q[2].k, q[3].k]24}2526fn abscissae(q: &Quad) -> [i64; 4] {27 [q[0].x, q[1].x, q[2].x, q[3].x]28}2930fn ordinates(q: &Quad) -> [i64; 4] {31 [q[0].y, q[1].y, q[2].y, q[3].y]32}3334fn descartes(q: &Quad) -> bool {35 let (k, x, y) = (curvatures(q), abscissae(q), ordinates(q));36 form(k, k) == 0 && form(k, x) == 0 && form(k, y) == 037}3839fn frame(q: &Quad) -> bool {40 let (x, y) = (abscissae(q), ordinates(q));41 form(x, x) == -4 && form(y, y) == -4 && form(x, y) == 042}4344fn reflect(q: &Quad, i: usize) -> Circle {45 let (mut k, mut x, mut y) = (0i64, 0i64, 0i64);46 for j in 0..4 {47 if j != i {48 k += q[j].k;49 x += q[j].x;50 y += q[j].y;51 }52 }53 Circle {54 k: 2 * k - q[i].k,55 x: 2 * x - q[i].x,56 y: 2 * y - q[i].y,57 }58}5960fn swap(q: &Quad, i: usize) -> Quad {61 let mut out = *q;62 out[i] = reflect(q, i);63 out64}6566fn strip_root() -> Quad {67 [68 Circle { k: 0, x: 0, y: -1 },69 Circle { k: 0, x: 0, y: 1 },70 Circle { k: 2, x: 0, y: 1 },71 Circle { k: 2, x: 2, y: 1 },72 ]73}7475fn bounded_root() -> Quad {76 [77 Circle { k: -1, x: 0, y: 0 },78 Circle { k: 2, x: 1, y: 0 },79 Circle { k: 2, x: -1, y: 0 },80 Circle { k: 3, x: 0, y: 2 },81 ]82}8384// THE CENSUS8586struct Run {87 circles: u64,88 quads: u64,89 broken: u64,90 backward: u64,91 strayed: u64,92 bins: Vec<u64>,93 res: [u64; 24],94 depth: usize,95 bottom: u64,96 top: u64,97 offford: u64,98}99100fn grow(101 root: Quad,102 first: &[usize],103 steps: &[i64],104 strip: bool,105 watch: Option<&mut HashSet<Circle>>,106) -> Run {107 let cap = *steps.last().unwrap();108 let mut run = Run {109 circles: 0,110 quads: 1,111 broken: 0,112 backward: 0,113 strayed: 0,114 bins: vec![0; steps.len()],115 res: [0; 24],116 depth: 0,117 bottom: 0,118 top: 0,119 offford: 0,120 };121 let mut seen = watch;122 if !descartes(&root) || !frame(&root) {123 run.broken += 1;124 }125 let mut stack: Vec<(Quad, usize, usize)> = Vec::new();126 for &i in first {127 let c = reflect(&root, i);128 if c.k <= cap {129 stack.push((swap(&root, i), i, 1));130 }131 }132 while let Some((q, last, depth)) = stack.pop() {133 run.quads += 1;134 run.depth = run.depth.max(depth);135 if !descartes(&q) || !frame(&q) {136 run.broken += 1;137 }138 let c = q[last];139 run.circles += 1;140 run.res[c.k.rem_euclid(24) as usize] += 1;141 let at = steps.iter().position(|&s| c.k <= s).unwrap();142 run.bins[at] += 1;143 if strip {144 if c.y == 1 {145 run.bottom += 1;146 if !is_ford(c.k, c.x) {147 run.offford += 1;148 }149 }150 if c.y == c.k - 1 {151 run.top += 1;152 }153 if c.x <= 0 || c.x >= c.k {154 run.strayed += 1;155 }156 }157 if let Some(set) = seen.as_mut() {158 if !set.insert(c) {159 run.backward += 1;160 }161 }162 for j in 0..4 {163 if j == last {164 continue;165 }166 let n = reflect(&q, j);167 if n.k > cap {168 continue;169 }170 if n.k <= q[j].k {171 run.backward += 1;172 continue;173 }174 stack.push((swap(&q, j), j, depth + 1));175 }176 }177 run178}179180fn totals(bins: &[u64]) -> Vec<u64> {181 let mut out = Vec::with_capacity(bins.len());182 let mut acc = 0u64;183 for &b in bins {184 acc += b;185 out.push(acc);186 }187 out188}189190fn is_ford(k: i64, x: i64) -> bool {191 if k % 2 != 0 {192 return false;193 }194 let s = k / 2;195 let b = (s as f64).sqrt().round() as i64;196 if b * b != s || b <= 0 {197 return false;198 }199 if x % (2 * b) != 0 {200 return false;201 }202 let a = x / (2 * b);203 a >= 0 && a <= b && gcd(a, b) == 1204}205206fn gcd(a: i64, b: i64) -> i64 {207 if b == 0 {208 a.abs()209 } else {210 gcd(b, a % b)211 }212}213214// THE FORD IDENTIFICATION215216struct Ford {217 nodes: u64,218 broken: u64,219 missed: u64,220 untangent: u64,221 bright: u128,222 deepest: i64,223}224225fn ford(top: i64) -> Ford {226 let mut f = Ford {227 nodes: 0,228 broken: 0,229 missed: 0,230 untangent: 0,231 bright: 0,232 deepest: 0,233 };234 let root = strip_root();235 let base = [root[0], root[2], root[3], root[1]];236 let mut stack = vec![(base, 0i64, 1i64, 1i64, 1i64)];237 while let Some((q, a, b, c, d)) = stack.pop() {238 let (p, r) = (a + c, b + d);239 if r > top {240 continue;241 }242 if (a * d - b * c) * (a * d - b * c) != 1 {243 f.untangent += 1;244 }245 let m = reflect(&q, 3);246 let child = [q[0], q[1], m, q[2]];247 if !descartes(&child) || !frame(&child) {248 f.broken += 1;249 }250 if m.k != 2 * r * r || m.x != 2 * p * r || m.y != 1 {251 f.missed += 1;252 }253 let lhs = (0 + 2 * b * b + 2 * d * d + 2 * r * r) as i128;254 let rhs: i128 = 2255 * ((2 * b * b) as i128 * (2 * b * b) as i128256 + (2 * d * d) as i128 * (2 * d * d) as i128257 + (2 * r * r) as i128 * (2 * r * r) as i128);258 if lhs * lhs != rhs {259 f.broken += 1;260 }261 if 2 * (2 * b * b + 2 * d * d) - 2 * (b - d) * (b - d) != 2 * r * r {262 f.broken += 1;263 }264 f.nodes += 1;265 f.deepest = f.deepest.max(r);266 f.bright += (top / r) as u128;267 stack.push(([q[0], q[1], m, q[2]], a, b, p, r));268 stack.push(([q[0], m, q[2], q[1]], p, r, c, d));269 }270 f.bright += (top / 1) as u128;271 f272}273274fn totients(top: usize) -> Vec<u64> {275 let mut phi: Vec<u64> = (0..=top as u64).collect();276 for i in 2..=top {277 if phi[i] == i as u64 {278 let mut j = i;279 while j <= top {280 phi[j] -= phi[j] / i as u64;281 j += i;282 }283 }284 }285 phi286}287288// THE VERBS289290const DELTA: f64 = 1.3056867280498771846459862068510;291292fn show(name: &str, run: &Run, steps: &[i64]) {293 let tot = totals(&run.bins);294 println!(295 "{name}: {} circles, {} quadruples, depth {}, broken {}, backward {}, strayed {}",296 run.circles, run.quads, run.depth, run.broken, run.backward, run.strayed297 );298 println!(" T N(T) log N / log T ratio");299 let mut last: Option<(i64, u64)> = None;300 for (i, &t) in steps.iter().enumerate() {301 let n = tot[i];302 if n == 0 {303 continue;304 }305 let e = (n as f64).ln() / (t as f64).ln();306 let r = match last {307 Some((pt, pn)) if pn > 0 => {308 format!(309 "{:.4}",310 ((n as f64 / pn as f64).ln()) / ((t as f64 / pt as f64).ln())311 )312 }313 _ => "-".to_string(),314 };315 println!(" {t} {n} {e:.4} {r}");316 last = Some((t, n));317 }318 let live: Vec<String> = (0..24)319 .filter(|&r| run.res[r] > 0)320 .map(|r| format!("{r}:{}", run.res[r]))321 .collect();322 println!(" mod 24 {}", live.join(" "));323}324325fn main() {326 let verb = env::args().nth(1).unwrap_or_else(|| "all".to_string());327 let all = verb == "all";328329 if all || verb == "ford" {330 for top in [50i64, 200, 1000, 4000] {331 let f = ford(top);332 let phi = totients(top as usize);333 let want: u64 = phi[1..=top as usize].iter().sum::<u64>() - 1;334 let bright = (top as u128) * (top as u128 + 1) / 2;335 println!(336 "ford b<={top}: nodes {} want {} broken {} missed {} untangent {} bright {} want {}",337 f.nodes, want, f.broken, f.missed, f.untangent, f.bright, bright338 );339 println!(340 " curvature at b={top} is {} deepest {}",341 2 * top * top,342 f.deepest343 );344 }345 }346347 if all || verb == "strip" {348 let steps = [349 2i64, 8, 32, 128, 512, 2048, 8192, 32768, 131072, 524288, 2097152,350 ];351 let mut set = HashSet::new();352 let small = grow(strip_root(), &[0, 1], &steps[..6], true, Some(&mut set));353 println!(354 "strip control T<=2048: {} circles, {} distinct, bottom {} top {}",355 small.circles,356 set.len(),357 small.bottom,358 small.top359 );360 let run = grow(strip_root(), &[0, 1], &steps, true, None);361 show("strip octaves", &run, &steps);362 let decades = [10i64, 100, 1000, 10000, 100000, 1000000];363 let dec = grow(strip_root(), &[0, 1], &decades, true, None);364 show("strip decades", &dec, &decades);365 println!(366 " bottom-tangent {} top-tangent {} off-Ford {}",367 run.bottom, run.top, run.offford368 );369 for q in [32i64, 181, 1024] {370 let f = grow(strip_root(), &[0, 1], &[2 * q * q], true, None);371 let phi = totients(q as usize);372 let want: u64 = phi[1..=q as usize].iter().sum::<u64>() - 1;373 println!(" line-tangent below 2*{q}^2: {} want {}", f.bottom, want);374 }375 }376377 if all || verb == "census" {378 let steps = [10i64, 100, 1000, 10000, 100000, 1000000, 10000000];379 let mut set = HashSet::new();380 let small = grow(381 bounded_root(),382 &[0, 1, 2, 3],383 &steps[..4],384 false,385 Some(&mut set),386 );387 println!(388 "bounded control T<=10000: {} circles, {} distinct",389 small.circles,390 set.len()391 );392 let run = grow(bounded_root(), &[0, 1, 2, 3], &steps, false, None);393 show("bounded (-1,2,2,3)", &run, &steps);394 }395396 if all || verb == "design" {397 let mut best: Vec<(f64, i64, f64)> = Vec::new();398 for q in 2i64..=100 {399 let v = (q as f64).powf(DELTA);400 let n = v.round();401 best.push(((v - n).abs(), q, v));402 }403 best.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());404 println!("design: q^delta against the integers, delta = {DELTA}");405 for &(g, q, v) in best.iter().take(6) {406 let d = v.round().ln() / (q as f64).ln();407 println!(408 " q {q} q^delta {v:.9} gap {g:.9} log N / log q {d:.9} off {:.9}",409 (d - DELTA).abs()410 );411 }412 let w = best.last().unwrap();413 println!(" worst gap {:.9} at q {}", w.0, w.1);414 let mut near: Vec<(f64, i64, f64)> = best415 .iter()416 .map(|&(_, q, v)| {417 let d = v.round().ln() / (q as f64).ln();418 ((d - DELTA).abs(), q, d)419 })420 .collect();421 near.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());422 let n = near[0];423 println!(424 " nearest design dimension {:.9} at q {} off {:.9}",425 n.2, n.1, n.0426 );427 }428}