main.rs
35.3 kB · rust · 1046 lines
1use mrlynum::factor::mobius_sieve;2use mrlynum::memory::{allowed_windows, kappa, perron, Rule};3use std::env;4use std::time::Instant;56// THE WINDOW PROFILE78fn induced_tables() -> (Vec<u8>, Vec<u8>) {9 let mut to2 = vec![0u8; 256];10 let mut to1 = vec![0u8; 256];11 for p in 0..256usize {12 let mut a = 0u8;13 let mut b = 0u8;14 for w in 0..8usize {15 if (p >> w) & 1 == 0 {16 continue;17 }18 a |= 1 << (w >> 1);19 a |= 1 << (w & 3);20 b |= 1 << ((w >> 2) & 1);21 b |= 1 << ((w >> 1) & 1);22 b |= 1 << (w & 1);23 }24 to2[p] = a;25 to1[p] = b;26 }27 (to2, to1)28}2930fn small_profiles(n: usize) -> (u8, u8, u8) {31 match n {32 1 => (0, 0, 2),33 2 => (0, 4, 3),34 3 => (0, 8, 2),35 _ => unreachable!(),36 }37}3839fn digit_profiles(n: usize) -> ([u8; 3], [u8; 3]) {40 let bits = usize::BITS as usize - n.leading_zeros() as usize;41 let digits: Vec<usize> = (0..bits).rev().map(|i| (n >> i) & 1).collect();42 let mut plain = [0u8; 3];43 let mut zeroed = [0u8; 3];44 let mut padded = vec![0usize];45 padded.extend_from_slice(&digits);46 for k in 1..=3usize {47 for run in digits.windows(k) {48 let mut w = 0usize;49 for &d in run {50 w = (w << 1) | d;51 }52 plain[k - 1] |= 1 << w;53 }54 for run in padded.windows(k) {55 let mut w = 0usize;56 for &d in run {57 w = (w << 1) | d;58 }59 zeroed[k - 1] |= 1 << w;60 }61 }62 (plain, zeroed)63}6465// THE ZETA TRANSFORM6667fn subset_sums_u64(source: &[u64], bits: u32) -> Vec<u64> {68 let mut out = source.to_vec();69 for b in 0..bits {70 for w in 0..out.len() {71 if (w >> b) & 1 == 1 {72 out[w] += out[w ^ (1 << b)];73 }74 }75 }76 out77}7879fn subset_sums_i64(source: &[i64], bits: u32) -> Vec<i64> {80 let mut out = source.to_vec();81 for b in 0..bits {82 for w in 0..out.len() {83 if (w >> b) & 1 == 1 {84 out[w] += out[w ^ (1 << b)];85 }86 }87 }88 out89}9091// THE PHASE GRID9293fn phase_grid(depth: u32, cap: usize) -> Vec<(String, usize)> {94 let steps = [95 1.0f64,96 1.189_207_115_002_721,97 1.414_213_562_373_095_1,98 1.681_792_830_507_429,99 ];100 let mut out: Vec<(String, usize)> = Vec::new();101 for level in 8..=depth {102 for (j, step) in steps.iter().enumerate() {103 let x = ((1u64 << level) as f64 * step).floor() as usize;104 if x > cap {105 continue;106 }107 let label = format!("{}.{:02}", level, j * 25);108 if out.last().map(|(_, v)| *v) == Some(x) {109 continue;110 }111 out.push((label, x));112 }113 }114 out115}116117// THE CLASS REPRESENTATIVES118119fn window_maps(width: usize) -> Vec<Vec<usize>> {120 let windows = 1usize << width;121 let flip: Vec<usize> = (0..windows).map(|w| windows - 1 - w).collect();122 let rev: Vec<usize> = (0..windows)123 .map(|w| {124 let mut out = 0usize;125 for j in 0..width {126 out |= ((w >> j) & 1) << (width - 1 - j);127 }128 out129 })130 .collect();131 let identity: Vec<usize> = (0..windows).collect();132 let mut maps = vec![identity.clone(), flip.clone(), rev.clone()];133 maps.push((0..windows).map(|w| flip[rev[w]]).collect());134 maps.sort();135 maps.dedup();136 maps137}138139fn representatives(width: usize) -> Vec<usize> {140 let windows = 1usize << width;141 let maps = window_maps(width);142 let codes = 1usize << windows;143 let mut rep = vec![usize::MAX; codes];144 for code in 0..codes {145 if rep[code] != usize::MAX {146 continue;147 }148 for map in &maps {149 let mut image = 0usize;150 for w in 0..windows {151 if (code >> w) & 1 == 1 {152 image |= 1 << map[w];153 }154 }155 rep[image] = code;156 }157 }158 rep159}160161// THE SWEEP162163struct Reading {164 mass: Vec<Vec<u64>>,165 meter: Vec<Vec<i64>>,166 peak: Vec<Vec<i64>>,167}168169struct Sweep {170 phases: Vec<(String, usize)>,171 read: [Reading; 3],172 mertens: Vec<i64>,173}174175fn sweep(mu: &[i8], profile: &[u8], phases: &[(String, usize)], anchors: &[usize]) -> Sweep {176 let (to2, to1) = induced_tables();177 let mut count = [vec![0u64; 4], vec![0u64; 16], vec![0u64; 256]];178 let mut mass = [vec![0i64; 4], vec![0i64; 16], vec![0i64; 256]];179 let mut current = [vec![0i64; 4], vec![0i64; 16], vec![0i64; 256]];180 let mut peak = [vec![0i64; 4], vec![0i64; 16], vec![0i64; 256]];181 let mut read = [182 Reading {183 mass: Vec::new(),184 meter: Vec::new(),185 peak: Vec::new(),186 },187 Reading {188 mass: Vec::new(),189 meter: Vec::new(),190 peak: Vec::new(),191 },192 Reading {193 mass: Vec::new(),194 meter: Vec::new(),195 peak: Vec::new(),196 },197 ];198 let mut mertens = Vec::new();199 let mut next_phase = 0usize;200 let mut next_anchor = 0usize;201 let cap = phases.last().map(|(_, x)| *x).unwrap_or(0);202 for n in 1..=cap {203 let (p1, p2, p3) = if n >= 4 {204 let p = profile[n] as usize;205 (to1[p] as usize, to2[p] as usize, p)206 } else {207 let (a, b, c) = small_profiles(n);208 (c as usize, b as usize, a as usize)209 };210 count[0][p1] += 1;211 count[1][p2] += 1;212 count[2][p3] += 1;213 let m = mu[n] as i64;214 if m != 0 {215 mass[0][p1] += m;216 mass[1][p2] += m;217 mass[2][p3] += m;218 for (slot, (seed, bits)) in [(p1, 3usize), (p2, 15), (p3, 255)].iter().enumerate() {219 let free = bits & !seed;220 let mut sub = free;221 loop {222 let rule = seed | sub;223 let value = current[slot][rule] + m;224 current[slot][rule] = value;225 let size = value.abs();226 if size > peak[slot][rule] {227 peak[slot][rule] = size;228 }229 if sub == 0 {230 break;231 }232 sub = (sub - 1) & free;233 }234 }235 }236 if next_anchor < anchors.len() && n == anchors[next_anchor] {237 mertens.push(current[0][3]);238 next_anchor += 1;239 }240 if next_phase < phases.len() && n == phases[next_phase].1 {241 for slot in 0..3 {242 let bits = [2u32, 4, 8][slot];243 let a = subset_sums_u64(&count[slot], bits);244 let z = subset_sums_i64(&mass[slot], bits);245 assert_eq!(z, current[slot], "zeta meter against the running meter");246 read[slot].mass.push(a);247 read[slot].meter.push(z);248 read[slot].peak.push(peak[slot].clone());249 }250 next_phase += 1;251 }252 }253 Sweep {254 phases: phases.to_vec(),255 read,256 mertens,257 }258}259260// THE MEMORYLESS CONTROLS261262fn design_meter(mu: &[i8], levels: usize, digits: &[u64]) -> (i64, i64, u64) {263 let mut powers = vec![1u64; levels + 1];264 for i in 1..=levels {265 powers[i] = powers[i - 1] * 3;266 }267 let lead: Vec<u64> = digits.iter().copied().filter(|&d| d != 0).collect();268 let base = digits.len();269 let mut meter = 0i64;270 let mut peak = 0i64;271 let mut mass = 0u64;272 for length in 1..=levels {273 let tails = base.pow((length - 1) as u32);274 for head in &lead {275 for tail in 0..tails {276 let mut value = head * powers[length - 1];277 let mut rest = tail;278 for place in 0..(length - 1) {279 value += digits[rest % base] * powers[place];280 rest /= base;281 }282 meter += mu[value as usize] as i64;283 if meter.abs() > peak {284 peak = meter.abs();285 }286 mass += 1;287 }288 }289 }290 (meter, peak, mass)291}292293// THE DIRECT RECOUNT294295fn recount(limit: usize) -> ([Vec<u64>; 3], [Vec<u64>; 3]) {296 let mut plain = [vec![0u64; 4], vec![0u64; 16], vec![0u64; 256]];297 let mut zeroed = [vec![0u64; 4], vec![0u64; 16], vec![0u64; 256]];298 for n in 1..=limit {299 let (p, z) = digit_profiles(n);300 for k in 1..=3usize {301 let bits = (1usize << (1 << k)) - 1;302 for target in [303 (p[k - 1] as usize, &mut plain[k - 1]),304 (z[k - 1] as usize, &mut zeroed[k - 1]),305 ] {306 let (seed, store) = target;307 let free = bits & !seed;308 let mut sub = free;309 loop {310 store[seed | sub] += 1;311 if sub == 0 {312 break;313 }314 sub = (sub - 1) & free;315 }316 }317 }318 }319 (plain, zeroed)320}321322fn accepts_agrees(limit: usize) {323 for n in 1..=limit {324 let bits = usize::BITS as usize - n.leading_zeros() as usize;325 let word: Vec<usize> = (0..bits).rev().map(|i| (n >> i) & 1).collect();326 let (p, _) = digit_profiles(n);327 for k in 1..=3usize {328 let codes = 1usize << (1 << k);329 for code in 0..codes {330 let rule = Rule::new(1, k, code as u64).unwrap();331 let want = rule.accepts(&word);332 let have = (p[k - 1] as usize) & !code == 0;333 assert_eq!(334 want, have,335 "profile membership against Rule::accepts at n {n} k {k} code {code}"336 );337 }338 }339 }340}341342// THE REPORT343344fn ratio(value: i64, mass: u64) -> String {345 if mass == 0 {346 return "-".to_string();347 }348 format!("{:.6}", value as f64 / (mass as f64).sqrt())349}350351struct Facts {352 windows: usize,353 rho: f64,354 coupling: f64,355 closed: bool,356 zero_mass: u64,357 mass: u64,358 meter: i64,359 peak: i64,360}361362fn print_rows(tag: &str, width: usize, code: usize, sweep: &Sweep) {363 let slot = width - 1;364 for (index, (label, _)) in sweep.phases.iter().enumerate() {365 let a = sweep.read[slot].mass[index][code];366 let m = sweep.read[slot].meter[index][code];367 let p = sweep.read[slot].peak[index][code];368 println!(369 "{tag} k={width} code={code} x={label} A={a} M={m} Mmax={p} r={} rmax={}",370 ratio(m, a),371 ratio(p, a)372 );373 }374}375376fn main() {377 let started = Instant::now();378 let depth: u32 = env::args()379 .nth(1)380 .and_then(|a| a.parse().ok())381 .unwrap_or(30);382 let cap = 1usize << depth;383 let mu = mobius_sieve(cap);384 let sieved = started.elapsed().as_secs_f64();385 let mut profile = vec![0u8; cap + 1];386 for n in 4..=cap {387 profile[n] = profile[n >> 1] | (1u8 << (n & 7));388 }389 let anchors: Vec<usize> = (1..=8)390 .map(|e| 10usize.pow(e))391 .filter(|&x| x <= cap)392 .collect();393 let phases = phase_grid(depth, cap);394 let walked = Instant::now();395 let sweep = sweep(&mu, &profile, &phases, &anchors);396 let walk = walked.elapsed().as_secs_f64();397 drop(profile);398399 let mertens_want = [-1i64, 1, 2, -23, -48, 212, 1037, 1928];400 assert_eq!(401 sweep.mertens,402 mertens_want[..anchors.len()].to_vec(),403 "the full line against A084237"404 );405 println!(406 "control mertens anchors 10^1..10^{} = {:?} (A084237)",407 anchors.len(),408 sweep.mertens409 );410411 let designs: [(&str, &[u64], &[(usize, i64, i64)]); 3] = [412 (413 "0 1",414 &[0, 1],415 &[(14, 11, 105), (16, 149, 173), (18, -30, 312)],416 ),417 ("1 2", &[1, 2], &[(18, -1461, 1582)]),418 ("0 2", &[0, 2], &[]),419 ];420 for (name, digits, pins) in designs {421 for levels in [14usize, 16, 18, 20] {422 let top = digits.iter().copied().max().unwrap() as usize423 * (3usize.pow(levels as u32) - 1)424 / 2;425 if top > cap {426 println!("control design base3 digits {name} L={levels} needs {top} which is past 2^{depth}, not pinned here");427 continue;428 }429 let (m, p, a) = design_meter(&mu, levels, digits);430 let pinned = pins.iter().find(|(l, _, _)| *l == levels);431 println!(432 "control design base3 digits {name} L={levels} A={a} M={m} Mmax={p} pinned={}",433 if pinned.is_some() { "yes" } else { "no" }434 );435 if let Some((_, want_m, want_p)) = pinned {436 assert_eq!(437 (m, p),438 (*want_m, *want_p),439 "base 3 design against lab/rs/mobius-designs"440 );441 }442 }443 }444445 let limit = (1usize << 20).min(cap);446 let (plain, zeroed) = recount(limit);447 accepts_agrees(1 << 12);448 let cut = phases.iter().position(|(_, x)| *x == limit).unwrap();449 for slot in 0..3 {450 assert_eq!(451 sweep.read[slot].mass[cut], plain[slot],452 "the profile recurrence against a direct digit recount at 2^20"453 );454 }455 println!("control recount 2^20 agrees with the profile recurrence on all 276 rules");456 println!("control Rule::accepts agrees with profile containment on all 276 rules below 2^12");457458 let last = sweep.phases.len() - 1;459 let mut facts: Vec<Vec<Facts>> = Vec::new();460 for width in 1..=3usize {461 let slot = width - 1;462 let codes = 1usize << (1 << width);463 let mut row = Vec::new();464 for code in 0..codes {465 let rule = Rule::new(1, width, code as u64).unwrap();466 row.push(Facts {467 windows: allowed_windows(&rule),468 rho: perron(&rule),469 coupling: kappa(&rule),470 closed: plain[slot][code] == zeroed[slot][code],471 zero_mass: zeroed[slot][code],472 mass: sweep.read[slot].mass[last][code],473 meter: sweep.read[slot].meter[last][code],474 peak: sweep.read[slot].peak[last][code],475 });476 }477 facts.push(row);478 }479480 let reps = representatives(3);481 let class_reps: Vec<usize> = (0..256).filter(|&c| reps[c] == c).collect();482 println!(483 "classes k=3 representatives={} (G_(1,3) orbits)",484 class_reps.len()485 );486487 let full: Vec<f64> = (0..sweep.phases.len())488 .map(|i| sweep.read[0].peak[i][3] as f64 / (sweep.read[0].mass[i][3] as f64).sqrt())489 .collect();490 let band_low = full.iter().cloned().fold(f64::INFINITY, f64::min);491 let band_high = full.iter().cloned().fold(f64::NEG_INFINITY, f64::max);492 let strays = |floor: u64| -> Vec<(usize, usize, u64)> {493 (1..=3usize)494 .flat_map(|w| (0..(1usize << (1 << w))).map(move |c| (w, c)))495 .filter(|&(w, c)| facts[w - 1][c].mass >= floor)496 .filter(|&(width, code)| {497 let slot = width - 1;498 (0..sweep.phases.len()).all(|index| {499 let a = sweep.read[slot].mass[index][code];500 if a < floor {501 return true;502 }503 let r = sweep.read[slot].peak[index][code] as f64 / (a as f64).sqrt();504 r < band_low || r > band_high505 })506 })507 .map(|(w, c)| (w, c, facts[w - 1][c].mass))508 .collect()509 };510511 let last_label = sweep.phases[last].0.clone();512 let score =513 |w: usize, c: usize| facts[w - 1][c].peak as f64 / (facts[w - 1][c].mass as f64).sqrt();514 let sweepwide = |w: usize, c: usize| {515 (0..sweep.phases.len())516 .filter(|&i| sweep.read[w - 1].mass[i][c] > 0)517 .map(|i| {518 (519 sweep.read[w - 1].peak[i][c] as f64520 / (sweep.read[w - 1].mass[i][c] as f64).sqrt(),521 i,522 )523 })524 .fold(525 (0.0f64, 0usize),526 |best, next| {527 if next.0 > best.0 {528 next529 } else {530 best531 }532 },533 )534 };535 let against = |w: usize, c: usize| {536 (0..sweep.phases.len())537 .filter(|&i| sweep.read[w - 1].mass[i][c] > 0)538 .map(|i| {539 let mine = sweep.read[w - 1].peak[i][c] as f64540 / (sweep.read[w - 1].mass[i][c] as f64).sqrt();541 let line =542 sweep.read[0].peak[i][3] as f64 / (sweep.read[0].mass[i][3] as f64).sqrt();543 (mine / line, i)544 })545 .fold(546 (0.0f64, 0usize),547 |best, next| {548 if next.0 > best.0 {549 next550 } else {551 best552 }553 },554 )555 };556 let rank = |floor: u64| {557 let mut out: Vec<(usize, usize)> = (1..=3usize)558 .flat_map(|w| (0..(1usize << (1 << w))).map(move |c| (w, c)))559 .filter(|&(w, c)| facts[w - 1][c].mass >= floor)560 .collect();561 out.sort_by(|&(wa, a), &(wb, b)| {562 score(wb, b)563 .partial_cmp(&score(wa, a))564 .unwrap()565 .then((wa, a).cmp(&(wb, b)))566 });567 out568 };569 let mut named: Vec<usize> = vec![23, 54, 126, 127, 255];570 named.extend(571 strays(10_000)572 .iter()573 .filter(|(w, _, _)| *w == 3)574 .map(|(_, c, _)| *c),575 );576 for floor in [64u64, 10_000] {577 let list = rank(floor);578 let high: Vec<(usize, usize)> = list.iter().take(10).copied().collect();579 let low: Vec<(usize, usize)> = list.iter().rev().take(10).copied().collect();580 named.extend(581 high.iter()582 .chain(low.iter())583 .filter(|(w, _)| *w == 3)584 .map(|(_, c)| *c),585 );586 println!(587 "top all widths rmax at phase {last_label} over {} rules of at least {floor} elements: {}",588 list.len(),589 high.iter()590 .map(|&(w, c)| format!("k{w}code{c}:{:.6}", score(w, c)))591 .collect::<Vec<_>>()592 .join(",")593 );594 println!(595 "bottom all widths rmax at phase {last_label} over {} rules of at least {floor} elements: {}",596 list.len(),597 low.iter()598 .map(|&(w, c)| format!("k{w}code{c}:{:.6}", score(w, c)))599 .collect::<Vec<_>>()600 .join(",")601 );602 let span_low = list603 .iter()604 .map(|&(w, c)| score(w, c))605 .fold(f64::INFINITY, f64::min);606 let span_high = list607 .iter()608 .map(|&(w, c)| score(w, c))609 .fold(f64::NEG_INFINITY, f64::max);610 println!(611 "span all widths floor={floor} rules={} rmax at phase {last_label} runs [{span_low:.6}, {span_high:.6}]",612 list.len()613 );614 let mut wide: Vec<(f64, usize, usize, usize)> = list615 .iter()616 .map(|&(w, c)| {617 let (v, i) = sweepwide(w, c);618 (v, w, c, i)619 })620 .collect();621 wide.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap());622 let wide_low = wide.iter().map(|r| r.0).fold(f64::INFINITY, f64::min);623 println!(624 "span sweepwide floor={floor} rules={} rmax over all {} phases runs [{wide_low:.6}, {:.6}], the top on k{} code {} at phase {}",625 wide.len(),626 sweep.phases.len(),627 wide[0].0,628 wide[0].1,629 wide[0].2,630 sweep.phases[wide[0].3].0631 );632 println!(633 "wide top ten sweepwide rmax floor={floor}: {}",634 wide.iter()635 .take(10)636 .map(|r| format!("k{}code{}:{:.6}@{}", r.1, r.2, r.0, sweep.phases[r.3].0))637 .collect::<Vec<_>>()638 .join(",")639 );640 let tail = sweep.phases.len() - sweep.phases.len() / 4;641 let late = wide.iter().filter(|r| r.3 >= tail).count();642 println!(643 "latepeak floor={floor} rules whose sweepwide rmax falls in the last quarter of the phases, from {}: {late} of {}",644 sweep.phases[tail].0,645 wide.len()646 );647648 let mut factors: Vec<(f64, usize, usize, usize)> = list649 .iter()650 .map(|&(w, c)| {651 let (v, i) = against(w, c);652 (v, w, c, i)653 })654 .collect();655 factors.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap());656 println!(657 "factor same phase against the full line, floor={floor}, largest over all phases: {}",658 factors659 .iter()660 .take(10)661 .map(|r| format!("k{}code{}:{:.6}@{}", r.1, r.2, r.0, sweep.phases[r.3].0))662 .collect::<Vec<_>>()663 .join(",")664 );665 let at_last: Vec<f64> = list666 .iter()667 .map(|&(w, c)| score(w, c) / score(1, 3))668 .collect();669 println!(670 "factor same phase at {last_label} floor={floor} runs [{:.6}, {:.6}] over {} rules",671 at_last.iter().cloned().fold(f64::INFINITY, f64::min),672 at_last.iter().cloned().fold(f64::NEG_INFINITY, f64::max),673 at_last.len()674 );675 }676677 for width in 1..=3usize {678 let slot = width - 1;679 let codes = 1usize << (1 << width);680 for code in 0..codes {681 let f = &facts[slot][code];682 let mut track = Vec::new();683 for (index, (label, _)) in sweep.phases.iter().enumerate() {684 if !label.ends_with(".00") {685 continue;686 }687 let level: u32 = label[..label.len() - 3].parse().unwrap();688 if level % 4 != 0 && level != depth {689 continue;690 }691 let a = sweep.read[slot].mass[index][code];692 let p = sweep.read[slot].peak[index][code];693 track.push(format!("{level}:{}", ratio(p, a)));694 }695 println!(696 "rule k={width} code={code} rep={} W={} rho={:.9} kappa={:.6} zeroclosed={} zeroA={} A={} M={} Mmax={} r={} rmax={} track={}",697 if width == 3 { reps[code] } else { code },698 f.windows,699 f.rho,700 f.coupling,701 if f.closed { "yes" } else { "no" },702 f.zero_mass,703 f.mass,704 f.meter,705 f.peak,706 ratio(f.meter, f.mass),707 ratio(f.peak, f.mass),708 track.join(",")709 );710 }711 }712713 for code in 0..4usize {714 print_rows("row", 1, code, &sweep);715 }716 for code in 0..16usize {717 print_rows("row", 2, code, &sweep);718 }719 named.sort();720 named.dedup();721 for code in &named {722 print_rows("row", 3, *code, &sweep);723 }724725 let fibbinary: Vec<u64> = vec![726 1, 2, 4, 5, 8, 9, 10, 16, 17, 18, 20, 21, 32, 33, 34, 36, 37, 40, 41, 42,727 ];728 let golden: Vec<u64> = (1..=64u64)729 .filter(|&n| Rule::new(1, 2, 7).unwrap().accepts(&binary_word(n)))730 .collect();731 assert_eq!(732 golden[..fibbinary.len()],733 fibbinary[..],734 "code 7 at k=2 against A003714"735 );736 println!(737 "control code 7 at k=2 opens {:?} which is A003714 without its zero",738 &golden[..12]739 );740 let mersenne: Vec<u64> = (1..=1024u64)741 .filter(|&n| Rule::new(1, 2, 11).unwrap().accepts(&binary_word(n)))742 .collect();743 assert_eq!(744 mersenne,745 vec![1, 3, 7, 15, 31, 63, 127, 255, 511, 1023],746 "code 11 at k=2 is A000225"747 );748 assert_eq!(749 facts[1][11].mass, depth as u64,750 "code 11 at k=2 holds one element per level"751 );752 println!(753 "control code 11 at k=2 opens the Mersenne numbers A000225, A=2^m-1 count {}",754 facts[1][11].mass755 );756757 let mut climbing: Vec<(usize, usize)> = Vec::new();758 let mut climbing_meter: Vec<(usize, usize)> = Vec::new();759 for width in 1..=3usize {760 let slot = width - 1;761 let codes = 1usize << (1 << width);762 for code in 0..codes {763 if facts[slot][code].mass < 1000 {764 continue;765 }766 let mut up = true;767 for index in (last - 7)..=last {768 let a = sweep.read[slot].mass[index][code] as f64;769 let p = sweep.read[slot].peak[index][code] as f64;770 let b = sweep.read[slot].mass[index - 1][code] as f64;771 let q = sweep.read[slot].peak[index - 1][code] as f64;772 if p / a.sqrt() <= q / b.sqrt() {773 up = false;774 break;775 }776 }777 if up {778 climbing.push((width, code));779 }780 let mut rises = true;781 for index in (last - 7)..=last {782 let a = sweep.read[slot].mass[index][code] as f64;783 let p = sweep.read[slot].meter[index][code].abs() as f64;784 let b = sweep.read[slot].mass[index - 1][code] as f64;785 let q = sweep.read[slot].meter[index - 1][code].abs() as f64;786 if p / a.sqrt() <= q / b.sqrt() {787 rises = false;788 break;789 }790 }791 if rises {792 climbing_meter.push((width, code));793 }794 }795 }796 println!(797 "climbing rmax rises at every one of the last eight phases on {} rules: {:?}",798 climbing.len(),799 climbing800 );801 println!(802 "climbing abs r rises at every one of the last eight phases on {} rules: {:?}",803 climbing_meter.len(),804 climbing_meter805 );806807 let line = &facts[0][3];808 println!(809 "band fullline k=1 code=3 A={} M={} Mmax={} r={} rmax={}",810 line.mass,811 line.meter,812 line.peak,813 ratio(line.meter, line.mass),814 ratio(line.peak, line.mass)815 );816 for width in 2..=3usize {817 let slot = width - 1;818 let codes = 1usize << (1 << width);819 let edges = [0.0f64, 1e-9, 0.1, 0.2, 0.3, 0.5, 0.8, 2.0];820 for floor in [64u64, 10_000] {821 for pair in edges.windows(2) {822 let (lo, hi) = (pair[0], pair[1]);823 let members: Vec<usize> = (0..codes)824 .filter(|&c| {825 let f = &facts[slot][c];826 f.mass >= floor && f.rho > 0.0 && f.coupling >= lo && f.coupling < hi827 })828 .collect();829 if members.is_empty() {830 continue;831 }832 let scores: Vec<f64> = members833 .iter()834 .map(|&c| facts[slot][c].peak as f64 / (facts[slot][c].mass as f64).sqrt())835 .collect();836 let mean = scores.iter().sum::<f64>() / scores.len() as f64;837 let least = scores.iter().cloned().fold(f64::INFINITY, f64::min);838 let most = scores.iter().cloned().fold(f64::NEG_INFINITY, f64::max);839 let argmax = members[scores.iter().position(|&s| s == most).unwrap()];840 let argmin = members[scores.iter().position(|&s| s == least).unwrap()];841 println!(842 "band k={width} floor={floor} kappa=[{lo},{hi}) rules={} rmax_mean={:.6} rmax_min={:.6} on code {argmin} rmax_max={:.6} on code {argmax}",843 members.len(),844 mean,845 least,846 most847 );848 }849 }850 let mut by_rho: Vec<(f64, usize)> = (0..codes)851 .filter(|&c| facts[slot][c].mass >= 64 && facts[slot][c].rho > 0.0)852 .map(|c| ((facts[slot][c].rho * 1e4).round() / 1e4, c))853 .collect();854 by_rho.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap().then(a.1.cmp(&b.1)));855 let mut seen: Vec<f64> = Vec::new();856 for (r, c) in &by_rho {857 if seen.contains(r) {858 continue;859 }860 seen.push(*r);861 let members: Vec<usize> = by_rho862 .iter()863 .filter(|(q, _)| q == r)864 .map(|(_, c)| *c)865 .collect();866 let scores: Vec<f64> = members867 .iter()868 .map(|&m| facts[slot][m].peak as f64 / (facts[slot][m].mass as f64).sqrt())869 .collect();870 let mean = scores.iter().sum::<f64>() / scores.len() as f64;871 println!(872 "rho k={width} rho={r:.4} rules={} least={c} rmax_mean={:.6} rmax_min={:.6} rmax_max={:.6}",873 members.len(),874 mean,875 scores.iter().cloned().fold(f64::INFINITY, f64::min),876 scores.iter().cloned().fold(f64::NEG_INFINITY, f64::max)877 );878 }879 }880881 let mut running = 0i64;882 let mut ceiling = 0i64;883 let mut early = Vec::new();884 for n in 1..=400usize {885 running += mu[n] as i64;886 if running.abs() > ceiling {887 ceiling = running.abs();888 }889 if [1usize, 5, 13, 31, 200, 256].contains(&n) {890 early.push(format!("{n}:{:.6}", ceiling as f64 / (n as f64).sqrt()));891 }892 }893 println!(894 "gridstart the full line rmax below the grid reads {}",895 early.join(",")896 );897898 let mut only_left = 0u64;899 let mut only_right = 0u64;900 let mut both = 0u64;901 let mut first_gap = 0usize;902 for n in 1..=(1usize << 20) {903 let (p, _) = digit_profiles(n);904 let left = (p[1] as usize) & !14usize == 0;905 let right = (p[2] as usize) & !126usize == 0;906 match (left, right) {907 (true, true) => both += 1,908 (true, false) => only_left += 1,909 (false, true) => {910 only_right += 1;911 if first_gap == 0 {912 first_gap = n;913 }914 }915 _ => {}916 }917 }918 println!(919 "pair k=2 code 14 against k=3 code 126 below 2^20: shared={both} only14={only_left} only126={only_right} symmetric={} least in 126 not 14 is {first_gap}",920 only_left + only_right921 );922923 for width in 1..=3usize {924 let codes = 1usize << (1 << width);925 let mut agree = Vec::new();926 for code in 0..codes {927 let rule = Rule::new(1, width, code as u64).unwrap();928 let mut same = true;929 'words: for length in 1..=14usize {930 for value in 0..(1usize << length) {931 let word: Vec<usize> = (0..length).rev().map(|i| (value >> i) & 1).collect();932 let trimmed: Vec<usize> =933 word.iter().copied().skip_while(|&d| d == 0).collect();934 if rule.accepts(&word) != rule.accepts(&trimmed) {935 same = false;936 break 'words;937 }938 }939 }940 if same {941 agree.push(code);942 }943 }944 println!(945 "reading k={width} the word language agrees with the integer set on {} of {codes} codes over every word to length 14: {:?}",946 agree.len(),947 agree948 );949 }950951 let edges = [0.0f64, 1e-9, 0.1, 0.2, 0.3, 0.5, 0.8, 2.0];952 for floor in [64u64, 10_000] {953 for only in [false, true] {954 for pair in edges.windows(2) {955 let (lo, hi) = (pair[0], pair[1]);956 let members: Vec<(usize, usize)> = (1..=3usize)957 .flat_map(|w| (0..(1usize << (1 << w))).map(move |c| (w, c)))958 .filter(|&(w, c)| {959 let f = &facts[w - 1][c];960 (!only || w == 3)961 && f.mass >= floor962 && f.rho > 0.0963 && f.coupling >= lo964 && f.coupling < hi965 })966 .collect();967 if members.is_empty() {968 continue;969 }970 let scores: Vec<f64> = members.iter().map(|&(w, c)| score(w, c)).collect();971 let mean = scores.iter().sum::<f64>() / scores.len() as f64;972 let least = scores.iter().cloned().fold(f64::INFINITY, f64::min);973 let most = scores.iter().cloned().fold(f64::NEG_INFINITY, f64::max);974 let argmax = members[scores.iter().position(|&s| s == most).unwrap()];975 let argmin = members[scores.iter().position(|&s| s == least).unwrap()];976 println!(977 "kappaband widths={} floor={floor} kappa=[{lo},{hi}) rules={} rmax_mean={mean:.6} at phase {last_label} rmax_min={least:.6} on k{} code {} rmax_max={most:.6} on k{} code {} members={}",978 if only { "3" } else { "1,2,3" },979 members.len(),980 argmin.0,981 argmin.1,982 argmax.0,983 argmax.1,984 members985 .iter()986 .map(|&(w, c)| format!("k{w}c{c}"))987 .collect::<Vec<_>>()988 .join(" ")989 );990 }991 }992 }993994 let top_at = full.iter().position(|&v| v == band_high).unwrap();995 let low_at = full.iter().position(|&v| v == band_low).unwrap();996 println!(997 "band control the full line rmax runs [{band_low:.6}, {band_high:.6}] over every phase, the top at phase {} and the floor at phase {}",998 sweep.phases[top_at].0, sweep.phases[low_at].0999 );1000 for floor in [64u64, 10_000] {1001 let census = (1..=3usize)1002 .flat_map(|w| (0..(1usize << (1 << w))).map(move |c| (w, c)))1003 .filter(|&(w, c)| facts[w - 1][c].mass >= floor)1004 .count();1005 let outside = strays(floor);1006 println!(1007 "band outside floor={floor} census={census} rules never inside the control band at any phase where they hold {floor} elements: {} {:?}",1008 outside.len(),1009 outside1010 );1011 }10121013 let closed_counts: Vec<usize> = (0..3)1014 .map(|slot| {1015 let codes = 1usize << (1 << (slot + 1));1016 (0..codes).filter(|&c| facts[slot][c].closed).count()1017 })1018 .collect();1019 println!("zeroclosed counts k=1,2,3 = {closed_counts:?} of 4, 16, 256 tested below 2^20");1020 for width in 1..=3usize {1021 let slot = width - 1;1022 let codes = 1usize << (1 << width);1023 for code in 0..codes {1024 let want = match width {1025 1 => code == 0 || code & 1 == 1,1026 2 => (code >> 1) & 1 == 1,1027 _ => (code >> 2) & 1 == 1 && (code >> 3) & 1 == 1,1028 };1029 assert_eq!(1030 facts[slot][code].closed, want,1031 "the zero-closed criterion at k {width} code {code}"1032 );1033 }1034 }1035 println!("zeroclosed criterion holds code for code: k=1 the codes allowing the digit 0 with the empty code, k=2 the codes allowing 01, k=3 the codes allowing both 010 and 011");1036 println!(1037 "run depth={depth} N={cap} sieve={sieved:.1}s sweep={walk:.1}s total={:.1}s phases={}",1038 started.elapsed().as_secs_f64(),1039 sweep.phases.len()1040 );1041}10421043fn binary_word(n: u64) -> Vec<usize> {1044 let bits = 64 - n.leading_zeros() as usize;1045 (0..bits).rev().map(|i| ((n >> i) & 1) as usize).collect()1046}