main.rs
17.4 kB · rust · 525 lines
1mod design;2mod gaussian;3mod mass;4mod orbit;5mod powder;6mod shadow;7mod shape;89use design::{bit_cells, canonical, carry, plane, square_group, BASE};10use mass::{distance, horizon, ripple, scaling_error, shells, Ripple};1112const SUBJECTS: [u128; 7] = [79, 95, 127, 239, 255, 495, 511];13const BINS: usize = 24;14const LOW: f64 = 27.0;1516fn dimension(code: u128) -> f64 {17 (code.count_ones() as f64).ln() / (BASE as f64).ln()18}1920fn list(values: &[f64], places: usize) -> String {21 values22 .iter()23 .map(|value| format!("{value:.places$}"))24 .collect::<Vec<_>>()25 .join(" ")26}2728fn corner_bit(table: &[(usize, usize)]) -> usize {29 table30 .iter()31 .position(|cell| *cell == (0, 0))32 .expect("a corner digit")33}3435fn centre_bit(table: &[(usize, usize)]) -> usize {36 table37 .iter()38 .position(|cell| *cell == (1, 1))39 .expect("a centre digit")40}4142fn read(43 code: u128,44 level: usize,45 digit: (usize, usize),46 table: &[(usize, usize)],47) -> (Ripple, f64, u64) {48 let grid = plane(code, BASE, level);49 let side = grid.shape[0] as f64;50 let bulk = shells(&grid, digit);51 let far = horizon(code, digit, table) * side;52 let error = scaling_error(&bulk, LOW, far, code.count_ones() as f64);53 let curve = ripple(&bulk, LOW, far, dimension(code), BINS);54 (curve, error, bulk.total())55}5657fn anchors(table: &[(usize, usize)]) {58 println!("ANCHORS");59 println!("bit to cell map: {table:?}");60 let group = square_group();61 let mut seen = vec![false; 512];62 let mut classes = 0;63 for code in 0..512u128 {64 if seen[code as usize] {65 continue;66 }67 classes += 1;68 for member in design::orbit(&group, code, table) {69 seen[member as usize] = true;70 }71 }72 println!(73 "square group order {} classes over 512 codes {classes}",74 group.len()75 );76 let grid = plane(127, BASE, 4);77 let bulk = shells(&grid, (1, 1));78 let side = grid.shape[0];79 let rings = mrlynum::spin::profile(&design::floats(&grid), side, 6000);80 let crate_mass = mrlynum::spin::mass(&rings, side);81 println!(82 "code 127 level 4: shell total {} crate profile mass {:.2} gap {:.2e}",83 bulk.total(),84 crate_mass,85 (crate_mass - bulk.total() as f64).abs() / bulk.total() as f6486 );87 let step = mrlynum::spin::reach(side) / 5999.0;88 for radius in [9.0f64, 27.0, 40.5] {89 let cut = (radius / step) as usize;90 let partial: f64 = (1..=cut)91 .map(|k| {92 let (a, b) = (k as f64 * step, (k - 1) as f64 * step);93 std::f64::consts::PI * (a + b) * (rings[k] as f64 + rings[k - 1] as f64) / 2.094 * step95 })96 .sum();97 println!(98 " radius {radius:>5.1} shell count {:>6} crate profile integral {:>9.2} ratio {:.4}",99 bulk.at(radius),100 partial,101 partial / bulk.at(radius) as f64102 );103 }104}105106fn holes(table: &[(usize, usize)]) {107 println!();108 println!("THE CENTRE HOLE, EXACT NEAREST FILLED CELL TO THE RASTER CENTRE");109 let centre = centre_bit(table);110 for code in SUBJECTS {111 let grid = plane(code, BASE, 5);112 let side = grid.shape[0];113 let four = mass::nearest_cell(&grid, (1, 1));114 let sixth = (BASE as u64).pow(4).pow(2);115 println!(116 " code {code:>3} fill {} centre digit {} four times the squared distance {four} against (side/3)^2 = {sixth} at least {} distance {:.6} side/6 {:.6}",117 code.count_ones(),118 if code >> centre & 1 == 1 { "filled" } else { "empty " },119 four >= sixth,120 (four as f64).sqrt() / 2.0,121 side as f64 / 6.0122 );123 }124}125126fn spin_dimension(table: &[(usize, usize)], level: usize) {127 println!();128 println!(129 "SPIN DIMENSION AT LEVEL {level}, FIXED POINT OF THE CORNER DIGIT AND OF THE CENTRE DIGIT"130 );131 let corner = corner_bit(table);132 let centre = centre_bit(table);133 for code in SUBJECTS {134 for (name, bit, digit) in [135 ("corner", corner, (0usize, 0usize)),136 ("centre", centre, (1, 1)),137 ] {138 if code >> bit & 1 == 0 {139 println!(" code {code:>3} {name}: digit empty, no fixed point");140 continue;141 }142 let (curve, error, total) = read(code, level, digit, table);143 println!(144 " code {code:>3} fill {} {name}: D exact {:.6} slope {:.6} gap {:.2e} periods {} scaling error {:.2e} ripple swing {:.5} drift {:.5} mass {total}",145 code.count_ones(),146 dimension(code),147 curve.slope,148 (curve.slope - dimension(code)).abs(),149 curve.periods,150 error,151 curve.swing,152 curve.drift153 );154 }155 }156}157158fn acid(table: &[(usize, usize)], level: usize) {159 println!();160 println!("THE EQUAL DIMENSION PAIRS AT LEVEL {level}, CORNER FIXED POINT");161 let corner = corner_bit(table);162 let mut kept: Vec<(u128, Ripple)> = Vec::new();163 for code in SUBJECTS {164 if code >> corner & 1 == 0 {165 continue;166 }167 let (curve, _, _) = read(code, level, (0, 0), table);168 println!(" code {code:>3} ripple {}", list(&curve.curve, 4));169 kept.push((code, curve));170 }171 for left in 0..kept.len() {172 for right in left + 1..kept.len() {173 if kept[left].0.count_ones() != kept[right].0.count_ones() {174 continue;175 }176 let gap = distance(&kept[left].1.curve, &kept[right].1.curve);177 let bar = kept[left].1.drift.max(kept[right].1.drift);178 println!(179 " {} against {}: fill {} ripple gap {:.5} drift bar {:.5} ratio {:.1}",180 kept[left].0,181 kept[right].0,182 kept[left].0.count_ones(),183 gap,184 bar,185 gap / bar186 );187 }188 }189}190191fn ripple_census(table: &[(usize, usize)], level: usize) {192 println!();193 println!("RIPPLE CENSUS AT LEVEL {level}: EVERY CODE WITH A FILLED CORNER DIGIT");194 let corner = corner_bit(table);195 let group = square_group();196 let transpose = group197 .iter()198 .find(|map| map[1] == BASE && map[BASE] == 1)199 .expect("the transpose")200 .clone();201 let flip: Vec<u128> = (0..512u128)202 .map(|code| carry(&transpose, code, table))203 .collect();204 let mut curves: Vec<(u128, Ripple)> = Vec::new();205 for code in 1..512u128 {206 if code >> corner & 1 == 0 {207 continue;208 }209 let (curve, _, _) = read(code, level, (0, 0), table);210 curves.push((code, curve));211 }212 let index = |code: u128| curves.iter().position(|row| row.0 == code);213 let mut mirror: f64 = 0.0;214 for (code, curve) in &curves {215 if let Some(at) = index(flip[*code as usize]) {216 mirror = mirror.max(distance(&curve.curve, &curves[at].1.curve));217 }218 }219 println!(220 " codes read {} transpose control, worst ripple gap {:.2e}",221 curves.len(),222 mirror223 );224 let stamp = |code: u128| code.min(flip[code as usize]);225 let mut collisions: Vec<(f64, f64, u128, u128)> = Vec::new();226 let mut closest = (f64::INFINITY, 0u128, 0u128);227 for left in 0..curves.len() {228 for right in left + 1..curves.len() {229 let (a, b) = (curves[left].0, curves[right].0);230 if a.count_ones() != b.count_ones() || stamp(a) == stamp(b) {231 continue;232 }233 if a != stamp(a) || b != stamp(b) {234 continue;235 }236 let gap = distance(&curves[left].1.curve, &curves[right].1.curve);237 let bar = curves[left].1.drift.max(curves[right].1.drift);238 if gap < closest.0 {239 closest = (gap, a, b);240 }241 if gap < bar {242 collisions.push((gap, bar, a, b));243 }244 }245 }246 collisions.sort_by(|x, y| x.0.partial_cmp(&y.0).expect("finite"));247 println!(248 " transpose class pairs of equal fill whose ripples sit inside their own drift bar: {}",249 collisions.len()250 );251 for (gap, bar, left, right) in &collisions {252 println!(253 " classes {left:>3} and {right:>3} fill {} gap {gap:.5} bar {bar:.5} ratio {:.2}",254 left.count_ones(),255 gap / bar256 );257 }258 let swing = |code: u128| curves[index(code).expect("a read code")].1.swing;259 println!(260 " closest distinct class pair: {} and {} fill {} gap {:.5} with ripple swings {:.5} and {:.5}",261 closest.1,262 closest.2,263 closest.1.count_ones(),264 closest.0,265 swing(closest.1),266 swing(closest.2)267 );268}269270fn powder_rings(level: usize, pad: usize, only: &[u128]) {271 println!();272 println!(273 "POWDER RINGS AT LEVEL {level}, PAD {pad}, RING AVERAGED POWER AGAINST THE FREQUENCY INDEX"274 );275 let side = (BASE as f64).powi(level as i32);276 let low = 3.0 * pad as f64 / side;277 let high = pad as f64 / 8.0;278 println!(" band {low:.1} to {high:.1} in frequency index, bins 240, phases {BINS}, slide window 3 periods by a quarter period");279 for code in SUBJECTS {280 if !only.is_empty() && !only.contains(&code) {281 continue;282 }283 let grid = plane(code, BASE, level);284 let read = powder::powder(&grid, pad, low, high, 240, BINS);285 let dim = dimension(code);286 println!(287 " code {code:>3} fill {} D {:.6} band slope {:.5} against -D {:.5} gap {:.4} slide {:.5} to {:.5} spread {:.4} porod -3 gap {:.4} log period swing {:.4}",288 code.count_ones(),289 dim,290 read.slope,291 -dim,292 (read.slope + dim).abs(),293 read.low,294 read.high,295 read.high - read.low,296 (read.slope + 3.0).abs(),297 read.swing298 );299 }300}301302fn spin_spectrum(table: &[(usize, usize)]) {303 println!();304 println!("SPIN SPECTRUM, P_m OVER ALL 512 BASE THREE CODES");305 let group = square_group();306 let stamp: Vec<u128> = (0..512u128)307 .map(|code| canonical(&group, code, table))308 .collect();309 let first = orbit::census(1, 1024, 12);310 let second = orbit::census(2, 768, 12);311 let mut pairs: Vec<(u128, u128, f64, f64)> = Vec::new();312 for left in 1..512u128 {313 for right in left + 1..512u128 {314 if left.count_ones() != right.count_ones() {315 continue;316 }317 if stamp[left as usize] == stamp[right as usize] {318 continue;319 }320 let one = orbit::gap(&first[left as usize], &first[right as usize]);321 let two = orbit::gap(&second[left as usize], &second[right as usize]);322 if one < 1e-9 && two < 1e-9 {323 pairs.push((left, right, one, two));324 }325 }326 }327 println!(328 " pairs outside one square class agreeing at levels 1 and 2 to 1e-9: {}",329 pairs.len()330 );331 let mut classes: Vec<(u128, u128)> = pairs332 .iter()333 .map(|(a, b, _, _)| (stamp[*a as usize], stamp[*b as usize]))334 .collect();335 classes.sort_unstable();336 classes.dedup();337 println!(" distinct class pairs among them: {}", classes.len());338 for (left, right) in classes.iter().take(8) {339 let third = orbit::gap(340 &orbit::spectrum(*left, 3, 1024, 24),341 &orbit::spectrum(*right, 3, 1024, 24),342 );343 println!(344 " classes {left} and {right} fill {} level 3 gap {:.2e} verdict {}",345 left.count_ones(),346 third,347 if third < 1e-9 {348 "isospectral"349 } else {350 "separated at level 3"351 }352 );353 }354 let mut buckets: Vec<Vec<u128>> = Vec::new();355 for code in 1..512u128 {356 let mut placed = false;357 for bucket in buckets.iter_mut() {358 let head = bucket[0] as usize;359 if first[head].len() == first[code as usize].len()360 && orbit::agree(&first[head], &first[code as usize], 1e-9)361 && orbit::agree(&second[head], &second[code as usize], 1e-9)362 {363 bucket.push(code);364 placed = true;365 break;366 }367 }368 if !placed {369 buckets.push(vec![code]);370 }371 }372 println!(373 " distinct spin spectra over the 511 nonempty codes at levels 1 and 2 together: {}",374 buckets.len()375 );376 let big = buckets.iter().map(|bucket| bucket.len()).max().unwrap_or(0);377 println!(" largest spectral bucket holds {big} codes");378}379380fn sponge_shadow(top: usize) {381 println!();382 println!("THE SPONGE SHADOW: LATTICE LINES IN DIRECTION (a,b,c) MEETING THE LEVEL L SPONGE, AND THE SAME FOR THE SOLID CUBE");383 let digits = shadow::digits();384 let cube = shadow::cube_digits();385 println!(386 " sponge digits {} cube digits {}",387 digits.len(),388 cube.len()389 );390 let mut full = Vec::new();391 for view in shadow::views(3) {392 let counts: Vec<usize> = (1..=top)393 .map(|level| shadow::shadow(level, view, &digits))394 .collect();395 let solid: Vec<usize> = (1..=top - 1)396 .map(|level| shadow::shadow(level, view, &cube))397 .collect();398 let share: Vec<f64> = solid399 .iter()400 .zip(&counts)401 .map(|(all, hit)| *hit as f64 / *all as f64)402 .collect();403 let opaque = share.iter().all(|value| (value - 1.0).abs() < 1e-12);404 if opaque {405 full.push(view);406 }407 println!(408 " view {view:?} sponge {counts:?} cube {solid:?} share {} {}",409 list(&share, 5),410 if opaque { "opaque" } else { "see through" }411 );412 }413 println!(414 " views whose lattice lines the sponge blocks completely to level {}: {full:?}",415 top - 1416 );417 let deep = shadow::shadow(top, [1, 1, 1], &cube);418 println!(419 " the space diagonal at level {top}: sponge {} cube {deep} equal {}",420 shadow::shadow(top, [1, 1, 1], &digits),421 shadow::shadow(top, [1, 1, 1], &digits) == deep422 );423}424425fn gaussian_farey(top: usize) {426 println!();427 println!("THE GAUSSIAN FAREY: RADII NEW AT SCALE n");428 let cap = 2 * top * top;429 let least = gaussian::least_factors(cap);430 let direct = gaussian::union_counts(top, false);431 let boxed = gaussian::union_counts(top, true);432 let rule: Vec<usize> = (1..=top)433 .map(|scale| gaussian::new_disc(scale, &least))434 .collect();435 let mismatch = direct436 .iter()437 .zip(&rule)438 .enumerate()439 .filter(|(_, (a, b))| a != b)440 .map(|(index, _)| index + 1)441 .collect::<Vec<_>>();442 println!(443 " disc reading, new radii at n = 1..20: {:?}",444 &direct[..20.min(direct.len())]445 );446 println!(447 " the square free rule reproduces the disc reading at every n up to {top}: {}",448 mismatch.is_empty()449 );450 if !mismatch.is_empty() {451 println!(452 " first mismatches at n = {:?}",453 &mismatch[..8.min(mismatch.len())]454 );455 }456 println!(457 " box reading, new radii at n = 1..20: {:?}",458 &boxed[..20.min(boxed.len())]459 );460 let split: Vec<usize> = (1..=top)461 .filter(|n| direct[n - 1] != boxed[n - 1])462 .collect();463 println!(464 " first scale where the box and the disc disagree: {:?}",465 split.first()466 );467 let primitives: Vec<usize> = (1..=20)468 .map(|scale| {469 (1..=2 * scale * scale)470 .filter(|norm| gaussian::primitive_norm(*norm, &least))471 .count()472 })473 .collect();474 println!(" norms below 2n^2 with a primitive representation, n = 1..20: {primitives:?}");475 let prefix = gaussian::two_square_prefix(cap, &least);476 let sieve: Vec<usize> = (1..=top)477 .map(|scale| gaussian::mobius_count(scale, &prefix, &least))478 .collect();479 let broken: Vec<usize> = (1..=top)480 .filter(|n| sieve[n - 1] != direct[n - 1])481 .collect();482 println!(" the Mobius identity new(n) = sum_d mu(d) B(2n^2/d^2) over d | rad(n) holds to n = {top}: {}", broken.is_empty());483 if !broken.is_empty() {484 println!(485 " it first fails at n = {:?}",486 &broken[..4.min(broken.len())]487 );488 }489 let reach = 192usize;490 let wide = gaussian::least_factors(2 * reach * reach);491 let long = gaussian::two_square_prefix(2 * reach * reach, &wide);492 println!(493 " the radical six family, new(n)/B(2n^2) climbing to the Jordan factor {:.5}:",494 gaussian::jordan(6, &wide)495 );496 for scale in [6usize, 12, 24, 48, 96, 192] {497 let count = gaussian::mobius_count(scale, &long, &wide);498 let all = long[2 * scale * scale];499 println!(500 " n {scale:>3} new {count:>6} all {all:>6} ratio {:.5}",501 count as f64 / all as f64502 );503 }504}505506fn main() {507 let table = bit_cells();508 if std::env::args().any(|a| a == "shape") {509 shape::control();510 shape::base_five();511 shape::cube_three();512 return;513 }514 anchors(&table);515 holes(&table);516 spin_dimension(&table, 6);517 acid(&table, 7);518 ripple_census(&table, 6);519 ripple_census(&table, 7);520 powder_rings(7, 4096, &[]);521 powder_rings(7, 8192, &[127, 255, 495]);522 spin_spectrum(&table);523 sponge_shadow(5);524 gaussian_farey(64);525}