main.rs
18.1 kB · rust · 700 lines
1use std::collections::HashSet;2use std::env;3use std::time::Instant;45// BITSET67struct Bits {8 top: u64,9 w: Vec<u64>,10}1112impl Bits {13 fn new(top: u64) -> Bits {14 let words = (top / 64 + 1) as usize;15 Bits {16 top,17 w: vec![0; words],18 }19 }2021 fn set(&mut self, x: u64) {22 self.w[(x / 64) as usize] |= 1u64 << (x % 64);23 }2425 fn get(&self, x: u64) -> bool {26 (self.w[(x / 64) as usize] >> (x % 64)) & 1 == 127 }2829 fn mask(&mut self) {30 let r = self.top % 64;31 let last = self.w.len() - 1;32 if r != 63 {33 self.w[last] &= (1u64 << (r + 1)) - 1;34 }35 }3637 fn or_shift(&mut self, sh: u64) {38 if sh > self.top {39 return;40 }41 let wsh = (sh / 64) as usize;42 let bsh = (sh % 64) as u32;43 let len = self.w.len();44 for i in (wsh..len).rev() {45 let src = i - wsh;46 let mut v = self.w[src] << bsh;47 if bsh > 0 && src > 0 {48 v |= self.w[src - 1] >> (64 - bsh);49 }50 self.w[i] |= v;51 }52 self.mask();53 }5455 fn count(&self) -> u64 {56 self.w.iter().map(|x| x.count_ones() as u64).sum()57 }5859 fn count_to(&self, d: u64) -> u64 {60 let full = (d / 64) as usize;61 let head: u64 = self.w[..full].iter().map(|x| x.count_ones() as u64).sum();62 let r = d % 64;63 let tail = if r == 63 {64 self.w[full]65 } else {66 self.w[full] & ((1u64 << (r + 1)) - 1)67 };68 head + tail.count_ones() as u6469 }70}7172// DESIGNS7374fn powers(base: u64, top: u64) -> Vec<u64> {75 let mut v = vec![1u64];76 while v77 .last()78 .unwrap()79 .checked_mul(base)80 .map_or(false, |p| p <= top)81 {82 v.push(v.last().unwrap() * base);83 }84 v85}8687fn members(base: u64, top: u64) -> Vec<u64> {88 let mut list = vec![0u64];89 for p in powers(base, top) {90 let n = list.len();91 for i in 0..n {92 let x = list[i] + p;93 if x <= top {94 list.push(x);95 }96 }97 }98 list.sort_unstable();99 list100}101102fn sumset(top: u64, direct: u64, shifted: u64) -> Bits {103 let mut s = Bits::new(top);104 for a in members(direct, top) {105 s.set(a);106 }107 for p in powers(shifted, top) {108 s.or_shift(p);109 }110 s111}112113fn sumset_direct(top: u64) -> Bits {114 let mut s = Bits::new(top);115 let a = members(3, top);116 let b = members(4, top);117 for &x in &a {118 for &y in &b {119 if x + y <= top {120 s.set(x + y);121 }122 }123 }124 s125}126127// SCAN128129#[derive(Clone, Copy)]130struct Ext {131 c: u64,132 x: u64,133}134135impl Ext {136 fn none() -> Ext {137 Ext { c: 0, x: 0 }138 }139140 fn above(&self, c: u64, x: u64) -> bool {141 self.x == 0 || (c as u128) * (self.x as u128) > (self.c as u128) * (x as u128)142 }143144 fn below(&self, c: u64, x: u64) -> bool {145 self.x == 0 || (c as u128) * (self.x as u128) < (self.c as u128) * (x as u128)146 }147}148149#[derive(Clone, Copy)]150struct Window {151 lo: u64,152 hi: u64,153 max: Ext,154 min: Ext,155}156157impl Window {158 fn new(lo: u64, hi: u64) -> Window {159 Window {160 lo,161 hi,162 max: Ext::none(),163 min: Ext::none(),164 }165 }166167 fn feed(&mut self, c: u64, x: u64, up: bool, down: bool) {168 if up && self.max.above(c, x) {169 self.max = Ext { c, x };170 }171 if down && self.min.below(c, x) {172 self.min = Ext { c, x };173 }174 }175}176177struct Scan {178 marks: Vec<(u64, u64)>,179 windows3: Vec<Window>,180 windows2: Vec<Window>,181 global: Window,182}183184fn window_list(base: u64, top: u64) -> Vec<Window> {185 let p = powers(base, top);186 let mut v = Vec::new();187 for i in 0..p.len() {188 let lo = p[i];189 let hi = if i + 1 < p.len() { p[i + 1] - 1 } else { top };190 v.push(Window::new(lo, hi));191 }192 v193}194195fn scan(s: &Bits, marks: &[u64]) -> Scan {196 let top = s.top;197 let windows3 = window_list(3, top);198 let windows2 = window_list(2, top);199 let mut special: HashSet<u64> = HashSet::new();200 for &m in marks {201 special.insert(m / 64);202 }203 for w in windows3.iter().chain(windows2.iter()) {204 special.insert(w.lo / 64);205 special.insert(w.hi / 64);206 if w.lo > 0 {207 special.insert((w.lo - 1) / 64);208 }209 }210 special.insert(0);211 special.insert(top / 64);212 let mut markset: Vec<u64> = marks.to_vec();213 markset.sort_unstable();214 markset.dedup();215 let mut sc = Scan {216 marks: Vec::new(),217 windows3,218 windows2,219 global: Window::new(1, top),220 };221 let mut i3 = 0usize;222 let mut i2 = 0usize;223 let mut mi = 0usize;224 let mut c: u64 = 0;225 for (wi, &word) in s.w.iter().enumerate() {226 let base = (wi as u64) * 64;227 if !special.contains(&(wi as u64)) {228 if word == 0 {229 let x = base + 63;230 sc.windows3[i3].feed(c, x, false, true);231 sc.windows2[i2].feed(c, x, false, true);232 sc.global.feed(c, x, false, true);233 } else if word == u64::MAX {234 c += 64;235 let x = base + 63;236 sc.windows3[i3].feed(c, x, true, false);237 sc.windows2[i2].feed(c, x, true, false);238 sc.global.feed(c, x, true, false);239 } else {240 for j in 0..64u64 {241 let x = base + j;242 let on = (word >> j) & 1 == 1;243 if on {244 c += 1;245 }246 sc.windows3[i3].feed(c, x, on, !on);247 sc.windows2[i2].feed(c, x, on, !on);248 sc.global.feed(c, x, on, !on);249 }250 }251 continue;252 }253 for j in 0..64u64 {254 let x = base + j;255 if x == 0 {256 continue;257 }258 if x > top {259 break;260 }261 if (word >> j) & 1 == 1 {262 c += 1;263 }264 while i3 + 1 < sc.windows3.len() && x >= sc.windows3[i3 + 1].lo {265 i3 += 1;266 }267 while i2 + 1 < sc.windows2.len() && x >= sc.windows2[i2 + 1].lo {268 i2 += 1;269 }270 sc.windows3[i3].feed(c, x, true, true);271 sc.windows2[i2].feed(c, x, true, true);272 sc.global.feed(c, x, true, true);273 while mi < markset.len() && markset[mi] == x {274 sc.marks.push((x, c));275 mi += 1;276 }277 }278 }279 sc280}281282// PRINT283284fn dens(c: u64, x: u64) -> String {285 let q = (c as u128) * 1_000_000u128 / (x as u128);286 format!("{}.{:06}", q / 1_000_000, q % 1_000_000)287}288289fn expo(c: u64, x: u64) -> String {290 if x < 2 || c == 0 {291 return "-".to_string();292 }293 let e = (c as f64).ln() / (x as f64).ln();294 let t = (e * 1_000_000.0).floor() / 1_000_000.0;295 format!("{:.6}", t)296}297298fn print_window(tag: &str, i: usize, w: &Window) {299 println!(300 "{} {:>2} [{}, {}] max {} at x = {} c = {} min {} at x = {} c = {}",301 tag,302 i,303 w.lo,304 w.hi,305 dens(w.max.c, w.max.x),306 w.max.x,307 w.max.c,308 dens(w.min.c, w.min.x),309 w.min.x,310 w.min.c311 );312}313314fn centres(top: u64) -> Vec<(u32, u32, u64, bool)> {315 let p3 = powers(3, top * 3);316 let p4 = powers(4, top * 4);317 let mut v = Vec::new();318 for s in 1..p4.len() {319 for r in 1..p3.len() {320 let d = (p3[r] - 1) / 2 + (p4[s] - 1) / 3;321 if d > top {322 continue;323 }324 let clean = p3[r] > d && p4[s] > d;325 if p3[r] < 3 * p4[s] && p4[s] < 3 * p3[r] {326 v.push((r as u32, s as u32, d, clean));327 }328 }329 }330 v.sort_by_key(|t| t.2);331 v332}333334fn density(k: u32) {335 let t0 = Instant::now();336 let top = 3u64.pow(k);337 let s = sumset(top, 3, 4);338 let built = t0.elapsed().as_secs_f64();339 let p3 = powers(3, top);340 let p4 = powers(4, top);341 let mut marks: Vec<u64> = Vec::new();342 marks.extend(p3.iter().copied());343 marks.extend(p4.iter().copied());344 marks.extend(p3.iter().map(|p| p / 2));345 marks.extend(p4.iter().map(|p| p / 3));346 let cs = centres(top);347 marks.extend(cs.iter().map(|t| t.2));348 marks.retain(|&x| x >= 1);349 let sc = scan(&s, &marks);350 let at = |x: u64| -> u64 { sc.marks.iter().find(|m| m.0 == x).map(|m| m.1).unwrap() };351 println!("top 3^{} = {}", k, top);352 println!("members of A up to top {}", members(3, top).len());353 println!("members of B up to top {}", members(4, top).len());354 println!("card(S meet [1, top]) {}", s.count() - 1);355 for (i, &x) in p3.iter().enumerate() {356 let c = at(x);357 println!(358 "3^{:<2} x = {} c = {} D = {} exponent {}",359 i,360 x,361 c,362 dens(c, x),363 expo(c, x)364 );365 }366 for (i, &x) in p4.iter().enumerate() {367 let c = at(x);368 println!(369 "4^{:<2} x = {} c = {} D = {} exponent {}",370 i,371 x,372 c,373 dens(c, x),374 expo(c, x)375 );376 }377 for (i, &p) in p3.iter().enumerate().skip(1) {378 let x = p / 2;379 let c = at(x);380 println!("3^{:<2}/2 x = {} c = {} D = {}", i, x, c, dens(c, x));381 }382 for (i, &p) in p4.iter().enumerate().skip(1) {383 let x = p / 3;384 let c = at(x);385 println!("4^{:<2}/3 x = {} c = {} D = {}", i, x, c, dens(c, x));386 }387 for &(r, sx, d, clean) in &cs {388 let c = at(d);389 println!(390 "centre r = {:<2} s = {:<2} 4^s/3^r = {} d = {} c = {} D = {} {}",391 r,392 sx,393 ratio(r, sx),394 d,395 c,396 dens(c, d),397 if clean { "clean" } else { "mixed" }398 );399 }400 for (i, w) in sc.windows3.iter().enumerate() {401 print_window("window3", i, w);402 }403 for (i, w) in sc.windows2.iter().enumerate() {404 print_window("window2", i, w);405 }406 print_window("global", 0, &sc.global);407 println!(408 "built in {:.2} s, scanned in {:.2} s",409 built,410 t0.elapsed().as_secs_f64() - built411 );412}413414// ENERGY415416fn zeros3(mut t: i64, k: u32) -> Option<u32> {417 let mut z = 0;418 for _ in 0..k {419 let r = t.rem_euclid(3);420 t = (t - r) / 3;421 if r == 2 {422 t += 1;423 } else if r == 0 {424 z += 1;425 }426 }427 if t == 0 {428 Some(z)429 } else {430 None431 }432}433434fn energy(k: u32, m: u32) -> u128 {435 let p4: Vec<i64> = (0..m).map(|i| 4i64.pow(i)).collect();436 let mut e: u128 = 0;437 let mut digits = vec![-1i64; m as usize];438 loop {439 let mut t = 0i64;440 let mut z4 = 0u32;441 for i in 0..m as usize {442 t += digits[i] * p4[i];443 if digits[i] == 0 {444 z4 += 1;445 }446 }447 if let Some(z3) = zeros3(t, k) {448 e += 1u128 << (z3 + z4);449 }450 let mut i = 0usize;451 loop {452 if i == m as usize {453 return e;454 }455 if digits[i] < 1 {456 digits[i] += 1;457 break;458 }459 digits[i] = -1;460 i += 1;461 }462 }463}464465#[cfg(test)]466fn energy_direct(k: u32, m: u32) -> u128 {467 let a = members(3, 3u64.pow(k) - 1);468 let b = members(4, 4u64.pow(m) - 1);469 let d = (3u64.pow(k) - 1) / 2 + (4u64.pow(m) - 1) / 3;470 let mut r = vec![0u64; d as usize + 1];471 for &x in &a {472 for &y in &b {473 r[(x + y) as usize] += 1;474 }475 }476 r.iter().map(|&v| (v as u128) * (v as u128)).sum()477}478479fn six(num: u128, den: u128) -> String {480 let q = num * 1_000_000 / den;481 format!("{}.{:06}", q / 1_000_000, q % 1_000_000)482}483484fn six_up(num: u128, den: u128) -> String {485 let q = (num * 1_000_000 + den - 1) / den;486 format!("{}.{:06}", q / 1_000_000, q % 1_000_000)487}488489fn ratio(r: u32, m: u32) -> String {490 six_up(4u128.pow(m) * 1_000_000, 3u128.pow(r) * 1_000_000)491}492493fn energies(kmax: u32) {494 let t0 = Instant::now();495 let top = 3u64.pow(kmax);496 let s = sumset(top, 3, 4);497 let mut first: Vec<Option<(u32, f64)>> = vec![None, None];498 let mut last: Option<(u32, f64)> = None;499 for (r, m, d, clean) in centres(top) {500 if r < 4 {501 continue;502 }503 let e = energy(r, m);504 let pairs = 1u128 << (r + m);505 let card = s.count_to(d) as u128;506 let range = d as u128 + 1;507 println!(508 "k = {:<2} m = {:<2} {} 4^m/3^k = {} d = {} card = {} fill {} energy {} random {} ratio {} bound {}",509 r,510 m,511 if clean { "clean" } else { "mixed" },512 ratio(r, m),513 d,514 card,515 six(card, range),516 e,517 six(pairs * pairs, range),518 six_up(e * range, pairs * pairs),519 six(pairs * pairs, e * range)520 );521 let q = (e as f64) * (range as f64) / (pairs as f64) / (pairs as f64);522 for (i, start) in [6u32, 11].iter().enumerate() {523 if r >= *start && first[i].is_none() {524 first[i] = Some((r, q));525 }526 }527 if r >= 6 {528 last = Some((r, q));529 }530 }531 for f in first {532 if let (Some((k0, q0)), Some((k1, q1))) = (f, last) {533 let eta = (q1 / q0).ln() / 3f64.ln() / ((k1 - k0) as f64);534 println!(535 "growth of Q from k = {} to k = {} as 3^(eta k): eta = {:.6}",536 k0,537 k1,538 (eta * 1_000_000.0).ceil() / 1_000_000.0539 );540 }541 }542 println!("{:.2} s", t0.elapsed().as_secs_f64());543}544545// CONTROL546547const OEIS_A367090: [u64; 58] = [548 62, 63, 143, 144, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221,549 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240,550 241, 242, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479,551 480,552];553554fn same(a: &Bits, b: &Bits) -> bool {555 a.top == b.top && a.w == b.w556}557558fn symmetric(s: &Bits, d: u64) -> bool {559 (0..=d).all(|x| s.get(x) == s.get(d - x))560}561562fn complement(s: &Bits, n: usize) -> Vec<u64> {563 (1..=s.top).filter(|&x| !s.get(x)).take(n).collect()564}565566fn control() {567 let t0 = Instant::now();568 let top = 3u64.pow(13);569 let direct = sumset_direct(top);570 let shifted = sumset(top, 3, 4);571 println!(572 "double loop against shift-or at 3^13: {}",573 if same(&direct, &shifted) {574 "agree"575 } else {576 "DIFFER"577 }578 );579 let top = 3u64.pow(17);580 let ab = sumset(top, 3, 4);581 let ba = sumset(top, 4, 3);582 println!(583 "A direct with B shifted against B direct with A shifted at 3^17: {}",584 if same(&ab, &ba) { "agree" } else { "DIFFER" }585 );586 let gaps = complement(&ab, OEIS_A367090.len());587 println!(588 "first {} non-members against A367090: {}",589 OEIS_A367090.len(),590 if gaps == OEIS_A367090 {591 "agree"592 } else {593 "DIFFER"594 }595 );596 for (r, s, d, clean) in centres(top) {597 println!(598 "centre r = {:<2} s = {:<2} d = {} {} symmetric {}",599 r,600 s,601 d,602 if clean { "clean" } else { "mixed" },603 symmetric(&ab, d)604 );605 }606 println!("{:.2} s", t0.elapsed().as_secs_f64());607}608609fn main() {610 let args: Vec<String> = env::args().collect();611 match args.get(1).map(|s| s.as_str()) {612 Some("density") => density(args.get(2).and_then(|s| s.parse().ok()).unwrap_or(17)),613 Some("control") => control(),614 Some("energy") => energies(args.get(2).and_then(|s| s.parse().ok()).unwrap_or(17)),615 _ => println!("verbs: density K, energy K, control"),616 }617}618619#[cfg(test)]620mod tests {621 use super::*;622623 #[test]624 fn shift_or_matches_double_loop() {625 let top = 3u64.pow(9);626 assert!(same(&sumset_direct(top), &sumset(top, 3, 4)));627 assert!(same(&sumset_direct(top), &sumset(top, 4, 3)));628 }629630 #[test]631 fn first_gaps_are_a367090() {632 let s = sumset(3u64.pow(7), 3, 4);633 assert_eq!(complement(&s, OEIS_A367090.len()), OEIS_A367090.to_vec());634 }635636 #[test]637 fn member_counts_are_powers_of_two() {638 assert_eq!(members(3, 3u64.pow(10) - 1).len(), 1 << 10);639 assert_eq!(members(4, 4u64.pow(6) - 1).len(), 1 << 6);640 }641642 #[test]643 fn scan_counts_match_popcount() {644 let top = 3u64.pow(9);645 let s = sumset(top, 3, 4);646 let sc = scan(&s, &[top]);647 assert_eq!(sc.marks, vec![(top, s.count() - 1)]);648 }649650 #[test]651 fn scan_extremes_match_direct_search() {652 let top = 3u64.pow(9);653 let s = sumset(top, 3, 4);654 let sc = scan(&s, &[]);655 for w in sc.windows3.iter().chain(sc.windows2.iter()) {656 let mut c = (1..w.lo).filter(|&x| s.get(x)).count() as u64;657 let mut best_max = Ext::none();658 let mut best_min = Ext::none();659 for x in w.lo..=w.hi {660 if s.get(x) {661 c += 1;662 }663 if best_max.above(c, x) {664 best_max = Ext { c, x };665 }666 if best_min.below(c, x) {667 best_min = Ext { c, x };668 }669 }670 assert_eq!((w.max.c, w.max.x), (best_max.c, best_max.x));671 assert_eq!((w.min.c, w.min.x), (best_min.c, best_min.x));672 }673 }674675 #[test]676 fn energy_matches_representation_histogram() {677 for (k, m) in [(4, 3), (6, 5), (9, 7)] {678 assert_eq!(energy(k, m), energy_direct(k, m));679 }680 }681682 #[test]683 fn prefix_count_matches_bit_scan() {684 let s = sumset(3u64.pow(9), 3, 4);685 for d in [0u64, 63, 64, 705, 15302, 19683] {686 assert_eq!(s.count_to(d), (0..=d).filter(|&x| s.get(x)).count() as u64);687 }688 }689690 #[test]691 fn clean_centres_are_symmetric() {692 let top = 3u64.pow(11);693 let s = sumset(top, 3, 4);694 for (_, _, d, clean) in centres(top) {695 if clean {696 assert!(symmetric(&s, d));697 }698 }699 }700}