main.rs
55.7 kB · rust · 1780 lines
1use mrlyrs::num::gauss::Ring;2use mrlyrs::num::radix::{flowsnake, koch, terdragon, twindragon, Base, Radix};3use std::collections::{HashMap, HashSet};4use std::time::Instant;56type Z = (i64, i64);78const BASES: [(Ring, Z); 7] = [9 (Ring::Gaussian, (2, 0)),10 (Ring::Gaussian, (1, 1)),11 (Ring::Gaussian, (2, 1)),12 (Ring::Eisenstein, (2, 0)),13 (Ring::Eisenstein, (2, 1)),14 (Ring::Eisenstein, (3, 0)),15 (Ring::Eisenstein, (3, 1)),16];1718// RING1920fn add(z: Z, w: Z) -> Z {21 (z.0 + w.0, z.1 + w.1)22}2324fn sub(z: Z, w: Z) -> Z {25 (z.0 - w.0, z.1 - w.1)26}2728#[derive(Clone)]29struct Lab {30 ring: Ring,31 base: Base,32 b: Z,33 n: usize,34 units: Vec<Z>,35 residues: Vec<Z>,36}3738impl Lab {39 fn new(ring: Ring, value: Z) -> Lab {40 let base = Base::new(ring, value).unwrap();41 Lab {42 ring,43 base,44 b: value,45 n: ring.units(),46 units: ring.associates(1, 0),47 residues: base.residues().unwrap(),48 }49 }50 fn mul(&self, z: Z, w: Z) -> Z {51 self.ring.mul(z, w)52 }53 fn norm(&self, z: Z) -> u64 {54 self.ring.norm(z.0, z.1)55 }56 fn turn(&self, z: Z, k: usize) -> Z {57 self.mul(z, self.units[k % self.n])58 }59 fn back(&self, k: usize) -> usize {60 (self.n - k % self.n) % self.n61 }62 fn conj(&self, z: Z) -> Z {63 self.ring.conjugate(z.0, z.1)64 }65 fn mirror(&self) -> Option<usize> {66 let c = self.conj(self.b);67 (0..self.n).find(|&k| self.turn(self.b, k) == c)68 }69 fn sum(&self, walk: &[usize]) -> Z {70 walk.iter().fold((0, 0), |acc, &k| add(acc, self.units[k]))71 }72 fn power(&self, level: usize) -> Z {73 self.base.power(level)74 }75 fn modulus(&self) -> f64 {76 (self.base.norm() as f64).sqrt()77 }78}7980fn mul_wide(ring: Ring, (a, b): (i128, i128), (c, d): (i128, i128)) -> (i128, i128) {81 match ring {82 Ring::Gaussian => (a * c - b * d, a * d + b * c),83 Ring::Eisenstein => (a * c - b * d, a * d + b * c - b * d),84 }85}8687fn norm_wide(ring: Ring, (a, b): (i128, i128)) -> i128 {88 match ring {89 Ring::Gaussian => a * a + b * b,90 Ring::Eisenstein => a * a - a * b + b * b,91 }92}9394fn wide(z: Z) -> (i128, i128) {95 (z.0 as i128, z.1 as i128)96}9798fn spell(ring: Ring, z: Z) -> String {99 let (re, im, unit) = match ring {100 Ring::Gaussian => (z.0, z.1, "i"),101 Ring::Eisenstein => (z.0 - z.1, z.1, "w"),102 };103 match (re, im) {104 (a, 0) => format!("{a}"),105 (0, 1) => unit.to_string(),106 (0, -1) => format!("-{unit}"),107 (0, b) => format!("{b}{unit}"),108 (a, 1) => format!("{a}+{unit}"),109 (a, -1) => format!("{a}-{unit}"),110 (a, b) if b > 0 => format!("{a}+{b}{unit}"),111 (a, b) => format!("{a}{b}{unit}"),112 }113}114115fn ring_word(ring: Ring) -> &'static str {116 match ring {117 Ring::Gaussian => "Z[i]",118 Ring::Eisenstein => "Z[w]",119 }120}121122fn spell_list(ring: Ring, list: &[Z]) -> String {123 list.iter()124 .map(|&z| spell(ring, z))125 .collect::<Vec<_>>()126 .join(" ")127}128129fn spell_turns(turns: &[usize]) -> String {130 turns131 .iter()132 .map(|k| k.to_string())133 .collect::<Vec<_>>()134 .join("")135}136137// DESIGN138139#[derive(Clone, Debug, PartialEq, Eq, Hash)]140struct Design {141 digits: Vec<Z>,142 twists: Vec<usize>,143}144145impl Design {146 fn radix(&self, lab: &Lab) -> Radix {147 let twists = self.twists.iter().map(|&k| lab.units[k]).collect();148 Radix::new(lab.base, self.digits.clone(), twists).unwrap()149 }150 fn size(&self) -> usize {151 self.digits.len()152 }153}154155fn canonical(lab: &Lab, code: usize, twists: &[usize]) -> Design {156 let digits: Vec<Z> = (0..lab.residues.len())157 .filter(|i| (code >> i) & 1 == 1)158 .map(|i| lab.residues[i])159 .collect();160 Design {161 digits,162 twists: twists.to_vec(),163 }164}165166fn walk_design(lab: &Lab, walk: &[usize]) -> Design {167 let mut digits = Vec::with_capacity(walk.len());168 let mut at = (0, 0);169 for &k in walk {170 digits.push(at);171 at = add(at, lab.units[k]);172 }173 Design {174 digits,175 twists: walk.to_vec(),176 }177}178179fn points(lab: &Lab, design: &Design, level: usize) -> Vec<Z> {180 let mut out = vec![(0i64, 0i64)];181 for step in 1..=level {182 let shift = lab.power(step - 1);183 let mut next = Vec::with_capacity(out.len() * design.size());184 for (&d, &t) in design.digits.iter().zip(design.twists.iter()) {185 let head = lab.mul(d, shift);186 for &p in &out {187 next.push(add(head, lab.turn(p, t)));188 }189 }190 out = next;191 }192 out193}194195fn incongruent(lab: &Lab, digits: &[Z]) -> bool {196 (0..digits.len()).all(|i| (0..i).all(|j| !lab.base.congruent(digits[i], digits[j])))197}198199// STEP200201fn junction_key(lab: &Lab, design: &Design, a: usize) -> Z {202 let k = design.size();203 let (d0, dk) = (design.digits[0], design.digits[k - 1]);204 let f0 = sub(lab.b, lab.units[design.twists[0]]);205 let fk = sub(lab.b, lab.units[design.twists[k - 1]]);206 let gap = sub(design.digits[a + 1], design.digits[a]);207 let one = lab.mul(lab.mul(gap, f0), fk);208 let two = lab.mul(lab.mul(lab.units[design.twists[a + 1]], d0), fk);209 let three = lab.mul(lab.mul(lab.units[design.twists[a]], dk), f0);210 sub(add(one, two), three)211}212213fn joints(lab: &Lab, design: &Design, top: usize) -> Vec<Vec<(i128, i128)>> {214 let ring = lab.ring;215 let k = design.size();216 let (d0, dk) = (wide(design.digits[0]), wide(design.digits[k - 1]));217 let (u0, uk) = (218 wide(lab.units[design.twists[0]]),219 wide(lab.units[design.twists[k - 1]]),220 );221 let b = wide(lab.b);222 let mut out = vec![Vec::with_capacity(top); k - 1];223 let (mut first, mut last, mut scale) = ((0i128, 0i128), (0i128, 0i128), (1i128, 0i128));224 for _ in 1..=top {225 for (a, row) in out.iter_mut().enumerate() {226 let gap = wide(sub(design.digits[a + 1], design.digits[a]));227 let one = mul_wide(ring, gap, scale);228 let two = mul_wide(ring, wide(lab.units[design.twists[a + 1]]), first);229 let three = mul_wide(ring, wide(lab.units[design.twists[a]]), last);230 row.push((one.0 + two.0 - three.0, one.1 + two.1 - three.1));231 }232 let f = mul_wide(ring, u0, first);233 let l = mul_wide(ring, uk, last);234 let df = mul_wide(ring, d0, scale);235 let dl = mul_wide(ring, dk, scale);236 first = (df.0 + f.0, df.1 + f.1);237 last = (dl.0 + l.0, dl.1 + l.1);238 scale = mul_wide(ring, scale, b);239 }240 out241}242243fn steps_by_level(lab: &Lab, design: &Design, top: usize) -> Vec<bool> {244 let table = joints(lab, design, top);245 (0..top)246 .map(|l| table.iter().all(|row| norm_wide(lab.ring, row[l]) == 1))247 .collect()248}249250fn lemma_step(lab: &Lab, design: &Design) -> bool {251 let k = design.size();252 k >= 2253 && (0..k - 1).all(|a| junction_key(lab, design, a) == (0, 0))254 && steps_by_level(lab, design, lab.n).iter().all(|&s| s)255}256257fn first_break(lab: &Lab, design: &Design) -> Option<usize> {258 if lemma_step(lab, design) {259 return None;260 }261 let top = 36;262 let found = steps_by_level(lab, design, top).iter().position(|&s| !s);263 assert!(264 found.is_some(),265 "a design outside the lemma kept unit steps to level {top}"266 );267 found.map(|l| l + 1)268}269270fn brute_steps(lab: &Lab, design: &Design, level: usize) -> bool {271 let p = points(lab, design, level);272 p.windows(2).all(|w| lab.norm(sub(w[1], w[0])) == 1)273}274275// MACHINE276277struct Machine {278 points: Vec<Z>,279 turn: Vec<u32>,280 near: Vec<bool>,281 zero: u32,282 k: usize,283 n: usize,284 head: Vec<u32>,285 edges: Vec<(u8, u8, u32)>,286 gaps: Vec<Option<u32>>,287 ends: Vec<Option<u32>>,288}289290fn machine(lab: &Lab, digits: &[Z]) -> Machine {291 let n = lab.n;292 let k = digits.len();293 let reach = digits.iter().map(|&d| lab.norm(d)).max().unwrap_or(0) as f64;294 let rho = 2.0 * reach.sqrt() / (lab.modulus() - 1.0);295 let cap = (rho * rho + 1e-9).max(1.0);296 let side = (2.0 * cap).sqrt().ceil() as i64 + 1;297 let mut points = Vec::new();298 for a in -side..=side {299 for b in -side..=side {300 let norm = lab.norm((a, b));301 if norm <= 1 || (norm as f64) < cap {302 points.push((a, b));303 }304 }305 }306 let width = 2 * side + 1;307 let mut grid = vec![u32::MAX; (width * width) as usize];308 for (i, &(a, b)) in points.iter().enumerate() {309 grid[((a + side) * width + b + side) as usize] = i as u32;310 }311 let find = |z: Z| -> Option<u32> {312 if z.0.abs() > side || z.1.abs() > side {313 return None;314 }315 let i = grid[((z.0 + side) * width + z.1 + side) as usize];316 (i != u32::MAX).then_some(i)317 };318 let mut turn = Vec::with_capacity(points.len() * n);319 for &z in &points {320 for j in 0..n {321 turn.push(find(lab.turn(z, j)).unwrap());322 }323 }324 let near = points.iter().map(|&z| lab.norm(z) <= 1).collect();325 let zero = find((0, 0)).unwrap();326 let mut head = vec![0u32];327 let mut edges = Vec::new();328 for &delta in &points {329 let scaled = lab.mul(lab.b, delta);330 for r in 0..n {331 for (e, &de) in digits.iter().enumerate() {332 for (f, &df) in digits.iter().enumerate() {333 let x = sub(add(scaled, de), lab.turn(df, r));334 if let Some(i) = find(x) {335 edges.push((e as u8, f as u8, i));336 }337 }338 }339 head.push(edges.len() as u32);340 }341 }342 let mut gaps = vec![None; k * k];343 for d in 0..k {344 for f in 0..k {345 if d != f {346 gaps[d * k + f] = find(sub(digits[d], digits[f]));347 }348 }349 }350 let ends = digits.iter().map(|&d| find(sub(lab.b, d))).collect();351 Machine {352 points,353 turn,354 near,355 zero,356 k,357 n,358 head,359 edges,360 gaps,361 ends,362 }363}364365impl Machine {366 fn states(&self) -> usize {367 self.points.len() * self.n368 }369 fn back(&self, k: usize) -> usize {370 (self.n - k % self.n) % self.n371 }372 fn hop(&self, t: &[usize], s: usize, edge: (u8, u8, u32)) -> usize {373 let (e, f, x) = (edge.0 as usize, edge.1 as usize, edge.2 as usize);374 let r = s % self.n;375 let p = self.turn[x * self.n + self.back(t[e])] as usize;376 p * self.n + (r + t[f] + self.n - t[e]) % self.n377 }378 fn start(&self, t: &[usize], d: usize, f: usize) -> Option<usize> {379 self.gaps[d * self.k + f].map(|g| {380 let p = self.turn[g as usize * self.n + self.back(t[d])] as usize;381 p * self.n + (t[f] + self.n - t[d]) % self.n382 })383 }384 fn end_start(&self, t: &[usize], f: usize) -> Option<usize> {385 self.ends[f].map(|g| g as usize * self.n + t[f] % self.n)386 }387 fn out(&self, s: usize) -> &[(u8, u8, u32)] {388 &self.edges[self.head[s] as usize..self.head[s + 1] as usize]389 }390 fn search(&self, t: &[usize], seeds: &[usize], maps: bool, tail: bool) -> Option<usize> {391 let mut dist = vec![u32::MAX; self.states()];392 let mut queue = std::collections::VecDeque::new();393 for &s in seeds {394 if dist[s] == u32::MAX {395 dist[s] = 0;396 queue.push_back(s);397 }398 }399 while let Some(s) = queue.pop_front() {400 let p = s / self.n;401 if p as u32 == self.zero && (!maps || s % self.n == 0) {402 return Some(dist[s] as usize + 1);403 }404 for &edge in self.out(s) {405 if tail && edge.0 != 0 {406 continue;407 }408 let next = self.hop(t, s, edge);409 if dist[next] == u32::MAX {410 dist[next] = dist[s] + 1;411 queue.push_back(next);412 }413 }414 }415 None416 }417 fn glue(&self, t: &[usize], maps: bool) -> Option<usize> {418 let mut seeds = Vec::new();419 for d in 0..self.k {420 for f in d + 1..self.k {421 if let Some(s) = self.start(t, d, f) {422 seeds.push(s);423 }424 }425 }426 self.search(t, &seeds, maps, false)427 }428 fn ending(&self, t: &[usize]) -> Option<usize> {429 let seeds: Vec<usize> = (0..self.k).filter_map(|f| self.end_start(t, f)).collect();430 self.search(t, &seeds, false, true)431 }432}433434fn linked(k: usize, edges: impl Iterator<Item = (usize, usize)>) -> bool {435 let mut root: Vec<usize> = (0..k).collect();436 fn find(root: &mut [usize], mut x: usize) -> usize {437 while root[x] != x {438 root[x] = root[root[x]];439 x = root[x];440 }441 x442 }443 let mut parts = k;444 for (a, b) in edges {445 let (x, y) = (find(&mut root, a), find(&mut root, b));446 if x != y {447 root[x] = y;448 parts -= 1;449 }450 }451 parts <= 1452}453454impl Machine {455 fn piece(&self, t: &[usize]) -> Option<usize> {456 let k = self.k;457 if k <= 1 {458 return None;459 }460 let mut pairs = Vec::new();461 for d in 0..k {462 for f in d + 1..k {463 if let Some(s) = self.start(t, d, f) {464 pairs.push((d, f, s));465 }466 }467 }468 let mut local = vec![u32::MAX; self.states()];469 let mut order: Vec<usize> = Vec::new();470 for &(_, _, s) in &pairs {471 if local[s] == u32::MAX {472 local[s] = order.len() as u32;473 order.push(s);474 }475 }476 let mut head = vec![0u32];477 let mut succ: Vec<u32> = Vec::new();478 let mut i = 0;479 while i < order.len() {480 let s = order[i];481 for &edge in self.out(s) {482 let next = self.hop(t, s, edge);483 if local[next] == u32::MAX {484 local[next] = order.len() as u32;485 order.push(next);486 }487 succ.push(local[next]);488 }489 head.push(succ.len() as u32);490 i += 1;491 }492 let size = order.len();493 let words = size.div_ceil(64);494 let bit = |w: &[u64], j: usize| (w[j / 64] >> (j % 64)) & 1 == 1;495 let mut w = vec![0u64; words];496 for (j, &s) in order.iter().enumerate() {497 if self.near[s / self.n] {498 w[j / 64] |= 1 << (j % 64);499 }500 }501 let mut history: Vec<Vec<u64>> = Vec::new();502 for level in 1.. {503 let edges = pairs504 .iter()505 .filter(|&&(_, _, s)| bit(&w, local[s] as usize))506 .map(|&(d, f, _)| (d, f));507 if !linked(k, edges) {508 return Some(level);509 }510 history.push(w.clone());511 let mut next = vec![0u64; words];512 for j in 0..size {513 let list = &succ[head[j] as usize..head[j + 1] as usize];514 if list.iter().any(|&x| bit(&w, x as usize)) {515 next[j / 64] |= 1 << (j % 64);516 }517 }518 if history.contains(&next) {519 return None;520 }521 assert!(history.len() < 4096, "the adjacency sets have not cycled");522 w = next;523 }524 None525 }526}527528fn brute_piece(lab: &Lab, design: &Design, level: usize) -> bool {529 let set: HashSet<Z> = points(lab, design, level).into_iter().collect();530 let start = match set.iter().next() {531 Some(&z) => z,532 None => return true,533 };534 let mut seen = HashSet::from([start]);535 let mut stack = vec![start];536 while let Some(z) = stack.pop() {537 for &u in &lab.units {538 let w = add(z, u);539 if set.contains(&w) && seen.insert(w) {540 stack.push(w);541 }542 }543 }544 seen.len() == set.len()545}546547fn brute_glue(lab: &Lab, design: &Design, level: usize) -> (bool, bool) {548 let mut place = HashSet::new();549 let mut map = HashSet::new();550 let mut word_twist = vec![0usize];551 for _ in 0..level {552 let mut next = Vec::with_capacity(word_twist.len() * design.size());553 for &t in &design.twists {554 for &w in &word_twist {555 next.push((t + w) % lab.n);556 }557 }558 word_twist = next;559 }560 let p = points(lab, design, level);561 for (z, &u) in p.iter().zip(word_twist.iter()) {562 place.insert(*z);563 map.insert((*z, u));564 }565 (place.len() < p.len(), map.len() < p.len())566}567568// GROUP569570#[derive(Clone)]571struct Symmetry {572 perm: Vec<usize>,573 twist: Vec<usize>,574 mirror: bool,575 unit: usize,576}577578fn residue_index(lab: &Lab) -> HashMap<Z, usize> {579 lab.residues580 .iter()581 .enumerate()582 .map(|(i, &z)| (z, i))583 .collect()584}585586fn stabiliser(lab: &Lab) -> Vec<Symmetry> {587 let index = residue_index(lab);588 let mut out = Vec::new();589 for unit in 0..lab.n {590 let perm: Option<Vec<usize>> = lab591 .residues592 .iter()593 .map(|&z| index.get(&lab.turn(z, unit)).copied())594 .collect();595 if let Some(perm) = perm {596 out.push(Symmetry {597 perm,598 twist: (0..lab.n).collect(),599 mirror: false,600 unit,601 });602 }603 }604 if let Some(eps) = lab.mirror() {605 for unit in 0..lab.n {606 let perm: Option<Vec<usize>> = lab607 .residues608 .iter()609 .map(|&z| index.get(&lab.turn(lab.conj(z), unit)).copied())610 .collect();611 if let Some(perm) = perm {612 let twist = (0..lab.n)613 .map(|k| (lab.back(k) + lab.back(eps)) % lab.n)614 .collect();615 out.push(Symmetry {616 perm,617 twist,618 mirror: true,619 unit,620 });621 }622 }623 }624 out625}626627fn act(sym: &Symmetry, code: usize, twists: &[usize]) -> (usize, Vec<usize>) {628 let q = sym.perm.len();629 let mut slot = vec![usize::MAX; q];630 let mut j = 0;631 let mut image = 0usize;632 for c in 0..q {633 if (code >> c) & 1 == 1 {634 slot[sym.perm[c]] = sym.twist[twists[j]];635 image |= 1 << sym.perm[c];636 j += 1;637 }638 }639 (640 image,641 slot.into_iter().filter(|&t| t != usize::MAX).collect(),642 )643}644645fn pack(n: usize, twists: &[usize]) -> u64 {646 twists647 .iter()648 .rev()649 .fold(0u64, |acc, &t| acc * n as u64 + t as u64)650}651652fn unpack(n: usize, k: usize, mut index: u64) -> Vec<usize> {653 let mut out = Vec::with_capacity(k);654 for _ in 0..k {655 out.push((index % n as u64) as usize);656 index /= n as u64;657 }658 out659}660661fn key(n: usize, code: usize, twists: &[usize]) -> u64 {662 ((code as u64) << 40) | pack(n, twists)663}664665fn fixed(lab: &Lab, sym: &Symmetry) -> u64 {666 let q = lab.residues.len();667 let mut total = 0u64;668 for code in 0..1usize << q {669 if code.count_ones() < 2 || act(sym, code, &vec![0; code.count_ones() as usize]).0 != code {670 continue;671 }672 let mut seen = 0usize;673 let mut count = 1u64;674 for c in 0..q {675 if (code >> c) & 1 == 0 || (seen >> c) & 1 == 1 {676 continue;677 }678 let mut len = 0;679 let mut at = c;680 while (seen >> at) & 1 == 0 {681 seen |= 1 << at;682 at = sym.perm[at];683 len += 1;684 }685 let stay = (0..lab.n)686 .filter(|&t| {687 let mut x = t;688 for _ in 0..len {689 x = sym.twist[x];690 }691 x == t692 })693 .count();694 count *= stay as u64;695 }696 total += count;697 }698 total699}700701fn burnside(lab: &Lab, group: &[Symmetry]) -> u64 {702 let total: u64 = group.iter().map(|sym| fixed(lab, sym)).sum();703 assert_eq!(total % group.len() as u64, 0, "Burnside is not an integer");704 total / group.len() as u64705}706707fn same_up_to_turn(lab: &Lab, left: &[Z], right: &[Z]) -> bool {708 let target: HashSet<Z> = right.iter().copied().collect();709 (0..lab.n).any(|v| left.iter().all(|&z| target.contains(&lab.turn(z, v))))710 && left.len() == right.len()711}712713// CENSUS714715#[derive(Clone, Default)]716struct Cell {717 pairs: u64,718 orbits: u64,719 step: u64,720 curve: u64,721 piece: u64,722 piece_orbits: u64,723 plane: u64,724 plane_orbits: u64,725 distinct: u64,726 distinct_orbits: u64,727}728729impl Cell {730 fn merge(&mut self, o: &Cell) {731 self.pairs += o.pairs;732 self.orbits += o.orbits;733 self.step += o.step;734 self.curve += o.curve;735 self.piece += o.piece;736 self.piece_orbits += o.piece_orbits;737 self.plane += o.plane;738 self.plane_orbits += o.plane_orbits;739 self.distinct += o.distinct;740 self.distinct_orbits += o.distinct_orbits;741 }742}743744type Mark = (usize, usize, Vec<usize>);745746#[derive(Clone, Default)]747struct Survey {748 cells: Vec<Cell>,749 steps: Vec<(usize, Vec<usize>, Option<usize>)>,750 held: Mark,751 split: Mark,752 overlap: Mark,753}754755fn raise(mark: &mut Mark, level: usize, code: usize, twists: &[usize]) {756 if level > mark.0 || (level == mark.0 && (code, twists) < (mark.1, &mark.2[..])) {757 *mark = (level, code, twists.to_vec());758 }759}760761fn merge_mark(mark: &mut Mark, other: &Mark) {762 if other.0 > 0 {763 raise(mark, other.0, other.1, &other.2);764 }765}766767fn unit_path(lab: &Lab, digits: &[Z]) -> bool {768 digits.windows(2).all(|w| lab.norm(sub(w[1], w[0])) == 1)769}770771fn unit_linked(lab: &Lab, digits: &[Z]) -> bool {772 let k = digits.len();773 let mut edges = Vec::new();774 for a in 0..k {775 for b in a + 1..k {776 if lab.norm(sub(digits[a], digits[b])) == 1 {777 edges.push((a, b));778 }779 }780 }781 linked(k, edges.into_iter())782}783784fn survey_chunk(lab: &Lab, group: &[Symmetry], code: usize, lo: u64, hi: u64) -> Survey {785 let q = lab.residues.len();786 let n = lab.n;787 let mut design = canonical(lab, code, &[]);788 let k = design.size();789 let m = machine(lab, &design.digits);790 let path = unit_path(lab, &design.digits);791 let joined = unit_linked(lab, &design.digits);792 let mut out = Survey {793 cells: vec![Cell::default(); q + 1],794 ..Survey::default()795 };796 let cell = &mut out.cells[k];797 for index in lo..hi {798 let t = unpack(n, k, index);799 cell.pairs += 1;800 design.twists = t.clone();801 if path {802 if lemma_step(lab, &design) {803 assert_eq!(804 lab.sum(&t),805 lab.b,806 "a unit-step pair whose twists do not sum to the base"807 );808 cell.step += 1;809 let glue = m.glue(&t, false);810 if glue.is_none() {811 cell.curve += 1;812 }813 out.steps.push((code, t.clone(), glue));814 } else {815 let level = first_break(lab, &design).unwrap();816 raise(&mut out.held, level - 1, code, &t);817 }818 }819 let own = key(n, code, &t);820 let mut images = Vec::with_capacity(group.len());821 let mut rep = true;822 for sym in group {823 let (c, u) = act(sym, code, &t);824 let image = key(n, c, &u);825 if image < own {826 rep = false;827 break;828 }829 images.push(image);830 }831 if !rep {832 continue;833 }834 images.sort_unstable();835 images.dedup();836 let size = images.len() as u64;837 cell.orbits += 1;838 if joined {839 match m.piece(&t) {840 None => {841 cell.piece += size;842 cell.piece_orbits += 1;843 }844 Some(level) => raise(&mut out.split, level, code, &t),845 }846 }847 if k == q {848 match m.glue(&t, true) {849 None => {850 cell.plane += size;851 cell.plane_orbits += 1;852 }853 Some(level) => raise(&mut out.overlap, level, code, &t),854 }855 if m.glue(&t, false).is_none() {856 cell.distinct += size;857 cell.distinct_orbits += 1;858 }859 }860 }861 out862}863864fn survey(lab: &Lab, group: &[Symmetry]) -> Survey {865 let q = lab.residues.len();866 let n = lab.n as u64;867 let chunk = 1u64 << 15;868 let mut work = Vec::new();869 for code in 0..1usize << q {870 let k = code.count_ones();871 if k < 2 {872 continue;873 }874 let total = n.pow(k);875 let mut lo = 0;876 while lo < total {877 work.push((code, lo, (lo + chunk).min(total)));878 lo += chunk;879 }880 }881 let next = std::sync::atomic::AtomicUsize::new(0);882 let threads = std::thread::available_parallelism()883 .map(|x| x.get())884 .unwrap_or(4);885 let parts: Vec<Survey> = std::thread::scope(|scope| {886 let handles: Vec<_> = (0..threads)887 .map(|_| {888 scope.spawn(|| {889 let mut acc = Survey {890 cells: vec![Cell::default(); q + 1],891 ..Survey::default()892 };893 loop {894 let i = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);895 if i >= work.len() {896 break;897 }898 let (code, lo, hi) = work[i];899 let part = survey_chunk(lab, group, code, lo, hi);900 for (a, b) in acc.cells.iter_mut().zip(part.cells.iter()) {901 a.merge(b);902 }903 acc.steps.extend(part.steps);904 merge_mark(&mut acc.held, &part.held);905 merge_mark(&mut acc.split, &part.split);906 merge_mark(&mut acc.overlap, &part.overlap);907 }908 acc909 })910 })911 .collect();912 handles.into_iter().map(|h| h.join().unwrap()).collect()913 });914 let mut out = Survey {915 cells: vec![Cell::default(); q + 1],916 ..Survey::default()917 };918 for part in parts {919 for (a, b) in out.cells.iter_mut().zip(part.cells.iter()) {920 a.merge(b);921 }922 out.steps.extend(part.steps);923 merge_mark(&mut out.held, &part.held);924 merge_mark(&mut out.split, &part.split);925 merge_mark(&mut out.overlap, &part.overlap);926 }927 out.steps.sort();928 out929}930931fn code_action(lab: &Lab) -> Vec<Symmetry> {932 let perms = lab.base.group().unwrap();933 let eps = lab.mirror();934 let mut out = Vec::new();935 let mut i = 0;936 for unit in 0..lab.n {937 for flip in [false, true] {938 if flip && eps.is_none() {939 continue;940 }941 let twist = match (flip, eps) {942 (true, Some(e)) => (0..lab.n)943 .map(|k| (lab.back(k) + lab.back(e)) % lab.n)944 .collect(),945 _ => (0..lab.n).collect(),946 };947 out.push(Symmetry {948 perm: perms[i].clone(),949 twist,950 mirror: flip,951 unit,952 });953 i += 1;954 }955 }956 out957}958959fn spell_symmetry(lab: &Lab, sym: &Symmetry) -> String {960 let v = spell(lab.ring, lab.units[sym.unit]);961 if sym.mirror {962 format!("x -> {v} conj(x)")963 } else {964 format!("x -> {v} x")965 }966}967968fn breach(969 lab: &Lab,970 action: &[Symmetry],971) -> Option<(usize, Vec<usize>, usize, Vec<usize>, String)> {972 let q = lab.residues.len();973 for code in 0..1usize << q {974 let k = code.count_ones() as usize;975 if k < 2 {976 continue;977 }978 let m = machine(lab, &canonical(lab, code, &[]).digits);979 for index in 0..(lab.n as u64).pow(k as u32) {980 let t = unpack(lab.n, k, index);981 let (piece, glue) = (m.piece(&t), m.glue(&t, false));982 for sym in action {983 let (c, u) = act(sym, code, &t);984 let other = machine(lab, &canonical(lab, c, &[]).digits);985 let (p2, g2) = (other.piece(&u), other.glue(&u, false));986 let by = spell_symmetry(lab, sym);987 if piece.is_none() != p2.is_none() {988 return Some((989 code,990 t,991 c,992 u,993 format!(994 "by {by}, one piece {} against {}",995 piece.is_none(),996 p2.is_none()997 ),998 ));999 }1000 if glue.is_none() != g2.is_none() {1001 return Some((1002 code,1003 t,1004 c,1005 u,1006 format!(1007 "by {by}, no two words on one point {} against {}",1008 glue.is_none(),1009 g2.is_none()1010 ),1011 ));1012 }1013 }1014 }1015 }1016 None1017}10181019fn run_census() {1020 println!("CENSUS canonical pairs (code, twist), |F| >= 2, digits in canonical order, every property decided at every level");1021 println!(" step: unit steps between consecutive words; curve: step and no two words on one point; piece: connected under unit adjacency; plane: |F| = N(b) and no two words with one place map; distinct: |F| = N(b) and no two words on one point");1022 for (ring, value) in BASES {1023 let clock = Instant::now();1024 let lab = Lab::new(ring, value);1025 let q = lab.residues.len();1026 let n = lab.n as u64;1027 let group = stabiliser(&lab);1028 let action = code_action(&lab);1029 let s = survey(&lab, &group);1030 let mut total = Cell::default();1031 for cell in &s.cells {1032 total.merge(cell);1033 }1034 let expect = (n + 1).pow(q as u32) - 1 - q as u64 * n;1035 assert_eq!(1036 total.pairs, expect,1037 "the pair count is not sum_(k>=2) binom(q,k) n^k"1038 );1039 assert_eq!(1040 total.orbits,1041 burnside(&lab, &group),1042 "orbit walk and Burnside disagree"1043 );1044 println!(1045 " {} base {} N(b) {q} residues {}",1046 ring_word(ring),1047 spell(ring, value),1048 spell_list(ring, &lab.residues)1049 );1050 println!(1051 " stabiliser of the residue system, order {}: {}",1052 group.len(),1053 group1054 .iter()1055 .map(|g| spell_symmetry(&lab, g))1056 .collect::<Vec<_>>()1057 .join(", ")1058 );1059 println!(1060 " cut: {} pairs, {} orbits under the stabiliser; the code action of order {} would give {} orbits",1061 total.pairs,1062 total.orbits,1063 action.len(),1064 burnside(&lab, &action)1065 );1066 println!(" |F| pairs orbits step curve piece piece-orbits plane plane-orbits distinct distinct-orbits");1067 for (k, c) in s.cells.iter().enumerate() {1068 if c.pairs == 0 {1069 continue;1070 }1071 println!(1072 " {k} {} {} {} {} {} {} {} {} {} {}",1073 c.pairs,1074 c.orbits,1075 c.step,1076 c.curve,1077 c.piece,1078 c.piece_orbits,1079 c.plane,1080 c.plane_orbits,1081 c.distinct,1082 c.distinct_orbits1083 );1084 }1085 println!(1086 " all {} {} {} {} {} {} {} {} {} {}",1087 total.pairs,1088 total.orbits,1089 total.step,1090 total.curve,1091 total.piece,1092 total.piece_orbits,1093 total.plane,1094 total.plane_orbits,1095 total.distinct,1096 total.distinct_orbits1097 );1098 let show = |m: &Mark| format!("level {} at code {} twists {}", m.0, m.1, spell_turns(&m.2));1099 println!(1100 " deepest unit steps outside the lemma: {}",1101 show(&s.held)1102 );1103 println!(" deepest first split: {}", show(&s.split));1104 println!(" deepest first shared place map: {}", show(&s.overlap));1105 for (code, t, glue) in &s.steps {1106 let d = canonical(&lab, *code, t);1107 println!(1108 " step pair: code {code} digits {} twists {} first glue {}",1109 spell_list(ring, &d.digits),1110 spell_turns(t),1111 glue.map_or("never".to_string(), |g| format!("level {g}"))1112 );1113 }1114 if q <= 5 {1115 match breach(&lab, &action) {1116 Some((c, t, c2, u, why)) => println!(1117 " the code action is not a symmetry: code {c} twists {} goes to code {c2} twists {}, {why}",1118 spell_turns(&t),1119 spell_turns(&u)1120 ),1121 None => println!(" the code action keeps every class on this base"),1122 }1123 }1124 println!(" {:.2} s", clock.elapsed().as_secs_f64());1125 }1126}11271128fn main() {1129 let verbs: Vec<String> = std::env::args().skip(1).collect();1130 let verbs = if verbs.is_empty() {1131 ["census"].iter().map(|s| s.to_string()).collect()1132 } else {1133 verbs1134 };1135 for verb in verbs {1136 let clock = Instant::now();1137 match verb.as_str() {1138 "census" => run_census(),1139 "curves" => run_curves(),1140 "junction" => run_junction(),1141 "named" => run_named(),1142 other => panic!("unknown verb {other}"),1143 }1144 println!(" {verb} {:.2} s", clock.elapsed().as_secs_f64());1145 }1146}11471148// WALKS11491150#[derive(Clone)]1151struct Walk {1152 turns: Vec<usize>,1153 radix: bool,1154 glue: Option<usize>,1155 end: Option<usize>,1156 plane: bool,1157 canonical: Option<usize>,1158}11591160impl Walk {1161 fn crossing(&self) -> Option<usize> {1162 match (self.glue, self.end) {1163 (Some(a), Some(b)) => Some(a.min(b)),1164 (a, b) => a.or(b),1165 }1166 }1167}11681169fn walks(lab: &Lab, k: usize) -> Vec<Vec<usize>> {1170 fn grow(lab: &Lab, k: usize, at: Z, turns: &mut Vec<usize>, out: &mut Vec<Vec<usize>>) {1171 let left = (k - turns.len()) as u64;1172 if lab.norm(sub(lab.b, at)) > left * left {1173 return;1174 }1175 if left == 0 {1176 out.push(turns.clone());1177 return;1178 }1179 for t in 0..lab.n {1180 turns.push(t);1181 grow(lab, k, add(at, lab.units[t]), turns, out);1182 turns.pop();1183 }1184 }1185 let mut out = Vec::new();1186 grow(lab, k, (0, 0), &mut Vec::new(), &mut out);1187 out1188}11891190fn images(lab: &Lab, turns: &[usize]) -> Vec<Vec<usize>> {1191 let mut out = vec![turns.to_vec()];1192 let mut back = turns.to_vec();1193 back.reverse();1194 out.push(back);1195 if let Some(eps) = lab.mirror() {1196 for w in out.clone() {1197 out.push(1198 w.iter()1199 .map(|&t| (lab.back(t) + lab.back(eps)) % lab.n)1200 .collect(),1201 );1202 }1203 }1204 out1205}12061207fn class_of(lab: &Lab, turns: &[usize]) -> Vec<usize> {1208 images(lab, turns).into_iter().min().unwrap()1209}12101211fn seen_by_census(lab: &Lab, design: &Design) -> Option<usize> {1212 let index = residue_index(lab);1213 for z in 0..lab.n {1214 let classes: Option<Vec<usize>> = design1215 .digits1216 .iter()1217 .map(|&d| index.get(&lab.turn(d, z)).copied())1218 .collect();1219 if let Some(c) = classes {1220 if c.windows(2).all(|w| w[0] < w[1]) {1221 return Some(c.iter().fold(0, |acc, &i| acc | 1 << i));1222 }1223 }1224 }1225 None1226}12271228fn classify(lab: &Lab, turns: &[usize]) -> Walk {1229 let design = walk_design(lab, turns);1230 assert!(lemma_step(lab, &design), "a walk outside the step lemma");1231 let m = machine(lab, &design.digits);1232 let plane = turns.len() == lab.residues.len() && m.glue(turns, true).is_none();1233 Walk {1234 turns: turns.to_vec(),1235 radix: incongruent(lab, &design.digits),1236 glue: m.glue(turns, false),1237 end: m.ending(turns),1238 plane,1239 canonical: seen_by_census(lab, &design),1240 }1241}12421243fn classify_all(lab: &Lab, list: &[Vec<usize>]) -> Vec<Walk> {1244 let next = std::sync::atomic::AtomicUsize::new(0);1245 let threads = std::thread::available_parallelism()1246 .map(|x| x.get())1247 .unwrap_or(4);1248 let chunk = 256;1249 let mut parts: Vec<Vec<(usize, Walk)>> = std::thread::scope(|scope| {1250 let handles: Vec<_> = (0..threads)1251 .map(|_| {1252 scope.spawn(|| {1253 let mut acc = Vec::new();1254 loop {1255 let i = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);1256 if i * chunk >= list.len() {1257 break;1258 }1259 for j in i * chunk..((i + 1) * chunk).min(list.len()) {1260 acc.push((j, classify(lab, &list[j])));1261 }1262 }1263 acc1264 })1265 })1266 .collect();1267 handles.into_iter().map(|h| h.join().unwrap()).collect()1268 });1269 let mut all: Vec<(usize, Walk)> = parts.drain(..).flatten().collect();1270 all.sort_by_key(|x| x.0);1271 all.into_iter().map(|x| x.1).collect()1272}12731274fn run_curves() {1275 println!("CURVES every walk of k unit steps from 0 to the base: digits the partial sums, twists the steps; every such design has unit steps at every level");1276 println!(" arc: the k^L + 1 vertices of level L distinct at every level; plane: k = N(b) and no two words with one place map; classes up to reversal and, where conj(b) is an associate, the mirror");1277 for (ring, value) in BASES {1278 let clock = Instant::now();1279 let lab = Lab::new(ring, value);1280 let q = lab.residues.len();1281 println!(1282 " {} base {} N(b) {q}",1283 ring_word(ring),1284 spell(ring, value)1285 );1286 println!(" k dim walks radix arcs arc-classes plane plane-arcs canonical deepest-crossing");1287 let mut listed = Vec::new();1288 for k in 2..=q {1289 let list = walks(&lab, k);1290 let all = classify_all(&lab, &list);1291 let radix: Vec<&Walk> = all.iter().filter(|w| w.radix).collect();1292 let arcs: Vec<&&Walk> = radix.iter().filter(|w| w.crossing().is_none()).collect();1293 let classes: HashSet<Vec<usize>> =1294 arcs.iter().map(|w| class_of(&lab, &w.turns)).collect();1295 let plane = radix.iter().filter(|w| w.plane).count();1296 let plane_arcs = arcs.iter().filter(|w| w.plane).count();1297 let canonical = radix.iter().filter(|w| w.canonical.is_some()).count();1298 let mut deep = (0usize, Vec::new());1299 for w in &radix {1300 if let Some(c) = w.crossing() {1301 if c > deep.0 || (c == deep.0 && w.turns < deep.1) {1302 deep = (c, w.turns.clone());1303 }1304 }1305 }1306 let other_arcs = all1307 .iter()1308 .filter(|w| !w.radix && w.crossing().is_none())1309 .count();1310 let widest = all.iter().filter_map(|w| w.crossing()).max().unwrap_or(0);1311 println!(1312 " {k} {:.4} {} {} {} {} {} {} {} level {} at {} (over every walk: {other_arcs} more arcs, deepest first revisit level {widest})",1313 2.0 * (k as f64).ln() / (q as f64).ln(),1314 all.len(),1315 radix.len(),1316 arcs.len(),1317 classes.len(),1318 plane,1319 plane_arcs,1320 canonical,1321 deep.0,1322 spell_turns(&deep.1)1323 );1324 let mut reps: Vec<Vec<usize>> = classes.into_iter().collect();1325 reps.sort();1326 for r in reps {1327 let w = classify(&lab, &r);1328 listed.push((k, r, w.plane, w.canonical));1329 }1330 }1331 for (k, r, plane, canonical) in &listed {1332 println!(1333 " arc class k {k}: turns {}{}{}",1334 spell_turns(r),1335 if *plane { ", plane-filling" } else { "" },1336 canonical.map_or(String::new(), |c| format!(1337 ", canonical code {c} after a turn"1338 ))1339 );1340 }1341 println!(" {:.2} s", clock.elapsed().as_secs_f64());1342 }1343}13441345// JUNCTION13461347const BRUTE: usize = 4096;13481349#[derive(Default)]1350struct Tests {1351 designs: u64,1352 checks: u64,1353 misses: u64,1354}13551356fn brute_top(k: usize) -> usize {1357 let mut top = 1;1358 while k.pow(top as u32 + 1) <= BRUTE {1359 top += 1;1360 }1361 top.max(2)1362}13631364fn brute_arc(lab: &Lab, design: &Design, level: usize) -> bool {1365 let mut p = points(lab, design, level);1366 p.push(lab.power(level));1367 p.iter().collect::<HashSet<_>>().len() == p.len()1368}13691370fn test_design(lab: &Lab, design: &Design, walk: bool, tests: &mut Tests) {1371 let t = &design.twists;1372 let m = machine(lab, &design.digits);1373 let broken = first_break(lab, design);1374 let (glue, maps, split) = (m.glue(t, false), m.glue(t, true), m.piece(t));1375 let end = if walk { m.ending(t) } else { None };1376 tests.designs += 1;1377 for level in 1..=brute_top(design.size()) {1378 let mut verdicts = vec![1379 (1380 brute_steps(lab, design, level),1381 broken.map_or(true, |b| level < b),1382 ),1383 (1384 brute_piece(lab, design, level),1385 split.map_or(true, |s| level < s),1386 ),1387 ];1388 let (glued, mapped) = brute_glue(lab, design, level);1389 verdicts.push((glued, glue.is_some_and(|g| g <= level)));1390 verdicts.push((mapped, maps.is_some_and(|g| g <= level)));1391 if walk {1392 let crossing = [glue, end].iter().flatten().min().copied();1393 verdicts.push((1394 brute_arc(lab, design, level),1395 crossing.map_or(true, |c| level < c),1396 ));1397 }1398 for (seen, said) in verdicts {1399 tests.checks += 1;1400 if seen != said {1401 tests.misses += 1;1402 }1403 }1404 }1405 if incongruent(lab, &design.digits) {1406 tests.checks += 1;1407 if design.radix(lab).words(3) != points(lab, design, 3) {1408 tests.misses += 1;1409 }1410 }1411}14121413fn test_symmetry(lab: &Lab, group: &[Symmetry], code: usize, t: &[usize]) -> bool {1414 let level = 3;1415 let mine = points(lab, &canonical(lab, code, t), level);1416 group.iter().all(|sym| {1417 let (c, u) = act(sym, code, t);1418 let theirs = points(lab, &canonical(lab, c, &u), level);1419 let moved: Vec<Z> = if sym.mirror {1420 mine.iter().map(|&z| lab.conj(z)).collect()1421 } else {1422 mine.clone()1423 };1424 same_up_to_turn(lab, &moved, &theirs)1425 })1426}14271428fn run_junction() {1429 println!("JUNCTION the step law, the difference automaton and the stabiliser against brute force over the words of each level");1430 println!(" per design and per level up to |F|^level <= {BRUTE}: unit steps, one piece, two words on one point, two words with one place map, and for walks the arc with its end vertex");1431 let mut total = Tests::default();1432 for (ring, value) in BASES {1433 let clock = Instant::now();1434 let lab = Lab::new(ring, value);1435 let q = lab.residues.len();1436 let group = stabiliser(&lab);1437 let mut all: Vec<(usize, u64)> = Vec::new();1438 for code in 0..1usize << q {1439 let k = code.count_ones();1440 if k >= 2 {1441 for index in 0..(lab.n as u64).pow(k) {1442 all.push((code, index));1443 }1444 }1445 }1446 let stride = (all.len() / 3000).max(1);1447 let mut tests = Tests::default();1448 let mut turned = 0u64;1449 for (i, &(code, index)) in all.iter().enumerate() {1450 let k = code.count_ones() as usize;1451 let t = unpack(lab.n, k, index);1452 let design = canonical(&lab, code, &t);1453 let path = unit_path(&lab, &design.digits) && lemma_step(&lab, &design);1454 if i % stride != 0 && !path {1455 continue;1456 }1457 test_design(&lab, &design, false, &mut tests);1458 if i % (stride * 4) == 0 {1459 assert!(1460 test_symmetry(&lab, &group, code, &t),1461 "the stabiliser moved a level set"1462 );1463 turned += 1;1464 }1465 }1466 let mut walked = 0;1467 for k in 2..=q {1468 for turns in walks(&lab, k) {1469 let design = walk_design(&lab, &turns);1470 if incongruent(&lab, &design.digits) {1471 test_design(&lab, &design, true, &mut tests);1472 walked += 1;1473 }1474 }1475 }1476 println!(1477 " {} base {}: {} canonical pairs of {} (every {stride}th and every step pair) and {walked} radix walks, {} level checks, {} misses; stabiliser images checked on {turned} pairs {:.2} s",1478 ring_word(ring),1479 spell(ring, value),1480 tests.designs - walked,1481 all.len(),1482 tests.checks,1483 tests.misses,1484 clock.elapsed().as_secs_f64()1485 );1486 assert_eq!(tests.misses, 0, "the lemma and the brute force disagree");1487 total.designs += tests.designs;1488 total.checks += tests.checks;1489 total.misses += tests.misses;1490 }1491 println!(1492 " all bases: {} designs, {} checks, {} misses",1493 total.designs, total.checks, total.misses1494 );1495 deep_witnesses();1496}14971498fn deep_witnesses() {1499 println!(" deep witnesses, brute force level by level");1500 let lab = Lab::new(Ring::Gaussian, (1, 1));1501 let held = canonical(&lab, 3, &[0, 2]);1502 let row: Vec<bool> = (1..=5).map(|l| brute_steps(&lab, &held, l)).collect();1503 println!(1504 " Z[i] base 1+i digits {} twists {}: unit steps at levels 1..5 {:?}, lemma first break level {:?}",1505 spell_list(lab.ring, &held.digits),1506 spell_turns(&held.twists),1507 row,1508 first_break(&lab, &held)1509 );1510 let lab = Lab::new(Ring::Gaussian, (2, 1));1511 let point = Design {1512 digits: vec![(0, 1), (1, 1), (1, 2)],1513 twists: vec![0, 1, 2],1514 };1515 let fixed = point1516 .digits1517 .iter()1518 .zip(point.twists.iter())1519 .all(|(&d, &t)| lab.mul(sub(lab.b, lab.units[t]), (1, 1)) == add(d, d));1520 println!(1521 " Z[i] base 2+i digits {} twists {}: incongruent {}, step law {}, unit steps at levels 1..8 {}, twists sum to {}, every map fixes (1+i)/2 {fixed}",1522 spell_list(lab.ring, &point.digits),1523 spell_turns(&point.twists),1524 incongruent(&lab, &point.digits),1525 lemma_step(&lab, &point),1526 (1..=8).all(|l| brute_steps(&lab, &point, l)),1527 spell(lab.ring, lab.sum(&point.twists))1528 );1529 let lab = Lab::new(Ring::Eisenstein, (2, 1));1530 let gloss = canonical(&lab, 6, &[0, 1]);1531 let level_two: Vec<u64> = points(&lab, &gloss, 2)1532 .windows(2)1533 .map(|w| lab.norm(sub(w[1], w[0])))1534 .collect();1535 println!(1536 " Z[w] base 1+w digits {} twists {}: junction identity {}, step norms at level 2 {:?}",1537 spell_list(lab.ring, &gloss.digits),1538 spell_turns(&gloss.twists),1539 junction_key(&lab, &gloss, 0) == (0, 0),1540 level_two1541 );1542 let middle = Design {1543 digits: vec![(-1, 0), (0, 0), (1, 0)],1544 twists: vec![0, 2, 0],1545 };1546 let m = machine(&lab, &middle.digits);1547 println!(1548 " Z[w] base 1+w digits {} twists {}: incongruent {}, step law {}, twists sum to {}, first revisit {:?}, first shared map {:?}, first split {:?}, unit steps and distinct points by brute force at levels 1..8 {}",1549 spell_list(lab.ring, &middle.digits),1550 spell_turns(&middle.twists),1551 incongruent(&lab, &middle.digits),1552 lemma_step(&lab, &middle),1553 spell(lab.ring, lab.sum(&middle.twists)),1554 m.glue(&middle.twists, false),1555 m.glue(&middle.twists, true),1556 m.piece(&middle.twists),1557 (1..=8).all(|l| brute_steps(&lab, &middle, l) && !brute_glue(&lab, &middle, l).0)1558 );1559 let walk = walk_design(&lab, &[0, 1]);1560 let row: Vec<bool> = (1..=6).map(|l| brute_arc(&lab, &walk, l)).collect();1561 let m = machine(&lab, &walk.digits);1562 println!(1563 " Z[w] base {} walk 01: arc at levels 1..6 {:?}, automaton first revisit level {:?}",1564 spell(lab.ring, lab.b),1565 row,1566 [m.glue(&walk.twists, false), m.ending(&walk.twists)]1567 .iter()1568 .flatten()1569 .min()1570 );1571 for (ring, value, code, turns, top) in [1572 (Ring::Gaussian, (2, 0), 14usize, vec![1usize, 3, 1], 6usize),1573 (Ring::Gaussian, (2, 1), 15, vec![0, 1, 1, 0], 7),1574 (Ring::Eisenstein, (2, 0), 11, vec![0, 3, 2], 10),1575 (Ring::Eisenstein, (2, 1), 6, vec![2, 4], 8),1576 (Ring::Eisenstein, (3, 0), 223, vec![0, 4, 3, 5, 2, 0, 4], 7),1577 (Ring::Eisenstein, (3, 1), 63, vec![4, 1, 1, 5, 2, 0], 7),1578 ] {1579 let lab = Lab::new(ring, value);1580 let design = canonical(&lab, code, &turns);1581 let m = machine(&lab, &design.digits);1582 let row: Vec<bool> = (1..=top).map(|l| brute_piece(&lab, &design, l)).collect();1583 println!(1584 " {} base {} code {code} twists {}: one piece at levels 1..{top} {:?}, automaton first split level {:?}",1585 ring_word(ring),1586 spell(ring, value),1587 spell_turns(&turns),1588 row,1589 m.piece(&turns)1590 );1591 }1592}15931594// NAMED15951596fn as_design(lab: &Lab, radix: &Radix) -> Design {1597 Design {1598 digits: radix.digits().to_vec(),1599 twists: radix1600 .twists()1601 .iter()1602 .map(|&u| lab.units.iter().position(|&v| v == u).unwrap())1603 .collect(),1604 }1605}16061607fn locate(name: &str, lab: &Lab, design: &Design) {1608 let t = &design.twists;1609 let m = machine(lab, &design.digits);1610 let step = lemma_step(lab, design);1611 let q = lab.residues.len();1612 let plane = design.size() == q && m.glue(t, true).is_none();1613 let distinct = m.glue(t, false).is_none();1614 let piece = m.piece(t).is_none();1615 let canonical = design.digits.iter().all(|d| lab.residues.contains(d));1616 let walk = step && design.digits[0] == (0, 0);1617 let crossing = if walk {1618 [m.glue(t, false), m.ending(t)]1619 .iter()1620 .flatten()1621 .min()1622 .copied()1623 } else {1624 None1625 };1626 println!(1627 " {name}: {} base {}, digits {}, twists {}, dim {:.6}",1628 ring_word(lab.ring),1629 spell(lab.ring, lab.b),1630 spell_list(lab.ring, &design.digits),1631 spell_turns(t),1632 2.0 * (design.size() as f64).ln() / (q as f64).ln()1633 );1634 println!(1635 " canonical digits {canonical}, seen by the canonical census {}, step {step}, piece {piece}, no two words on one point {distinct}, plane {plane}{}",1636 canonical && design.digits.windows(2).all(|w| {1637 lab.residues.iter().position(|&r| r == w[0]) < lab.residues.iter().position(|&r| r == w[1])1638 }),1639 if walk {1640 format!(1641 ", walk {} with first revisit {}",1642 spell_turns(t),1643 crossing.map_or("never, an arc".to_string(), |c| format!("at level {c}"))1644 )1645 } else if step {1646 String::new()1647 } else {1648 format!(", unit steps fail at level {}", first_break(lab, design).unwrap())1649 }1650 );1651}16521653fn hilbert(order: usize) -> Vec<Z> {1654 let side = 1i64 << order;1655 (0..side * side)1656 .map(|d| {1657 let (mut x, mut y, mut t) = (0i64, 0i64, d);1658 let mut s = 1i64;1659 while s < side {1660 let rx = 1 & (t / 2);1661 let ry = 1 & (t ^ rx);1662 if ry == 0 {1663 if rx == 1 {1664 x = s - 1 - x;1665 y = s - 1 - y;1666 }1667 std::mem::swap(&mut x, &mut y);1668 }1669 x += s * rx;1670 y += s * ry;1671 t /= 4;1672 s *= 2;1673 }1674 (x, y)1675 })1676 .collect()1677}16781679fn steps_of(list: &[Z]) -> Vec<Z> {1680 list.windows(2).map(|w| sub(w[1], w[0])).collect()1681}16821683fn run_named() {1684 println!("NAMED the classical curves located in the census; w = e^(i pi/3) in every printed element");1685 let lab = Lab::new(Ring::Eisenstein, (3, 0));1686 let k = as_design(&lab, &koch().unwrap());1687 assert_eq!(1688 k,1689 walk_design(&lab, &[0, 1, 5, 0]),1690 "the Koch design is not the walk 0150"1691 );1692 locate("Koch", &lab, &k);1693 let lab = Lab::new(Ring::Eisenstein, (2, 1));1694 let t = as_design(&lab, &terdragon().unwrap());1695 assert_eq!(1696 t,1697 canonical(&lab, 7, &[0, 2, 0]),1698 "the terdragon is not code 7 twisted 020"1699 );1700 assert_eq!(1701 t,1702 walk_design(&lab, &[0, 2, 0]),1703 "the terdragon is not the walk 020"1704 );1705 locate("terdragon", &lab, &t);1706 let lab = Lab::new(Ring::Gaussian, (1, 1));1707 let d = as_design(&lab, &twindragon().unwrap());1708 assert_eq!(1709 d,1710 canonical(&lab, 3, &[0, 0]),1711 "the twindragon is not code 3 untwisted"1712 );1713 locate("twindragon", &lab, &d);1714 locate("walk 01 at 1+i", &lab, &walk_design(&lab, &[0, 1]));1715 let lab = Lab::new(Ring::Eisenstein, (3, 1));1716 let d = as_design(&lab, &flowsnake().unwrap());1717 assert_eq!(1718 d,1719 canonical(&lab, 127, &[0; 7]),1720 "the flowsnake is not code 127 untwisted"1721 );1722 locate("flowsnake", &lab, &d);1723 let gosper = [0usize, 1, 3, 2, 0, 0, 5];1724 assert_eq!(1725 lab.sum(&gosper),1726 lab.b,1727 "the mirrored Gosper generator does not reach the base"1728 );1729 locate(1730 "Gosper generator 0132005 with every flag F",1731 &lab,1732 &walk_design(&lab, &gosper),1733 );1734 let s: Vec<Z> = gosper.iter().map(|&g| lab.units[g]).collect();1735 let flags = [false, true, true, false, false, false, true];1736 let copy = |j: usize| -> Vec<Z> {1737 let order: Vec<usize> = if flags[j] {1738 (0..7).rev().collect()1739 } else {1740 (0..7).collect()1741 };1742 order[..6].iter().map(|&i| lab.mul(s[j], s[i])).collect()1743 };1744 let inner: Vec<Z> = s[..6].to_vec();1745 let fits: Vec<usize> = (0..7)1746 .filter(|&j| {1747 (0..lab.n).any(|v| inner.iter().map(|&z| lab.turn(z, v)).collect::<Vec<_>>() == copy(j))1748 })1749 .collect();1750 println!(1751 " Gosper with flags FRRFFFR: the copies whose inner steps are a turn of the level-1 inner steps are {:?} of 0..6; a radix word order needs all seven",1752 fits1753 );1754 let lab = Lab::new(Ring::Gaussian, (2, 0));1755 let one: Vec<Z> = hilbert(1);1756 let two: Vec<Z> = hilbert(2);1757 let first = steps_of(&one);1758 let quarter = steps_of(&two[..4]);1759 let turns: Vec<usize> = (0..lab.n)1760 .filter(|&v| first.iter().map(|&z| lab.turn(z, v)).collect::<Vec<_>>() == quarter)1761 .collect();1762 println!(1763 " Hilbert: level 1 {} has steps {}, the first quarter of level 2 has steps {}; units carrying one to the other: {:?}",1764 spell_list(lab.ring, &one),1765 spell_list(lab.ring, &first),1766 spell_list(lab.ring, &quarter),1767 turns1768 );1769 let walk: Vec<usize> = first1770 .iter()1771 .map(|&z| lab.units.iter().position(|&u| u == z).unwrap())1772 .chain(std::iter::once(0))1773 .collect();1774 assert_eq!(lab.sum(&walk), lab.b);1775 locate(1776 "Hilbert level 1 as a walk with every flag F",1777 &lab,1778 &walk_design(&lab, &walk),1779 );1780}