main.rs
17.5 kB · rust · 600 lines
1use mrlynum::memory::{allowed_windows, kappa, perron, transfer, Rule};2use std::collections::BTreeSet;3use std::time::Instant;45// THE STEP67fn free(n: u64) -> u64 {8 n ^ (n << 1) ^ 19}1011fn defect(n: u64) -> u64 {12 (3 * n + 1) ^ free(n)13}1415fn d_loc(n: u64) -> u32 {16 defect(n).count_ones()17}1819fn carries(n: u64, bits: usize) -> u64 {20 let mut q = 1u64;21 let mut prev = 0u64;22 let mut out = 0u64;23 for i in 0..bits {24 let b = (n >> i) & 1;25 let next = u64::from(b + prev + q >= 2);26 if next == 1 {27 out |= 1 << (i + 1);28 }29 prev = b;30 q = next;31 }32 out33}3435fn t_free(n: u64) -> u64 {36 if n % 2 == 0 {37 n / 238 } else {39 free(n) / 240 }41}4243// THE GENERAL MULTIPLIER4445fn support(q: u64) -> Vec<usize> {46 (0..64).filter(|j| (q >> j) & 1 == 1).collect()47}4849fn differences(supp: &[usize]) -> BTreeSet<usize> {50 let mut out = BTreeSet::new();51 for a in 0..supp.len() {52 for b in (a + 1)..supp.len() {53 out.insert(supp[b] - supp[a]);54 }55 }56 out57}5859fn skeleton(n: u64, supp: &[usize]) -> u64 {60 let mut x = 1u64;61 for &j in supp {62 x ^= n << j;63 }64 x65}6667fn zero_code(supp: &[usize]) -> (usize, u64) {68 let width = supp[supp.len() - 1] + 1;69 let gaps = differences(supp);70 let mut code = 0u64;71 for w in 0..(1usize << width) {72 let mut ok = true;73 for a in 0..width {74 for b in (a + 1)..width {75 let hit = (w >> (width - 1 - a)) & 1 == 1 && (w >> (width - 1 - b)) & 1 == 1;76 if hit && gaps.contains(&(b - a)) {77 ok = false;78 }79 }80 }81 if ok {82 code |= 1 << w;83 }84 }85 (width, code)86}8788fn word(n: u64, length: usize) -> Vec<usize> {89 (0..length).rev().map(|i| ((n >> i) & 1) as usize).collect()90}9192// THE STATISTICS9394fn longest_run(n: u64) -> u32 {95 let mut best = 0;96 let mut run = 0;97 for i in 0..64 {98 if (n >> i) & 1 == 1 {99 run += 1;100 if run > best {101 best = run;102 }103 } else {104 run = 0;105 }106 }107 best108}109110fn v2(n: u64) -> u32 {111 if n == 0 {112 64113 } else {114 n.trailing_zeros()115 }116}117118fn digits(n: u64) -> u32 {119 64 - n.leading_zeros()120}121122// THE DEPTH123124fn depth_of(bits: &[u8]) -> usize {125 let mut q = 1u8;126 let mut prev = 0u8;127 let mut run = 0usize;128 let mut best = 0usize;129 let step = |b: u8, q: &mut u8, prev: &mut u8, run: &mut usize, best: &mut usize| {130 let next = u8::from(b + *prev + *q >= 2);131 if next == 1 {132 *run += 1;133 if *run > *best {134 *best = *run;135 }136 } else {137 *run = 0;138 }139 *prev = b;140 *q = next;141 };142 for &b in bits {143 step(b, &mut q, &mut prev, &mut run, &mut best);144 }145 for _ in 0..2 {146 step(0, &mut q, &mut prev, &mut run, &mut best);147 }148 best149}150151fn splitmix(state: &mut u64) -> u64 {152 *state = state.wrapping_add(0x9e37_79b9_7f4a_7c15);153 let mut z = *state;154 z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);155 z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);156 z ^ (z >> 31)157}158159// THE REPORTS160161fn report_identity() {162 let start = Instant::now();163 let limit = 1u64 << 18;164 let mut bad_carry = 0u64;165 let mut bad_conj = 0u64;166 let mut bad_rule = 0u64;167 for n in 0..limit {168 if carries(n, 22) != defect(n) {169 bad_carry += 1;170 }171 let m = 2 * n + 1;172 if m ^ (2 * m) != 2 * free(n) + 1 {173 bad_conj += 1;174 }175 let mut sixty = 0u64;176 for i in 0..24 {177 let hi = (m >> i) & 1;178 let lo = if i == 0 { 0 } else { (m >> (i - 1)) & 1 };179 sixty |= (hi ^ lo) << i;180 }181 if sixty != (m ^ (2 * m)) {182 bad_rule += 1;183 }184 }185 println!("identity limit 2^18");186 println!("identity carry_word_mismatches {bad_carry}");187 println!("identity conjugacy_mismatches {bad_conj}");188 println!("identity rule60_mismatches {bad_rule}");189 let mut worst = 0usize;190 for length in 2..=40usize {191 let mut u = 1u64;192 for j in 1..length {193 if j % 2 == 0 {194 u |= 1 << j;195 }196 }197 let v = u ^ 1;198 let diff = (3 * u + 1) ^ (3 * v + 1);199 let all = (2..=length).all(|i| (diff >> i) & 1 == 1);200 assert!(all, "radius witness fails at length {length}");201 worst = length;202 }203 println!("identity radius_witness_max_digit {worst}");204 println!("identity seconds {:.2}", start.elapsed().as_secs_f64());205}206207fn report_rules() {208 let start = Instant::now();209 println!("rule q supp gaps width code windows rho kappa");210 for q in [3u64, 5, 7, 9, 11, 15] {211 let supp = support(q);212 let gaps = differences(&supp);213 let (width, code) = zero_code(&supp);214 let rule = Rule::new(1, width, code).expect("rule in range");215 let supp_text: Vec<String> = supp.iter().map(|j| j.to_string()).collect();216 let gap_text: Vec<String> = gaps.iter().map(|j| j.to_string()).collect();217 println!(218 "rule {q} {} {} {width} {code} {} {:.6} {:.6}",219 supp_text.join(","),220 gap_text.join(","),221 allowed_windows(&rule),222 perron(&rule),223 kappa(&rule)224 );225 let limit = 1u64 << 16;226 let mut mismatches = 0u64;227 for n in 0..limit {228 let free_carry = (q * n + 1) ^ skeleton(n, &supp) == 0;229 let padded = word(n, 16 + width);230 let accepted = n % 2 == 0 && rule.accepts(&padded);231 if free_carry != accepted {232 mismatches += 1;233 }234 }235 println!("rule {q} set_equality_mismatches_below_2^16 {mismatches}");236 }237 let golden = Rule::new(1, 2, 7).expect("golden in range");238 let supergolden = Rule::new(1, 3, 23).expect("supergolden in range");239 println!("kappa golden {:.6}", kappa(&golden));240 println!("kappa supergolden {:.6}", kappa(&supergolden));241 println!("kappa golden_transfer {:?}", transfer(&golden));242 println!("rule seconds {:.2}", start.elapsed().as_secs_f64());243}244245fn report_density() {246 let start = Instant::now();247 let pi = [2i64, 1, 1, 2];248 let balance = [249 2 * pi[0] - (pi[0] + pi[1] + pi[2]),250 2 * pi[1] - pi[3],251 2 * pi[2] - pi[0],252 2 * pi[3] - (pi[1] + pi[2] + pi[3]),253 ];254 println!("density stationary_sixths {pi:?}");255 println!("density balance_residuals {balance:?}");256 println!(257 "density carry_on_mass {}/{}",258 pi[1] + pi[3],259 pi.iter().sum::<i64>()260 );261 println!("density digits mean mean_minus_half_L mean_over_L exact_residual");262 for l in 8..=22usize {263 let limit = 1u64 << l;264 let mut total = 0u128;265 for n in 0..limit {266 total += u128::from(d_loc(n));267 }268 let mean = total as f64 / limit as f64;269 let sign = if l % 2 == 0 { 1i128 } else { -1 };270 let closed = 3 * l as i128 * (limit / 2) as i128 + limit as i128 - sign;271 println!(272 "density {l} {mean:.6} {:.6} {:.6} {}",273 mean - l as f64 / 2.0,274 mean / l as f64,275 3 * total as i128 - closed276 );277 }278 println!("density seconds {:.2}", start.elapsed().as_secs_f64());279}280281fn report_free() {282 let start = Instant::now();283 let limit = 1u64 << 20;284 let mut grew = 0u64;285 let mut unreached = 0u64;286 for n in 1..limit {287 if digits(t_free(n)) > digits(n) {288 grew += 1;289 }290 let mut x = n;291 let mut steps = 0;292 while x != 1 && steps < 4096 {293 x = t_free(x);294 steps += 1;295 }296 if x != 1 {297 unreached += 1;298 }299 }300 println!("free limit 2^20");301 println!("free digit_count_increases {grew}");302 println!("free values_not_reaching_one {unreached}");303 println!("free fixed_point_check {}", t_free(1));304 println!("free seconds {:.2}", start.elapsed().as_secs_f64());305}306307fn report_refutations() {308 let start = Instant::now();309 let limit = 1u64 << 16;310 let stats: Vec<(&str, Box<dyn Fn(u64) -> u64>)> = vec![311 ("popcount", Box::new(|n: u64| u64::from(n.count_ones()))),312 ("longest_run", Box::new(|n: u64| u64::from(longest_run(n)))),313 ("v2", Box::new(|n: u64| u64::from(v2(n)))),314 ("digit_count", Box::new(|n: u64| u64::from(digits(n)))),315 (316 "popcount_and_v2",317 Box::new(|n: u64| u64::from(n.count_ones()) * 128 + u64::from(v2(n))),318 ),319 (320 "all_four",321 Box::new(|n: u64| {322 ((u64::from(n.count_ones()) * 128 + u64::from(longest_run(n))) * 128323 + u64::from(v2(n)))324 * 128325 + u64::from(digits(n))326 }),327 ),328 ];329 println!("refute statistic least_pair d_loc_pair disagreeing_pairs_below_2^16");330 for (name, key) in &stats {331 let mut first: std::collections::HashMap<u64, (u64, u32)> =332 std::collections::HashMap::new();333 let mut groups: std::collections::HashMap<u64, std::collections::HashMap<u32, u64>> =334 std::collections::HashMap::new();335 let mut witness: Option<(u64, u64, u32, u32)> = None;336 for n in 1..limit {337 let k = key(n);338 let d = d_loc(n);339 *groups.entry(k).or_default().entry(d).or_insert(0) += 1;340 match first.get(&k) {341 None => {342 first.insert(k, (n, d));343 }344 Some(&(m, e)) => {345 if e != d && witness.is_none() {346 witness = Some((m, n, e, d));347 }348 }349 }350 }351 let mut pairs = 0u64;352 for counts in groups.values() {353 let held: u64 = counts.values().sum();354 let agreeing: u64 = counts.values().map(|c| c * (c - 1) / 2).sum();355 pairs += held * (held - 1) / 2 - agreeing;356 }357 match witness {358 Some((m, n, e, d)) => println!("refute {name} {m},{n} {e},{d} {pairs}"),359 None => println!("refute {name} none none {pairs}"),360 }361 }362 println!("refute seconds {:.2}", start.elapsed().as_secs_f64());363}364365fn report_depth() {366 let start = Instant::now();367 for l in [4usize, 8, 16, 32, 40] {368 let bits = vec![1u8; l];369 assert_eq!(depth_of(&bits), l + 1, "worst case fails at {l}");370 }371 println!("depth worst_case_is_L_plus_1_at_L 4,8,16,32,40");372 let golden = Rule::new(1, 2, 7).expect("golden in range");373 let phi = perron(&golden);374 let base = 2.0 / phi;375 println!("depth phi {phi:.9} two_over_phi {base:.9}");376 println!(377 "depth L samples mean_depth standard_error log_base_L offset increment_per_quadrupling"378 );379 let mut state = 0x5eed_1234_9abc_def0u64;380 let mut previous: Option<(usize, f64)> = None;381 for l in [16usize, 64, 256, 1024, 4096, 16384] {382 let samples = 200_000usize;383 let mut total = 0u64;384 let mut squares = 0u128;385 let mut bits = vec![0u8; l];386 for _ in 0..samples {387 let mut i = 0;388 while i < l {389 let chunk = splitmix(&mut state);390 let take = std::cmp::min(64, l - i);391 for j in 0..take {392 bits[i + j] = ((chunk >> j) & 1) as u8;393 }394 i += take;395 }396 let depth = depth_of(&bits) as u64;397 total += depth;398 squares += u128::from(depth * depth);399 }400 let mean = total as f64 / samples as f64;401 let second = squares as f64 / samples as f64;402 let error = ((second - mean * mean) / (samples as f64 - 1.0)).sqrt();403 let predicted = (l as f64).ln() / base.ln();404 let increment = match previous {405 Some((pl, pm)) if l == 4 * pl => format!("{:.3}", mean - pm),406 _ => "-".to_string(),407 };408 println!(409 "depth {l} {samples} {mean:.3} {error:.3} {predicted:.3} {:.3} {increment}",410 mean - predicted411 );412 previous = Some((l, mean));413 }414 println!(415 "depth predicted_increment_per_quadrupling {:.6}",416 4.0f64.ln() / base.ln()417 );418 println!("depth seconds {:.2}", start.elapsed().as_secs_f64());419}420421fn main() {422 let whole = Instant::now();423 report_identity();424 report_rules();425 report_density();426 report_free();427 report_refutations();428 report_depth();429 println!("run seconds {:.2}", whole.elapsed().as_secs_f64());430}431432// THE TESTS433434#[cfg(test)]435mod tests {436 use super::*;437438 #[test]439 fn the_carry_word_is_the_defect() {440 for n in 0..(1u64 << 16) {441 assert_eq!(carries(n, 20), defect(n));442 }443 }444445 #[test]446 fn the_substitution_clears_the_constant() {447 for n in 0..(1u64 << 16) {448 let m = 2 * n + 1;449 assert_eq!(m ^ (2 * m), 2 * free(n) + 1);450 }451 }452453 #[test]454 fn the_skeleton_is_rule_sixty() {455 for n in 0..(1u64 << 16) {456 let m = 2 * n + 1;457 let mut sixty = 0u64;458 for i in 0..20 {459 let hi = (m >> i) & 1;460 let lo = if i == 0 { 0 } else { (m >> (i - 1)) & 1 };461 sixty |= (hi ^ lo) << i;462 }463 assert_eq!(sixty, m ^ (2 * m));464 }465 }466467 #[test]468 fn one_digit_reaches_every_digit() {469 for length in 2..=40usize {470 let mut u = 1u64;471 for j in 1..length {472 if j % 2 == 0 {473 u |= 1 << j;474 }475 }476 let diff = (3 * u + 1) ^ (3 * (u ^ 1) + 1);477 assert!((2..=length).all(|i| (diff >> i) & 1 == 1));478 }479 }480481 #[test]482 fn the_zero_carry_codes_are_pinned() {483 let pinned = [484 (3u64, 2usize, 7u64),485 (5, 3, 95),486 (7, 3, 23),487 (9, 4, 22015),488 (11, 4, 279),489 (15, 4, 279),490 ];491 for (q, width, code) in pinned {492 assert_eq!(zero_code(&support(q)), (width, code));493 }494 }495496 #[test]497 fn the_zero_carry_set_is_the_rule() {498 for q in [3u64, 5, 7, 9, 11, 15] {499 let supp = support(q);500 let (width, code) = zero_code(&supp);501 let rule = Rule::new(1, width, code).unwrap();502 for n in 0..(1u64 << 14) {503 let carry_free = (q * n + 1) ^ skeleton(n, &supp) == 0;504 assert_eq!(carry_free, n % 2 == 0 && rule.accepts(&word(n, 14 + width)));505 }506 }507 }508509 #[test]510 fn the_two_named_couplings_come_back() {511 let golden = Rule::new(1, 2, 7).unwrap();512 let supergolden = Rule::new(1, 3, 23).unwrap();513 assert!((kappa(&golden) - 0.098_239).abs() < 5e-7);514 assert!((kappa(&supergolden) - 0.115_204).abs() < 5e-7);515 assert!((perron(&golden) - 1.618_033_988_749_895).abs() < 1e-9);516 assert!((perron(&supergolden) - 1.465_571_231_876_768).abs() < 1e-9);517 }518519 #[test]520 fn the_stationary_vector_balances() {521 let pi = [2i64, 1, 1, 2];522 assert_eq!(2 * pi[0], pi[0] + pi[1] + pi[2]);523 assert_eq!(2 * pi[1], pi[3]);524 assert_eq!(2 * pi[2], pi[0]);525 assert_eq!(2 * pi[3], pi[1] + pi[2] + pi[3]);526 assert_eq!(2 * (pi[1] + pi[3]), pi.iter().sum::<i64>());527 }528529 #[test]530 fn the_carry_density_is_one_half() {531 for l in [16usize, 18, 20] {532 let limit = 1u64 << l;533 let mut total = 0u128;534 for n in 0..limit {535 total += u128::from(d_loc(n));536 }537 let sign = if l % 2 == 0 { 1i128 } else { -1 };538 let closed = 3 * l as i128 * (limit / 2) as i128 + limit as i128 - sign;539 assert_eq!(3 * total as i128, closed);540 }541 }542543 #[test]544 fn the_carry_free_map_has_one_cycle() {545 for n in 1..(1u64 << 18) {546 assert!(digits(t_free(n)) <= digits(n));547 let mut x = n;548 let mut steps = 0;549 while x != 1 && steps < 4096 {550 x = t_free(x);551 steps += 1;552 }553 assert_eq!(x, 1);554 }555 }556557 #[test]558 fn the_defect_is_no_function_of_its_statistics() {559 let pinned: [(&str, fn(u64) -> u64, u64, u64, u32, u32); 6] = [560 ("popcount", |n| u64::from(n.count_ones()), 1, 2, 2, 0),561 ("longest_run", |n| u64::from(longest_run(n)), 1, 2, 2, 0),562 ("v2", |n| u64::from(v2(n)), 1, 3, 2, 3),563 ("digit_count", |n| u64::from(digits(n)), 2, 3, 0, 3),564 (565 "popcount_and_v2",566 |n| u64::from(n.count_ones()) * 128 + u64::from(v2(n)),567 3,568 5,569 3,570 4,571 ),572 (573 "all_four",574 |n| {575 ((u64::from(n.count_ones()) * 128 + u64::from(longest_run(n))) * 128576 + u64::from(v2(n)))577 * 128578 + u64::from(digits(n))579 },580 19,581 25,582 3,583 4,584 ),585 ];586 for (_, key, m, n, dm, dn) in pinned {587 assert_eq!(key(m), key(n));588 assert_eq!(d_loc(m), dm);589 assert_eq!(d_loc(n), dn);590 assert_ne!(dm, dn);591 }592 }593594 #[test]595 fn the_worst_depth_is_the_digit_count_plus_one() {596 for l in [1usize, 2, 4, 8, 16, 32, 40] {597 assert_eq!(depth_of(&vec![1u8; l]), l + 1);598 }599 }600}