shape.rs
20.6 kB · rust · 702 lines
1use crate::design::{plane, BASE};2use crate::orbit;3use std::collections::HashMap;45pub struct Classes {6 pub index: Vec<u32>,7 pub count: usize,8 pub cells: usize,9}1011pub fn d4(n: usize) -> Vec<Vec<usize>> {12 let mut out: Vec<Vec<usize>> = Vec::new();13 for flip in 0..2 {14 for turn in 0..4 {15 let map: Vec<usize> = (0..n * n)16 .map(|f| {17 let (mut r, mut c) = (f / n, f % n);18 if flip == 1 {19 std::mem::swap(&mut r, &mut c);20 }21 for _ in 0..turn {22 let next = (c, n - 1 - r);23 r = next.0;24 c = next.1;25 }26 r * n + c27 })28 .collect();29 if !out.contains(&map) {30 out.push(map);31 }32 }33 }34 out35}3637pub fn oct(n: usize) -> Vec<Vec<usize>> {38 let axes = [39 [0, 1, 2],40 [0, 2, 1],41 [1, 0, 2],42 [1, 2, 0],43 [2, 0, 1],44 [2, 1, 0],45 ];46 let mut out: Vec<Vec<usize>> = Vec::new();47 for axis in axes {48 for signs in 0..8u32 {49 let map: Vec<usize> = (0..n * n * n)50 .map(|f| {51 let p = [f / (n * n), (f / n) % n, f % n];52 let mut q = [p[axis[0]], p[axis[1]], p[axis[2]]];53 for k in 0..3 {54 if signs >> k & 1 == 1 {55 q[k] = n - 1 - q[k];56 }57 }58 (q[0] * n + q[1]) * n + q[2]59 })60 .collect();61 if !out.contains(&map) {62 out.push(map);63 }64 }65 }66 out67}6869pub fn classes(cells: usize, group: &[Vec<usize>]) -> Classes {70 let mut index = vec![u32::MAX; cells * cells];71 let mut count = 0u32;72 for j in 0..cells {73 for l in j..cells {74 if index[j * cells + l] != u32::MAX {75 continue;76 }77 for map in group {78 let (a, b) = (map[j], map[l]);79 index[a * cells + b] = count;80 index[b * cells + a] = count;81 }82 count += 1;83 }84 }85 Classes {86 index,87 count: count as usize,88 cells,89 }90}9192pub fn burnside(cells: usize, group: &[Vec<usize>]) -> u128 {93 let mut total = 0u128;94 for map in group {95 let mut seen = vec![false; cells];96 let mut cycles = 0u32;97 for j in 0..cells {98 if seen[j] {99 continue;100 }101 cycles += 1;102 let mut x = j;103 while !seen[x] {104 seen[x] = true;105 x = map[x];106 }107 }108 total += 1u128 << cycles;109 }110 total / group.len() as u128111}112113pub fn census(points: &[usize], table: &Classes) -> Vec<u32> {114 let mut out = vec![0u32; table.count];115 for a in 0..points.len() {116 for b in a..points.len() {117 out[table.index[points[a] * table.cells + points[b]] as usize] += 1;118 }119 }120 out121}122123pub fn canon(bits: u128, cells: usize, group: &[Vec<usize>]) -> u128 {124 let mut best = u128::MAX;125 for map in group {126 let mut x = 0u128;127 for j in 0..cells {128 if bits >> j & 1 == 1 {129 x |= 1u128 << map[j];130 }131 }132 best = best.min(x);133 }134 best135}136137fn cells_of(bits: u128, cells: usize) -> Vec<usize> {138 (0..cells).filter(|j| bits >> j & 1 == 1).collect()139}140141fn kron_plane(bits: u128, q: usize) -> Vec<usize> {142 let base = cells_of(bits, q * q);143 let side = q * q;144 let mut out = Vec::with_capacity(base.len() * base.len());145 for &p in &base {146 for &s in &base {147 let r = (p / q) * q + s / q;148 let c = (p % q) * q + s % q;149 out.push(r * side + c);150 }151 }152 out.sort_unstable();153 out154}155156fn picture(bits: u128) -> String {157 (0..3)158 .map(|r| {159 (0..3)160 .map(|c| {161 if bits >> (r * 3 + c) & 1 == 1 {162 '#'163 } else {164 '.'165 }166 })167 .collect::<String>()168 })169 .collect::<Vec<_>>()170 .join("/")171}172173fn solve(matrix: &mut Vec<Vec<f64>>, rhs: &mut Vec<f64>) -> Option<Vec<f64>> {174 let n = rhs.len();175 for col in 0..n {176 let mut pivot = col;177 for row in col..n {178 if matrix[row][col].abs() > matrix[pivot][col].abs() {179 pivot = row;180 }181 }182 if matrix[pivot][col].abs() < 1e-9 {183 return None;184 }185 matrix.swap(col, pivot);186 rhs.swap(col, pivot);187 for row in 0..n {188 if row == col {189 continue;190 }191 let f = matrix[row][col] / matrix[col][col];192 for k in col..n {193 matrix[row][k] -= f * matrix[col][k];194 }195 rhs[row] -= f * rhs[col];196 }197 }198 Some((0..n).map(|i| rhs[i] / matrix[i][i]).collect())199}200201fn rank(rows: &[Vec<f64>], tolerance: f64) -> usize {202 let mut work: Vec<Vec<f64>> = rows.to_vec();203 let width = work[0].len();204 let mut got = 0usize;205 for col in 0..width {206 let mut pivot = None;207 for row in got..work.len() {208 if work[row][col].abs() > tolerance {209 pivot = Some(row);210 break;211 }212 }213 let Some(pivot) = pivot else { continue };214 work.swap(got, pivot);215 for row in 0..work.len() {216 if row == got {217 continue;218 }219 let f = work[row][col] / work[got][col];220 for k in col..width {221 work[row][k] -= f * work[got][k];222 }223 }224 got += 1;225 if got == work.len() {226 break;227 }228 }229 got230}231232pub fn control() {233 println!();234 println!("SHAPE READING, THE 511 AS CONTROL");235 let g1 = d4(BASE);236 let g2 = d4(BASE * BASE);237 let t1 = classes(9, &g1);238 let t2 = classes(81, &g2);239 println!(240 " pair classes level 1 {} level 2 {}; Burnside orbits {} of which {} nonempty",241 t1.count,242 t2.count,243 burnside(9, &g1),244 burnside(9, &g1) - 1245 );246 let mut mismatch = 0usize;247 for code in 1..512u128 {248 let grid = plane(code, BASE, 2);249 let live: Vec<usize> = grid250 .bytes()251 .iter()252 .enumerate()253 .filter(|(_, b)| **b != 0)254 .map(|(f, _)| f)255 .collect();256 if live != kron_plane(code, BASE) {257 mismatch += 1;258 }259 }260 println!(" level-2 renders disagreeing with the Kronecker square: {mismatch}");261 let mut reps: Vec<u128> = (1..512u128).map(|c| canon(c, 9, &g1)).collect();262 reps.sort_unstable();263 reps.dedup();264 println!(" nonempty orbits by canonical form: {}", reps.len());265 let mut first: HashMap<Vec<u32>, Vec<u128>> = HashMap::new();266 let mut second: HashMap<Vec<u32>, Vec<u128>> = HashMap::new();267 for &code in &reps {268 first269 .entry(census(&cells_of(code, 9), &t1))270 .or_default()271 .push(code);272 second273 .entry(census(&kron_plane(code, BASE), &t2))274 .or_default()275 .push(code);276 }277 println!(278 " distinct pair censuses: level 1 {} level 2 {}",279 first.len(),280 second.len()281 );282 let mut clashes: Vec<Vec<u128>> = first.values().filter(|v| v.len() > 1).cloned().collect();283 clashes.sort();284 let spin1 = orbit::census(1, 1024, 12);285 let spin2 = orbit::census(2, 768, 12);286 for group in &clashes {287 let (a, b) = (group[0], group[1]);288 println!(289 " level-1 clash {a} {} against {b} {} fill {} spectrum gap level 1 {:.2e} level 2 {:.2e}",290 picture(a),291 picture(b),292 a.count_ones(),293 orbit::gap(&spin1[a as usize], &spin1[b as usize]),294 orbit::gap(&spin2[a as usize], &spin2[b as usize])295 );296 }297 let mut buckets: Vec<usize> = Vec::new();298 for code in 1..512usize {299 if !buckets300 .iter()301 .any(|&head| orbit::agree(&spin1[head], &spin1[code], 1e-9))302 {303 buckets.push(code);304 }305 }306 println!(307 " spectra from level 1 alone: {} against {} pair censuses",308 buckets.len(),309 first.len()310 );311 let rows: Vec<Vec<f64>> = (1..512usize)312 .map(|code| {313 census(&cells_of(code as u128, 9), &t1)314 .iter()315 .map(|&v| v as f64)316 .collect()317 })318 .collect();319 let mut chosen: Vec<usize> = Vec::new();320 for i in 0..rows.len() {321 let mut trial: Vec<Vec<f64>> = chosen.iter().map(|&j| rows[j].clone()).collect();322 trial.push(rows[i].clone());323 if rank(&trial, 1e-9) == trial.len() {324 chosen.push(i);325 }326 if chosen.len() == t1.count {327 break;328 }329 }330 let mut weights: Vec<Vec<f64>> = Vec::new();331 let mut worst = 0.0f64;332 for m in 0..13 {333 let mut matrix: Vec<Vec<f64>> = chosen.iter().map(|&j| rows[j].clone()).collect();334 let mut rhs: Vec<f64> = chosen.iter().map(|&j| spin1[j + 1][m]).collect();335 let Some(w) = solve(&mut matrix, &mut rhs) else {336 println!(" the level-1 system is singular at order {m}");337 return;338 };339 let scale = (1..512usize)340 .map(|c| spin1[c][m].abs())341 .fold(0.0f64, f64::max)342 .max(1e-300);343 for (i, row) in rows.iter().enumerate() {344 let predicted: f64 = row.iter().zip(&w).map(|(a, b)| a * b).sum();345 worst = worst.max((predicted - spin1[i + 1][m]).abs() / scale);346 }347 let top = w.iter().fold(0.0f64, |a, b| a.max(b.abs())).max(1e-300);348 weights.push(w.iter().map(|v| v / top).collect());349 }350 let odd: Vec<Vec<f64>> = weights351 .iter()352 .enumerate()353 .filter(|(m, _)| m % 2 == 1)354 .map(|(_, w)| w.clone())355 .collect();356 let even: Vec<Vec<f64>> = weights357 .iter()358 .enumerate()359 .filter(|(m, _)| m % 2 == 0)360 .map(|(_, w)| w.clone())361 .collect();362 println!(363 " the level-1 spectrum is a linear functional of the pair census: independent censuses solved {}, worst relative residual over 511 codes by 13 orders {:.2e}, coefficient rank {} of {}",364 chosen.len(),365 worst,366 rank(&weights, 1e-8),367 t1.count368 );369 println!(370 " half-turning one member of a pair caps the {} odd orders at the 3-dimensional antisymmetric part and they reach {}; the {} even orders reach {}",371 odd.len(),372 rank(&odd, 1e-8),373 even.len(),374 rank(&even, 1e-8)375 );376}377378fn chunk_tables(cells: usize, group: &[Vec<usize>]) -> Vec<Vec<Vec<u32>>> {379 let chunks = cells.div_ceil(8);380 group381 .iter()382 .map(|map| {383 (0..chunks)384 .map(|c| {385 (0..256u32)386 .map(|v| {387 let mut x = 0u32;388 for b in 0..8 {389 let j = c * 8 + b;390 if j < cells && v >> b & 1 == 1 {391 x |= 1u32 << map[j];392 }393 }394 x395 })396 .collect()397 })398 .collect()399 })400 .collect()401}402403fn image(code: u32, tabs: &[Vec<u32>]) -> u32 {404 let mut x = 0u32;405 for (c, tab) in tabs.iter().enumerate() {406 x |= tab[(code >> (8 * c) & 255) as usize];407 }408 x409}410411fn pack(counts: &[u32], width: u32) -> (u128, u128) {412 let per = 128 / width as usize;413 let (mut lo, mut hi) = (0u128, 0u128);414 for (c, &v) in counts.iter().enumerate() {415 assert!(v < 1 << width, "a class count overflows the key");416 if c < per {417 lo |= (v as u128) << (c as u32 * width);418 } else {419 hi |= (v as u128) << ((c - per) as u32 * width);420 }421 }422 (lo, hi)423}424425fn pair_list(table: &Classes) -> Vec<(u32, u16)> {426 let mut out = Vec::new();427 for j in 0..table.cells {428 for l in j + 1..table.cells {429 out.push((430 (1u32 << j) | (1u32 << l),431 table.index[j * table.cells + l] as u16,432 ));433 }434 }435 out436}437438fn level_one_sweep(439 cells: usize,440 group: &[Vec<usize>],441 table: &Classes,442 width: u32,443) -> Vec<(u128, u128, u32)> {444 let tabs = chunk_tables(cells, group);445 let pairs = pair_list(table);446 let diag: Vec<u16> = (0..cells)447 .map(|j| table.index[j * cells + j] as u16)448 .collect();449 let mut out: Vec<(u128, u128, u32)> = Vec::new();450 let mut counts = vec![0u32; table.count];451 for code in 1..(1u32 << cells) {452 if tabs.iter().any(|t| image(code, t) < code) {453 continue;454 }455 counts.iter_mut().for_each(|v| *v = 0);456 for j in 0..cells {457 if code >> j & 1 == 1 {458 counts[diag[j] as usize] += 1;459 }460 }461 for &(mask, class) in &pairs {462 if code & mask == mask {463 counts[class as usize] += 1;464 }465 }466 let (lo, hi) = pack(&counts, width);467 out.push((lo, hi, code));468 }469 out.sort_unstable();470 out471}472473fn groups(sweep: &[(u128, u128, u32)]) -> Vec<Vec<u32>> {474 let mut out: Vec<Vec<u32>> = Vec::new();475 let mut k = 0usize;476 while k < sweep.len() {477 let mut j = k + 1;478 while j < sweep.len() && sweep[j].0 == sweep[k].0 && sweep[j].1 == sweep[k].1 {479 j += 1;480 }481 if j - k > 1 {482 out.push(sweep[k..j].iter().map(|e| e.2).collect());483 }484 k = j;485 }486 out487}488489fn mix(x: u64) -> u64 {490 let mut z = x.wrapping_add(0x9E3779B97F4A7C15);491 z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);492 z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);493 z ^ (z >> 31)494}495496fn digest(points: &[usize], table: &Classes, scratch: &mut [u32], touched: &mut Vec<u32>) -> u128 {497 touched.clear();498 for a in 0..points.len() {499 for b in a..points.len() {500 let c = table.index[points[a] * table.cells + points[b]] as usize;501 if scratch[c] == 0 {502 touched.push(c as u32);503 }504 scratch[c] += 1;505 }506 }507 let (mut one, mut two) = (0u64, 0u64);508 for &c in touched.iter() {509 let n = scratch[c as usize];510 scratch[c as usize] = 0;511 let k = mix(((c as u64) << 34) ^ n as u64);512 one = one.wrapping_add(k);513 two ^= mix(k ^ 0x51_7C_C1_B7_27_22_0A_95);514 }515 (one as u128) << 64 | two as u128516}517518fn exact(points: &[usize], table: &Classes) -> Vec<(u32, u32)> {519 let mut seen: Vec<(u32, u32)> = Vec::new();520 let mut tally = vec![0u32; table.count];521 for a in 0..points.len() {522 for b in a..points.len() {523 tally[table.index[points[a] * table.cells + points[b]] as usize] += 1;524 }525 }526 for (c, &n) in tally.iter().enumerate() {527 if n > 0 {528 seen.push((c as u32, n));529 }530 }531 seen532}533534fn kron_cube(bits: u32, q: usize) -> Vec<usize> {535 let base: Vec<usize> = (0..q * q * q).filter(|j| bits >> j & 1 == 1).collect();536 let side = q * q;537 let mut out = Vec::with_capacity(base.len() * base.len());538 for &p in &base {539 for &s in &base {540 let a = (p / (q * q)) * q + s / (q * q);541 let b = ((p / q) % q) * q + (s / q) % q;542 let c = (p % q) * q + s % q;543 out.push((a * side + b) * side + c);544 }545 }546 out.sort_unstable();547 out548}549550fn level_two_window(551 clashes: &[Vec<u32>],552 table: &Classes,553 render: impl Fn(u32) -> Vec<usize>,554 budget: u128,555) -> (u32, usize, usize, Vec<(u32, u32)>) {556 let mut strata: Vec<Vec<&Vec<u32>>> = vec![Vec::new(); 64];557 for group in clashes {558 let weight = group.iter().map(|c| c.count_ones()).max().unwrap_or(0);559 strata[weight as usize].push(group);560 }561 let (mut top, mut done, mut spent) = (0u32, 0usize, 0u128);562 let mut survivors: Vec<(u32, u32)> = Vec::new();563 let mut scratch = vec![0u32; table.count];564 let mut touched: Vec<u32> = Vec::new();565 for weight in 0..strata.len() {566 let stratum = &strata[weight];567 if stratum.is_empty() {568 continue;569 }570 let bill: u128 = stratum571 .iter()572 .map(|group| {573 group574 .iter()575 .map(|c| {576 let n = c.count_ones() as u128 * c.count_ones() as u128;577 n * (n + 1) / 2578 })579 .sum::<u128>()580 })581 .sum();582 if spent + bill > budget {583 break;584 }585 spent += bill;586 top = weight as u32;587 done += stratum.len();588 for group in stratum {589 let seen: Vec<u128> = group590 .iter()591 .map(|&c| digest(&render(c), table, &mut scratch, &mut touched))592 .collect();593 for a in 0..group.len() {594 for b in a + 1..group.len() {595 if seen[a] == seen[b]596 && exact(&render(group[a]), table) == exact(&render(group[b]), table)597 {598 survivors.push((group[a], group[b]));599 }600 }601 }602 }603 }604 (top, done, clashes.len(), survivors)605}606607pub fn base_five() {608 println!();609 println!("SHAPE READING, BASE FIVE PLANE");610 let g1 = d4(5);611 let g2 = d4(25);612 let t1 = classes(25, &g1);613 let t2 = classes(625, &g2);614 let orbits = burnside(25, &g1);615 println!(616 " pair classes level 1 {} level 2 {}; Burnside orbits {} of which {} nonempty",617 t1.count,618 t2.count,619 orbits,620 orbits - 1621 );622 let sweep = level_one_sweep(25, &g1, &t1, 4);623 println!(624 " canonical codes swept over all 2^25: {}; distinct level-1 pair censuses: {}",625 sweep.len(),626 {627 let mut k = 0usize;628 for j in 0..sweep.len() {629 if j == 0 || (sweep[j].0, sweep[j].1) != (sweep[j - 1].0, sweep[j - 1].1) {630 k += 1;631 }632 }633 k634 }635 );636 let clashes = groups(&sweep);637 let orbits_in: usize = clashes.iter().map(|g| g.len()).sum();638 println!(639 " level-1 collisions: {} censuses shared by {} orbits, the heaviest group holding {}",640 clashes.len(),641 orbits_in,642 clashes.iter().map(|g| g.len()).max().unwrap_or(0)643 );644 let (top, done, all, alive) =645 level_two_window(&clashes, &t2, |c| kron_plane(c as u128, 5), 24_000_000_000);646 println!(647 " level-2 window: every colliding census of weight at most {top}, {done} of {all}; pairs surviving the level-2 pair census: {}",648 alive.len()649 );650 for (a, b) in alive.iter().take(4) {651 println!(" surviving pair {a} against {b}");652 }653}654655pub fn cube_three() {656 println!();657 println!("SHAPE READING, BASE THREE AT D = 3");658 let g1 = oct(3);659 let g2 = oct(9);660 let t1 = classes(27, &g1);661 let t2 = classes(729, &g2);662 let orbits = burnside(27, &g1);663 println!(664 " cube group order {}; pair classes level 1 {} level 2 {}; Burnside orbits {} of which {} nonempty",665 g1.len(),666 t1.count,667 t2.count,668 orbits,669 orbits - 1670 );671 let sweep = level_one_sweep(27, &g1, &t1, 6);672 println!(673 " canonical codes swept over all 2^27: {}; distinct level-1 pair censuses: {}",674 sweep.len(),675 {676 let mut k = 0usize;677 for j in 0..sweep.len() {678 if j == 0 || (sweep[j].0, sweep[j].1) != (sweep[j - 1].0, sweep[j - 1].1) {679 k += 1;680 }681 }682 k683 }684 );685 let clashes = groups(&sweep);686 let orbits_in: usize = clashes.iter().map(|g| g.len()).sum();687 println!(688 " level-1 collisions: {} censuses shared by {} orbits, the heaviest group holding {}",689 clashes.len(),690 orbits_in,691 clashes.iter().map(|g| g.len()).max().unwrap_or(0)692 );693 let (top, done, all, alive) =694 level_two_window(&clashes, &t2, |c| kron_cube(c, 3), 200_000_000_000);695 println!(696 " level-2 window: every colliding census of weight at most {top}, {done} of {all}; pairs surviving the level-2 pair census: {}",697 alive.len()698 );699 for (a, b) in alive.iter().take(4) {700 println!(" surviving pair {a} against {b}");701 }702}