main.rs
38.0 kB · rust · 1176 lines
1use mrlymath::bang::factory::create;2use mrlynum::gauss::Ring;3use mrlynum::radix::{flowsnake, gasket, koch, terdragon, tile, twindragon, Base, Radix};4use std::collections::{HashMap, HashSet};5use std::time::Instant;67type Point = (f64, f64);8type Map = (Point, Point);910const BUDGET: usize = 200_000;11const COMPARE: usize = 20_000;1213fn ring_word(ring: Ring) -> &'static str {14 match ring {15 Ring::Gaussian => "Z[i]",16 Ring::Eisenstein => "Z[omega]",17 }18}1920fn spell(ring: Ring, (a, b): (i64, i64)) -> String {21 let unit = match ring {22 Ring::Gaussian => "i",23 Ring::Eisenstein => "w",24 };25 match (a, b) {26 (a, 0) => format!("{a}"),27 (0, 1) => unit.to_string(),28 (0, -1) => format!("-{unit}"),29 (0, b) => format!("{b}{unit}"),30 (a, 1) => format!("{a}+{unit}"),31 (a, -1) => format!("{a}-{unit}"),32 (a, b) if b > 0 => format!("{a}+{b}{unit}"),33 (a, b) => format!("{a}{b}{unit}"),34 }35}3637fn complex_mul(z: (f64, f64), w: (f64, f64)) -> (f64, f64) {38 (z.0 * w.0 - z.1 * w.1, z.0 * w.1 + z.1 * w.0)39}4041fn complex_div(z: (f64, f64), w: (f64, f64)) -> (f64, f64) {42 let den = w.0 * w.0 + w.1 * w.1;43 ((z.0 * w.0 + z.1 * w.1) / den, (z.1 * w.0 - z.0 * w.1) / den)44}4546// IFS4748fn iterate(maps: &[Map], level: usize) -> Vec<(f64, f64)> {49 let mut out = vec![(0.0f64, 0.0f64)];50 for _ in 0..level {51 let mut next = Vec::with_capacity(out.len() * maps.len());52 for &(scale, shift) in maps {53 for &p in &out {54 let t = complex_mul(scale, p);55 next.push((t.0 + shift.0, t.1 + shift.1));56 }57 }58 out = next;59 }60 out61}6263fn koch_maps() -> Vec<Map> {64 let root = 3f64.sqrt();65 vec![66 ((1.0 / 3.0, 0.0), (0.0, 0.0)),67 ((1.0 / 6.0, root / 6.0), (1.0 / 3.0, 0.0)),68 ((1.0 / 6.0, -root / 6.0), (0.5, root / 6.0)),69 ((1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0)),70 ]71}7273fn gasket_maps() -> Vec<Map> {74 let root = 3f64.sqrt();75 let corners = [(1.0, 1.0), (3.0, 1.0), (2.0, 1.0 + root)];76 corners77 .iter()78 .map(|&(x, y)| ((0.5, 0.0), (x / 2.0, y / 2.0)))79 .collect()80}8182fn twindragon_maps() -> Vec<Map> {83 let scale = complex_div((1.0, 0.0), (1.0, 1.0));84 vec![(scale, (0.0, 0.0)), (scale, scale)]85}8687fn tile_maps(base: Base) -> Vec<Map> {88 let ring = base.ring();89 let b = ring.place(base.value().0, base.value().1);90 let scale = complex_div((1.0, 0.0), b);91 base.residues()92 .into_iter()93 .map(|d| {94 let place = ring.place(d.0, d.1);95 (scale, complex_div(place, b))96 })97 .collect()98}99100fn pin(got: &[(f64, f64)], want: &[(f64, f64)]) -> (f64, (f64, f64), f64) {101 let last = got.len() - 1;102 let dp = (got[last].0 - got[0].0, got[last].1 - got[0].1);103 let dt = (want[last].0 - want[0].0, want[last].1 - want[0].1);104 let den = dp.0 * dp.0 + dp.1 * dp.1;105 let scale = (dt.0 * dp.0 + dt.1 * dp.1) / den;106 let turn = (dt.1 * dp.0 - dt.0 * dp.1) / den;107 let shift = (want[0].0 - scale * got[0].0, want[0].1 - scale * got[0].1);108 (scale, shift, turn)109}110111fn deviation(got: &[(f64, f64)], want: &[(f64, f64)], scale: f64, shift: (f64, f64)) -> f64 {112 got.iter()113 .zip(want.iter())114 .map(|(a, b)| {115 let x = scale * a.0 + shift.0 - b.0;116 let y = scale * a.1 + shift.1 - b.1;117 (x * x + y * y).sqrt()118 })119 .fold(0.0f64, f64::max)120}121122fn top_level(size: usize, budget: usize) -> usize {123 let mut top = 1;124 while size.pow(top as u32 + 1) <= budget {125 top += 1;126 }127 top128}129130fn face_off(name: &str, design: &Radix, maps: &[Map], free: bool, note: &str) {131 let top = top_level(design.size(), COMPARE);132 let got = design.plane(top);133 let want = iterate(maps, top);134 assert_eq!(got.len(), want.len(), "{name} point counts differ");135 let (scale, shift, turn) = if free {136 pin(&got, &want)137 } else {138 (1.0, (0.0, 0.0), 0.0)139 };140 let worst = deviation(&got, &want, scale, shift);141 let motion = if free {142 format!(143 "scale {scale:.6} shift ({:.6}, {:.6}) turn {turn:.1e}",144 shift.0, shift.1145 )146 } else {147 "none".to_string()148 };149 println!(150 "{name} level {top} points {} worst {worst:.3e}",151 got.len()152 );153 println!(" motion allowed: {motion}");154 println!(" {note}");155}156157// TERDRAGON158159fn lsystem(level: usize) -> Vec<u8> {160 let mut out = b"F".to_vec();161 for _ in 0..level {162 let mut next = Vec::with_capacity(out.len() * 5);163 for &c in &out {164 if c == b'F' {165 next.extend_from_slice(b"F+F-F");166 } else {167 next.push(c);168 }169 }170 out = next;171 }172 out173}174175fn turtle(level: usize) -> Vec<(f64, f64)> {176 let root = 3f64.sqrt();177 let step = [(1.0, 0.0), (-0.5, root / 2.0), (-0.5, -root / 2.0)];178 let mut heading = 0usize;179 let mut at = (0.0f64, 0.0f64);180 let mut starts = Vec::new();181 for c in lsystem(level) {182 match c {183 b'F' => {184 starts.push(at);185 at = (at.0 + step[heading].0, at.1 + step[heading].1);186 }187 b'+' => heading = (heading + 1) % 3,188 _ => heading = (heading + 2) % 3,189 }190 }191 starts.into_iter().map(|p| complex_div(p, at)).collect()192}193194fn terdragon_reading(name: &str, design: &Radix, level: usize) -> f64 {195 let got = design.plane(level);196 let want = turtle(level);197 assert_eq!(got.len(), want.len(), "{name} point counts differ");198 let worst = deviation(&got, &want, 1.0, (0.0, 0.0));199 println!(200 "{name} twists {} level {level} points {} worst {worst:.3e}",201 design202 .twists()203 .iter()204 .map(|&u| spell(design.ring(), u))205 .collect::<Vec<_>>()206 .join(" "),207 got.len()208 );209 worst210}211212fn run_compare() {213 println!("COMPARE each design against an independent f64 iterated function system, seed 0");214 face_off(215 "gasket ",216 &gasket(),217 &gasket_maps(),218 true,219 "independent: the three ratio 1/2 similarities fixing an equilateral triangle placed at (1,1), (3,1), (2,1+sqrt 3); a translation and a positive scaling are allowed because that statement fixes the gasket only up to similarity, and the turn residual printed above is the check that no rotation was needed",220 );221 face_off(222 "koch ",223 &koch(),224 &koch_maps(),225 false,226 "independent in f64 only: z/3, e^(i pi/3) z/3 + 1/3, e^(-i pi/3) z/3 + 1/2 + i sqrt(3)/6, z/3 + 2/3 are the crate maps coefficient for coefficient, so this is a float self-check of exact ring arithmetic and no motion is allowed",227 );228 face_off(229 "twindragon",230 &twindragon(),231 &twindragon_maps(),232 false,233 "self-check by definition: z/(1+i) and (z+1)/(1+i) are the radix maps of base 1+i on its two residues, so no motion is allowed and the number is f64 round-off",234 );235 let seven = flowsnake();236 face_off(237 "tile7 ",238 &seven,239 &tile_maps(seven.base()),240 false,241 "self-check by definition: the seven maps (z + d)/(3+w) over the canonical residues are the radix maps of that base, so no motion is allowed; the flowsnake name is not tested here",242 );243 println!(" terdragon against the L-system F -> F + F - F at 120 degrees, three segments, unsourced reading");244 let level = 8;245 let twisted = terdragon_reading(" twisted ", &terdragon(), level);246 let plain = terdragon_reading(" untwisted", &terdragon().with_twists(&[0, 0, 0]), level);247 println!(248 " verdict: the twisted reading {} and the untwisted code 7 {}",249 if twisted < 1e-9 { "matches" } else { "misses" },250 if plain < 1e-9 { "matches" } else { "misses" }251 );252 assert!(twisted < 1e-9, "the twisted reading misses by {twisted}");253 assert!(plain > 1e-3, "the untwisted code 7 matches after all");254}255256// KOCH257258fn run_koch() {259 let design = koch();260 println!(261 "KOCH base {} on {}",262 spell(design.ring(), design.base().value()),263 ring_word(design.ring())264 );265 println!(266 " digits {} twists {}",267 design268 .digits()269 .iter()270 .map(|&d| spell(design.ring(), d))271 .collect::<Vec<_>>()272 .join(" "),273 design274 .twists()275 .iter()276 .map(|&u| spell(design.ring(), u))277 .collect::<Vec<_>>()278 .join(" ")279 );280 println!(281 " canonical digits {} class code {}",282 design.canonical(),283 design.code()284 );285 println!(" the four maps are phi_d coefficient for coefficient, so the column below is a float self-check");286 println!(" level words fill distinct worst");287 let maps = koch_maps();288 for level in 1..=5 {289 let got = design.plane(level);290 let want = iterate(&maps, level);291 let worst = deviation(&got, &want, 1.0, (0.0, 0.0));292 println!(293 " {level:>5} {:>5} {:>4} {:>8} {worst:.3e}",294 got.len(),295 design.fill(level),296 design.distinct(level)297 );298 assert_eq!(got.len() as u128, design.fill(level));299 assert_eq!(design.fill(level), 4u128.pow(level as u32));300 assert!(worst < 1e-9, "level {level} misses by {worst}");301 }302}303304// NAMED305306fn report(name: &str, design: &Radix) {307 let ring = design.ring();308 let q = design.base().norm();309 let size = design.size();310 let top = top_level(size, BUDGET);311 let counts: Vec<String> = (1..=top)312 .map(|level| design.distinct(level).to_string())313 .collect();314 let pairwise = (1..=top).all(|level| design.distinct(level) as u128 == design.fill(level));315 println!(316 "{name} {} base {} q {q} code {} |F| {size} fill |F|^L dim {:.6}",317 ring_word(ring),318 spell(ring, design.base().value()),319 design.code(),320 design.dimension()321 );322 println!(323 " digits {}",324 design325 .digits()326 .iter()327 .map(|&d| spell(ring, d))328 .collect::<Vec<_>>()329 .join(" ")330 );331 println!(" distinct levels 1..{top}: {}", counts.join(" "));332 println!(" pairwise distinct {pairwise}");333}334335fn run_named() {336 println!("NAMED codes only; the names are tested by the verb compare and nowhere here");337 report("code 7 at 2 ", &gasket());338 report("code 3 at 1+i ", &twindragon());339 report("code 7 at 2+w ", &terdragon().with_twists(&[0, 0, 0]));340 report("code 127 at 3+w ", &flowsnake());341 report("code 147 at 3 ", &koch());342 let glue = Radix::from_code(Base::new(Ring::Gaussian, (2, 0)), 3).with_twists(&[0, 2]);343 println!(344 "twist glue Z[i] base 2 q 4 code 3 |F| 2 twists 1 -1 dim {:.6}",345 glue.dimension()346 );347 println!(" words scaled by b^2: {:?}", glue.words(2));348 let top = 16;349 println!(350 " fill levels 1..{top}: {}",351 (1..=top)352 .map(|l| glue.fill(l).to_string())353 .collect::<Vec<_>>()354 .join(" ")355 );356 let counts: Vec<String> = (1..=top)357 .map(|level| glue.distinct(level).to_string())358 .collect();359 println!(" distinct levels 1..{top}: {}", counts.join(" "));360 for level in 1..=top {361 assert_eq!(362 glue.distinct(level),363 (1usize << (level - 1)) + 1,364 "the distinct count leaves 2^(L-1) + 1 at level {level}"365 );366 }367 println!(" distinct equals 2^(L-1) + 1 at every level to {top}");368}369370// TODAY371372fn run_today() {373 let level = 2;374 let mut checked = 0usize;375 let mut missed = 0usize;376 let mut first: Option<String> = None;377 for m in [2u64, 3u64] {378 let cells = (m * m) as usize;379 for code in 0..(1u128 << cells) {380 let design = tile(m, code);381 let got: HashSet<(i64, i64)> = design.words(level).into_iter().collect();382 let tensor = create(code, m as usize, 2, m as usize, level).unwrap();383 let side = m.pow(level as u32) as usize;384 let mut want = HashSet::new();385 for row in 0..side {386 for col in 0..side {387 if tensor.bytes()[row * side + col] == 1 {388 want.insert((col as i64, row as i64));389 }390 }391 }392 checked += 1;393 assert_eq!(design.words(level).len() as u128, design.fill(level));394 if got != want {395 missed += 1;396 if first.is_none() {397 first = Some(format!("m {m} code {code}"));398 }399 }400 }401 }402 println!("TODAY level {level}, bases 2 and 3, every plane code");403 println!(" codes checked {checked} mismatches {missed}");404 println!(405 " the code is read in box row-major order, bit r m + c, not in canonical residue order"406 );407 for m in [2u64, 3u64] {408 let full = tile(m, (1u128 << (m * m)) - 1);409 let base = Base::new(Ring::Gaussian, (m as i64, 0));410 println!(411 " m {m}: box residues are the canonical system {}",412 full.canonical()413 );414 println!(415 " box {}",416 full.digits()417 .iter()418 .map(|&d| spell(Ring::Gaussian, d))419 .collect::<Vec<_>>()420 .join(" ")421 );422 println!(423 " canonical {}",424 base.residues()425 .iter()426 .map(|&z| spell(Ring::Gaussian, z))427 .collect::<Vec<_>>()428 .join(" ")429 );430 }431 match first {432 None => println!(" first mismatch none"),433 Some(where_) => println!(" first mismatch {where_}"),434 }435}436437// CENSUS438439fn cycles(map: &[usize]) -> u32 {440 let mut seen = vec![false; map.len()];441 let mut count = 0;442 for start in 0..map.len() {443 if seen[start] {444 continue;445 }446 count += 1;447 let mut at = start;448 while !seen[at] {449 seen[at] = true;450 at = map[at];451 }452 }453 count454}455456fn burnside(group: &[Vec<usize>]) -> u128 {457 let total: u128 = group.iter().map(|map| 1u128 << cycles(map)).sum();458 assert_eq!(total % group.len() as u128, 0, "Burnside is not an integer");459 total / group.len() as u128460}461462fn image(group: &[Vec<usize>]) -> usize {463 group.iter().collect::<HashSet<_>>().len()464}465466fn act(map: &[usize], code: usize) -> usize {467 let mut out = 0usize;468 for (i, &j) in map.iter().enumerate() {469 if (code >> i) & 1 == 1 {470 out |= 1 << j;471 }472 }473 out474}475476fn orbit(group: &[Vec<usize>], code: usize) -> Vec<usize> {477 let mut seen = HashSet::new();478 let mut stack = vec![code];479 seen.insert(code);480 while let Some(at) = stack.pop() {481 for map in group {482 let next = act(map, at);483 if seen.insert(next) {484 stack.push(next);485 }486 }487 }488 let mut out: Vec<usize> = seen.into_iter().collect();489 out.sort_unstable();490 out491}492493fn walk(group: &[Vec<usize>], q: usize) -> u128 {494 let mut seen = vec![false; 1usize << q];495 let mut orbits = 0u128;496 for code in 0..(1usize << q) {497 if seen[code] {498 continue;499 }500 orbits += 1;501 for member in orbit(group, code) {502 seen[member] = true;503 }504 }505 orbits506}507508fn run_census() {509 let bases = [510 (Ring::Gaussian, (2i64, 0i64)),511 (Ring::Gaussian, (1, 1)),512 (Ring::Gaussian, (2, 1)),513 (Ring::Eisenstein, (2, 0)),514 (Ring::Eisenstein, (2, 1)),515 (Ring::Eisenstein, (3, 0)),516 (Ring::Eisenstein, (3, 1)),517 ];518 println!(519 "CENSUS classes of digit CODES under the residue action, never designs up to similarity"520 );521 println!(" ring base q abstract image mirror codes classes walk");522 for (ring, value) in bases {523 let base = Base::new(ring, value);524 let q = base.norm() as usize;525 let group = base.group();526 let count = burnside(&group);527 let seen = walk(&group, q);528 assert_eq!(count, seen, "Burnside and the orbit walk disagree");529 println!(530 " {:>9} {:>4} {q} {:>8} {:>5} {:>6} {:>5} {count:>7} {seen:>4}",531 ring_word(ring),532 spell(ring, value),533 group.len(),534 image(&group),535 base.mirrored(),536 1u128 << q537 );538 println!(539 " residues {}",540 base.residues()541 .iter()542 .map(|&z| spell(ring, z))543 .collect::<Vec<_>>()544 .join(" ")545 );546 }547 println!(" abstract is the order of R^* semidirect <conj>, image the order it acts through on the residues");548 let base = Base::new(Ring::Eisenstein, (3, 0));549 let q = base.norm() as usize;550 let twisted: u128 = (0..1u128 << q)551 .map(|code| 6u128.pow(code.count_ones()))552 .sum();553 println!(554 " base 3 on Z[omega]: {} codes in {} classes of codes; over the codes, not the classes, the twist vectors number sum_k binom(9,k) 6^k = {twisted} = 7^{q}",555 1u128 << q,556 burnside(&base.group())557 );558 assert_eq!(twisted, 7u128.pow(q as u32));559 let mut sizes: HashMap<usize, usize> = HashMap::new();560 for code in 0..1u128 << q {561 *sizes.entry(code.count_ones() as usize).or_default() += 1;562 }563 let mut keys: Vec<usize> = sizes.keys().copied().collect();564 keys.sort_unstable();565 println!(566 " codes by |F|: {}",567 keys.iter()568 .map(|k| format!("{k}:{}", sizes[k]))569 .collect::<Vec<_>>()570 .join(" ")571 );572}573574// AFFINE575576#[derive(Clone, Copy, Debug, PartialEq, Eq)]577struct Alg {578 p: i128,579 q: i128,580 r: i128,581 m: i128,582}583584fn gcd(a: i128, b: i128) -> i128 {585 if b == 0 {586 a.abs().max(1)587 } else {588 gcd(b, a % b)589 }590}591592fn twist_of(ring: Ring) -> i128 {593 match ring {594 Ring::Gaussian => 0,595 Ring::Eisenstein => 1,596 }597}598599impl Alg {600 fn new(m: i128, p: i128, q: i128, r: i128) -> Alg {601 assert!(r != 0, "a rational needs a nonzero denominator");602 let (p, q, r) = if r < 0 { (-p, -q, -r) } else { (p, q, r) };603 if p == 0 && q == 0 {604 return Alg {605 p: 0,606 q: 0,607 r: 1,608 m,609 };610 }611 let g = gcd(gcd(p, q), r);612 Alg {613 p: p / g,614 q: q / g,615 r: r / g,616 m,617 }618 }619 fn whole(m: i128, p: i64, q: i64) -> Alg {620 Alg::new(m, p as i128, q as i128, 1)621 }622 fn sub(self, other: Alg) -> Alg {623 Alg::new(624 self.m,625 self.p * other.r - other.p * self.r,626 self.q * other.r - other.q * self.r,627 self.r * other.r,628 )629 }630 fn add(self, other: Alg) -> Alg {631 Alg::new(632 self.m,633 self.p * other.r + other.p * self.r,634 self.q * other.r + other.q * self.r,635 self.r * other.r,636 )637 }638 fn mul(self, other: Alg) -> Alg {639 Alg::new(640 self.m,641 self.p * other.p - self.q * other.q,642 self.p * other.q + self.q * other.p - self.m * self.q * other.q,643 self.r * other.r,644 )645 }646 fn conjugate(self) -> Alg {647 Alg::new(self.m, self.p - self.m * self.q, -self.q, self.r)648 }649 fn div(self, other: Alg) -> Alg {650 let n = other.p * other.p - self.m * other.p * other.q + other.q * other.q;651 assert!(n != 0, "no division by zero");652 self.mul(Alg::new(653 self.m,654 other.r * (other.p - self.m * other.q),655 -other.r * other.q,656 n,657 ))658 }659}660661fn direct(left: &[Alg], right: &[Alg]) -> bool {662 if left.len() != right.len() {663 return false;664 }665 if left.len() < 2 {666 return true;667 }668 let span = left[1].sub(left[0]);669 for first in 0..right.len() {670 for second in 0..right.len() {671 if first == second {672 continue;673 }674 let v = right[second].sub(right[first]).div(span);675 let t = right[first].sub(v.mul(left[0]));676 if left.iter().all(|&d| right.contains(&v.mul(d).add(t))) {677 return true;678 }679 }680 }681 false682}683684fn conjugate_sets(left: &[Alg], right: &[Alg], mirror: bool) -> bool {685 if direct(left, right) {686 return true;687 }688 if !mirror {689 return false;690 }691 let flipped: Vec<Alg> = left.iter().map(|d| d.conjugate()).collect();692 direct(&flipped, right)693}694695fn digits_of(ring: Ring, residues: &[(i64, i64)], code: usize) -> Vec<Alg> {696 let m = twist_of(ring);697 residues698 .iter()699 .enumerate()700 .filter(|(i, _)| (code >> i) & 1 == 1)701 .map(|(_, &z)| Alg::whole(m, z.0, z.1))702 .collect()703}704705fn similar_classes(706 ring: Ring,707 residues: &[(i64, i64)],708 codes: &[usize],709 mirror: bool,710) -> Vec<usize> {711 let sets: Vec<Vec<Alg>> = codes712 .iter()713 .map(|&c| digits_of(ring, residues, c))714 .collect();715 let mut label = vec![usize::MAX; codes.len()];716 let mut classes = 0;717 for i in 0..codes.len() {718 if label[i] != usize::MAX {719 continue;720 }721 label[i] = classes;722 for j in i + 1..codes.len() {723 if label[j] == usize::MAX && conjugate_sets(&sets[i], &sets[j], mirror) {724 label[j] = classes;725 }726 }727 classes += 1;728 }729 label730}731732fn orbit_labels(group: &[Vec<usize>], codes: &[usize]) -> (Vec<usize>, usize) {733 let mut seat: HashMap<usize, usize> = HashMap::new();734 let mut orbits = 0usize;735 for &code in codes {736 if seat.contains_key(&code) {737 continue;738 }739 for member in orbit(group, code) {740 seat.insert(member, orbits);741 }742 orbits += 1;743 }744 (codes.iter().map(|c| seat[c]).collect(), orbits)745}746747fn crossing(748 orbits: &[usize],749 classes: &[usize],750 codes: &[usize],751) -> (Option<(usize, usize)>, Option<(usize, usize)>) {752 let mut split = None;753 let mut merge = None;754 for i in 0..codes.len() {755 for j in i + 1..codes.len() {756 let same_orbit = orbits[i] == orbits[j];757 let same_class = classes[i] == classes[j];758 if same_orbit && !same_class && split.is_none() {759 split = Some((codes[i], codes[j]));760 }761 if !same_orbit && same_class && merge.is_none() {762 merge = Some((codes[i], codes[j]));763 }764 }765 }766 (split, merge)767}768769fn pair_word(pair: Option<(usize, usize)>) -> String {770 match pair {771 Some((a, b)) => format!("{a}/{b}"),772 None => "-".to_string(),773 }774}775776fn spell_code(ring: Ring, residues: &[(i64, i64)], code: usize) -> String {777 format!(778 "{code} ({})",779 residues780 .iter()781 .enumerate()782 .filter(|(i, _)| (code >> i) & 1 == 1)783 .map(|(_, &z)| spell(ring, z))784 .collect::<Vec<_>>()785 .join(", ")786 )787}788789fn grid(z: (f64, f64)) -> (i64, i64) {790 ((z.0 * 1e9).round() as i64, (z.1 * 1e9).round() as i64)791}792793fn shape(points: &[(f64, f64)]) -> Vec<(i64, i64)> {794 if points.len() < 2 {795 return Vec::new();796 }797 let mut best: Option<Vec<(i64, i64)>> = None;798 for first in 0..points.len() {799 for second in 0..points.len() {800 if first == second {801 continue;802 }803 let span = (804 points[second].0 - points[first].0,805 points[second].1 - points[first].1,806 );807 let mut key: Vec<(i64, i64)> = points808 .iter()809 .map(|p| {810 grid(complex_div(811 (p.0 - points[first].0, p.1 - points[first].1),812 span,813 ))814 })815 .collect();816 key.sort_unstable();817 if best.as_ref().is_none_or(|held| key < *held) {818 best = Some(key);819 }820 }821 }822 best.expect("two digits fix a shape")823}824825fn float_similar(ring: Ring, residues: &[(i64, i64)], codes: &[usize], mirror: bool) -> usize {826 let mut seen = HashSet::new();827 for &code in codes {828 let points: Vec<(f64, f64)> = residues829 .iter()830 .enumerate()831 .filter(|(i, _)| (code >> i) & 1 == 1)832 .map(|(_, &z)| ring.place(z.0, z.1))833 .collect();834 let mut key = shape(&points);835 if mirror {836 let flipped: Vec<(f64, f64)> = points.iter().map(|&(x, y)| (x, -y)).collect();837 key = key.min(shape(&flipped));838 }839 seen.insert(key);840 }841 seen.len()842}843844#[derive(Clone, Copy, Debug, Hash, PartialEq, Eq, PartialOrd, Ord)]845struct Frac {846 n: i128,847 d: i128,848}849850impl Frac {851 fn new(n: i128, d: i128) -> Frac {852 assert!(d != 0, "a fraction needs a nonzero denominator");853 let (n, d) = if d < 0 { (-n, -d) } else { (n, d) };854 if n == 0 {855 return Frac { n: 0, d: 1 };856 }857 let g = gcd(n, d);858 Frac { n: n / g, d: d / g }859 }860}861862fn det(u: (i64, i64), v: (i64, i64)) -> i128 {863 u.0 as i128 * v.1 as i128 - u.1 as i128 * v.0 as i128864}865866fn lattice_of(residues: &[(i64, i64)], code: usize) -> Vec<(i64, i64)> {867 residues868 .iter()869 .enumerate()870 .filter(|(i, _)| (code >> i) & 1 == 1)871 .map(|(_, &z)| z)872 .collect()873}874875fn affine_key(points: &[(i64, i64)]) -> Vec<(Frac, Frac)> {876 let n = points.len();877 if n < 2 {878 return Vec::new();879 }880 let sub = |a: (i64, i64), b: (i64, i64)| (a.0 - b.0, a.1 - b.1);881 let mut best: Option<Vec<(Frac, Frac)>> = None;882 for i in 0..n {883 for j in 0..n {884 for k in 0..n {885 if i == j || i == k || j == k {886 continue;887 }888 let e1 = sub(points[j], points[i]);889 let e2 = sub(points[k], points[i]);890 let base = det(e1, e2);891 if base == 0 {892 continue;893 }894 let mut key: Vec<(Frac, Frac)> = points895 .iter()896 .map(|&p| {897 let u = sub(p, points[i]);898 (Frac::new(det(u, e2), base), Frac::new(det(e1, u), base))899 })900 .collect();901 key.sort_unstable();902 if best.as_ref().is_none_or(|held| key < *held) {903 best = Some(key);904 }905 }906 }907 }908 if let Some(key) = best {909 return key;910 }911 for i in 0..n {912 for j in 0..n {913 if i == j {914 continue;915 }916 let v = sub(points[j], points[i]);917 let mut key: Vec<(Frac, Frac)> = points918 .iter()919 .map(|&p| {920 let u = sub(p, points[i]);921 let along = if v.0 != 0 {922 Frac::new(u.0 as i128, v.0 as i128)923 } else {924 Frac::new(u.1 as i128, v.1 as i128)925 };926 (along, Frac::new(0, 1))927 })928 .collect();929 key.sort_unstable();930 if best.as_ref().is_none_or(|held| key < *held) {931 best = Some(key);932 }933 }934 }935 best.expect("two distinct digits fix a line")936}937938fn plane_classes(residues: &[(i64, i64)], codes: &[usize]) -> Vec<usize> {939 let mut seat: HashMap<Vec<(Frac, Frac)>, usize> = HashMap::new();940 let mut out = Vec::with_capacity(codes.len());941 for &code in codes {942 let key = affine_key(&lattice_of(residues, code));943 let next = seat.len();944 out.push(*seat.entry(key).or_insert(next));945 }946 out947}948949fn carry(matrix: [[i64; 2]; 2], points: &[(i64, i64)]) -> Vec<(i64, i64)> {950 points951 .iter()952 .map(|&(a, b)| {953 (954 matrix[0][0] * a + matrix[0][1] * b,955 matrix[1][0] * a + matrix[1][1] * b,956 )957 })958 .collect()959}960961fn same_set(left: &[(i64, i64)], right: &[(i64, i64)]) -> bool {962 left.len() == right.len() && left.iter().all(|z| right.contains(z))963}964965fn run_affine() {966 let mut below = 0usize;967 let mut equal = 0usize;968 let mut above = 0usize;969 let mut identical = 0usize;970 let mut crossed = 0usize;971 let mut affine_below = 0usize;972 let mut affine_equal = 0usize;973 let mut affine_above = 0usize;974 let mut affine_crossed = 0usize;975 let bases = [976 (Ring::Gaussian, (2i64, 0i64)),977 (Ring::Gaussian, (1, 1)),978 (Ring::Gaussian, (2, 1)),979 (Ring::Eisenstein, (2, 0)),980 (Ring::Eisenstein, (2, 1)),981 (Ring::Eisenstein, (3, 0)),982 (Ring::Eisenstein, (3, 1)),983 ];984 println!("AFFINE the untwisted canonical digit sets of every base up to SIMILARITY and up to AFFINE conjugacy, beside the code census, in exact arithmetic");985 println!(" conjugating the place maps by an invertible real affine h(x) = H x + s gives (y + H d + s (b - 1))/b, again an untwisted base-b place map exactly when H commutes with multiplication by 1/b, and s (b - 1) sweeps the plane since N(b) >= 2 forces b != 1");986 println!(" at a non-real base the centraliser of 1/b in the two by two real matrices is C, so the conjugacy group is the similarity group x -> v x + t; at a real base 1/b is the scalar (1/b) I, it commutes with every H, and the conjugacy group is the whole real affine group GL_2 semidirect R^2");987 println!(" the mirror x -> v conj(x) + t preserves the untwisted base-b family exactly when conj(b) = b, since a direct conjugacy keeps the derivative 1/b and a mirror one sends it to 1/conj(b); at a real base it is one element of the full affine group and not the only new one");988 println!(" simil is the class count under the similarity group, affine the class count under the conjugacy group, equal to simil at the four non-real bases by the centraliser lemma and computed over GL_2(Q) semidirect Q^2 at the three real bases");989 println!(" orbits counts digit CODES under the unit group joined by conjugation where conj(b) is an associate of b; ssplit and smerge witness the crossing of orbits with simil, asplit and amerge the crossing of orbits with affine");990 println!(991 " ring base q |F| codes orbits simil affine ssplit smerge asplit amerge"992 );993 for (ring, value) in bases {994 let base = Base::new(ring, value);995 let residues = base.residues();996 let q = residues.len();997 let group = base.group();998 let real = value.1 == 0;999 let mut codes_seen = 0u128;1000 let mut orbits_seen = 0usize;1001 let mut simil_seen = 0usize;1002 let mut affine_seen = 0usize;1003 for size in 0..=q {1004 let codes: Vec<usize> = (0..1usize << q)1005 .filter(|c| c.count_ones() as usize == size)1006 .collect();1007 let (orbits, orbit_count) = orbit_labels(&group, &codes);1008 let simil = similar_classes(ring, &residues, &codes, real);1009 let simil_count = simil.iter().copied().max().unwrap() + 1;1010 assert_eq!(1011 float_similar(ring, &residues, &codes, real),1012 simil_count,1013 "the float rerun of the normal form disagrees with the exact similarity classes"1014 );1015 let affine = if real {1016 plane_classes(&residues, &codes)1017 } else {1018 simil.clone()1019 };1020 let affine_count = affine.iter().copied().max().unwrap() + 1;1021 for i in 0..codes.len() {1022 for j in i + 1..codes.len() {1023 assert!(1024 simil[i] != simil[j] || affine[i] == affine[j],1025 "the similarity classes do not refine the affine classes"1026 );1027 }1028 }1029 let (split, merge) = crossing(&orbits, &simil, &codes);1030 let (asplit, amerge) = crossing(&orbits, &affine, &codes);1031 println!(1032 " {:>9} {:>4} {q} {size:>3} {:>5} {orbit_count:>6} {simil_count:>5} {affine_count:>6} {:>7} {:>7} {:>7} {:>7}",1033 ring_word(ring),1034 spell(ring, value),1035 codes.len(),1036 pair_word(split),1037 pair_word(merge),1038 pair_word(asplit),1039 pair_word(amerge)1040 );1041 match simil_count.cmp(&orbit_count) {1042 std::cmp::Ordering::Less => below += 1,1043 std::cmp::Ordering::Equal => equal += 1,1044 std::cmp::Ordering::Greater => above += 1,1045 }1046 match affine_count.cmp(&orbit_count) {1047 std::cmp::Ordering::Less => affine_below += 1,1048 std::cmp::Ordering::Equal => affine_equal += 1,1049 std::cmp::Ordering::Greater => affine_above += 1,1050 }1051 if split.is_none() && merge.is_none() {1052 identical += 1;1053 }1054 if split.is_some() && merge.is_some() {1055 crossed += 1;1056 }1057 if asplit.is_some() && amerge.is_some() {1058 affine_crossed += 1;1059 }1060 codes_seen += codes.len() as u128;1061 orbits_seen += orbit_count;1062 simil_seen += simil_count;1063 affine_seen += affine_count;1064 }1065 println!(1066 " totals: codes {codes_seen} = 2^{q}, code classes {orbits_seen}, similarity classes {simil_seen}, affine classes {affine_seen}, real base {real}, mirror in the code group {}",1067 base.mirrored()1068 );1069 assert_eq!(codes_seen, 1u128 << q, "the sizes do not exhaust the codes");1070 assert_eq!(1071 orbits_seen as u128,1072 burnside(&group),1073 "the per-size orbit counts do not sum to the census, which controls the code column alone"1074 );1075 }1076 println!(1077 " over the {} cells the similarity count is below the code count in {below}, equal in {equal} and above it in {above}, the two partitions identical in {identical} and crossing in {crossed}",1078 below + equal + above1079 );1080 println!(1081 " over the same cells the affine count is below the code count in {affine_below}, equal in {affine_equal} and above it in {affine_above}, and the two partitions cross in {affine_crossed}"1082 );1083 let base = Base::new(Ring::Eisenstein, (3, 0));1084 let residues = base.residues();1085 let group = base.group();1086 println!(1087 " the fixed witnesses at base 3 on Z[omega], |F| = 3, untwisted canonical digit sets"1088 );1089 for (a, b, note) in [1090 (1091 131usize,1092 137usize,1093 "share a census orbit and are not similar",1094 ),1095 (7, 42, "are similar and sit in different census orbits"),1096 ] {1097 let codes = vec![a, b];1098 let (orbits, _) = orbit_labels(&group, &codes);1099 let simil = similar_classes(Ring::Eisenstein, &residues, &codes, true);1100 println!(1101 " codes {} and {} {note}: same orbit {}, same similarity class {}",1102 spell_code(Ring::Eisenstein, &residues, a),1103 spell_code(Ring::Eisenstein, &residues, b),1104 orbits[0] == orbits[1],1105 simil[0] == simil[1]1106 );1107 assert_eq!(orbits[0] == orbits[1], note.starts_with("share"));1108 assert_eq!(simil[0] == simil[1], !note.starts_with("share"));1109 }1110 for (a, b, matrix) in [1111 (131usize, 137usize, [[0i64, 2], [1, -1]]),1112 (7, 131, [[1, 1], [0, 1]]),1113 ] {1114 let left = lattice_of(&residues, a);1115 let right = lattice_of(&residues, b);1116 let moved = carry(matrix, &left);1117 let turn = det((matrix[0][0], matrix[1][0]), (matrix[0][1], matrix[1][1]));1118 println!(1119 " codes {} and {} are affinely conjugate by H = [[{}, {}], [{}, {}]] of determinant {turn}, same affine key {}",1120 spell_code(Ring::Eisenstein, &residues, a),1121 spell_code(Ring::Eisenstein, &residues, b),1122 matrix[0][0],1123 matrix[0][1],1124 matrix[1][0],1125 matrix[1][1],1126 affine_key(&left) == affine_key(&right)1127 );1128 assert!(turn != 0, "a conjugacy needs an invertible H");1129 assert!(same_set(&moved, &right), "H does not carry the digit set");1130 assert_eq!(1131 affine_key(&left),1132 affine_key(&right),1133 "the affine normal form misses a witnessed conjugacy"1134 );1135 }1136 let three: Vec<usize> = (0..512usize).filter(|c| c.count_ones() == 3).collect();1137 let (_, orbit_count) = orbit_labels(&group, &three);1138 let simil = similar_classes(Ring::Eisenstein, &residues, &three, true);1139 let affine = plane_classes(&residues, &three);1140 assert_eq!(three.len(), 84);1141 assert_eq!(orbit_count, 13);1142 assert_eq!(simil.iter().copied().max().unwrap() + 1, 9);1143 assert_eq!(affine.iter().copied().max().unwrap() + 1, 2);1144 let bare = similar_classes(Ring::Eisenstein, &residues, &three, false);1145 println!(1146 " the mirror is load-bearing for similarity: the 84 three-digit codes fall in {} direct classes and {} once the mirror joins, and in {} affine classes, the collinear triples against the rest",1147 bare.iter().copied().max().unwrap() + 1,1148 simil.iter().copied().max().unwrap() + 1,1149 affine.iter().copied().max().unwrap() + 11150 );1151}11521153fn main() {1154 let verbs: Vec<String> = std::env::args().skip(1).collect();1155 let verbs = if verbs.is_empty() {1156 ["koch", "compare", "named", "today", "census", "affine"]1157 .iter()1158 .map(|s| s.to_string())1159 .collect()1160 } else {1161 verbs1162 };1163 for verb in verbs {1164 let clock = Instant::now();1165 match verb.as_str() {1166 "koch" => run_koch(),1167 "compare" => run_compare(),1168 "named" => run_named(),1169 "today" => run_today(),1170 "census" => run_census(),1171 "affine" => run_affine(),1172 other => panic!("unknown verb {other}"),1173 }1174 println!(" {verb} {:.2} s", clock.elapsed().as_secs_f64());1175 }1176}