main.rs
20.9 kB · rust · 693 lines
1use std::cmp::Ordering;2use std::env;3use std::time::Instant;45// PACKING67fn binom(n: usize, k: usize) -> u64 {8 if k > n {9 return 0;10 }11 let k = k.min(n - k);12 let mut r: u128 = 1;13 for i in 0..k {14 r = r * (n - i) as u128 / (i + 1) as u128;15 }16 r as u6417}1819struct Pack {20 off: Vec<u32>,21 wid: Vec<u32>,22 w: usize,23 bits: u32,24}2526fn pack_for(d: usize) -> Pack {27 let mut off = Vec::with_capacity(d + 1);28 let mut wid = Vec::with_capacity(d + 1);29 let mut o = 0u32;30 for i in 0..=d {31 let b = (64 - binom(d, i).leading_zeros()).max(1);32 off.push(o);33 wid.push(b);34 o += b;35 }36 Pack {37 off,38 wid,39 w: (o as usize + 63) / 64,40 bits: o,41 }42}4344fn put(key: &mut [u64], i: usize, pk: &Pack, v: u32) {45 let bit = pk.off[i] as usize;46 let b = pk.wid[i] as usize;47 let w = bit >> 6;48 let off = bit & 63;49 key[w] |= (v as u64) << off;50 if off + b > 64 {51 key[w + 1] |= (v as u64) >> (64 - off);52 }53}5455fn get(key: &[u64], i: usize, pk: &Pack) -> u32 {56 let bit = pk.off[i] as usize;57 let b = pk.wid[i] as usize;58 let w = bit >> 6;59 let off = bit & 63;60 let mut v = key[w] >> off;61 if off + b > 64 {62 v |= key[w + 1] << (64 - off);63 }64 (v & ((1u64 << b) - 1)) as u3265}6667fn addk(dst: &mut [u64], src: &[u64]) {68 let mut c = 0u64;69 for i in 0..dst.len() {70 let (a, o1) = dst[i].overflowing_add(src[i]);71 let (a, o2) = a.overflowing_add(c);72 dst[i] = a;73 c = (o1 as u64) | (o2 as u64);74 }75}7677fn subk(dst: &mut [u64], src: &[u64]) {78 let mut c = 0u64;79 for i in 0..dst.len() {80 let (a, o1) = dst[i].overflowing_sub(src[i]);81 let (a, o2) = a.overflowing_sub(c);82 dst[i] = a;83 c = (o1 as u64) | (o2 as u64);84 }85}8687fn cmpk(a: &[u64], b: &[u64]) -> Ordering {88 for i in (0..a.len()).rev() {89 if a[i] != b[i] {90 return a[i].cmp(&b[i]);91 }92 }93 Ordering::Equal94}9596fn hashk(k: &[u64]) -> u64 {97 let mut h: u64 = 0xcbf2_9ce4_8422_2325;98 for &w in k {99 h = (h ^ w).wrapping_mul(0x100_0000_01b3);100 }101 h ^ (h >> 31)102}103104// TABLE105106struct Table {107 w: usize,108 cap: usize,109 keys: Vec<u64>,110 vals: Vec<i64>,111 len: usize,112}113114impl Table {115 fn new(w: usize, cap: usize) -> Table {116 Table {117 w,118 cap,119 keys: vec![0u64; w * cap],120 vals: vec![0i64; cap],121 len: 0,122 }123 }124125 fn bytes(&self) -> usize {126 self.keys.len() * 8 + self.vals.len() * 8127 }128129 fn empty(&self, s: usize) -> bool {130 self.keys[s * self.w..(s + 1) * self.w]131 .iter()132 .all(|&x| x == 0)133 }134135 fn key(&self, s: usize) -> &[u64] {136 &self.keys[s * self.w..(s + 1) * self.w]137 }138139 fn grow(&mut self) {140 let cap = self.cap * 2;141 let mut keys = vec![0u64; self.w * cap];142 let mut vals = vec![0i64; cap];143 for s in 0..self.cap {144 if self.empty(s) {145 continue;146 }147 let mut t = (hashk(self.key(s)) as usize) & (cap - 1);148 while keys[t * self.w..(t + 1) * self.w].iter().any(|&x| x != 0) {149 t = (t + 1) & (cap - 1);150 }151 keys[t * self.w..(t + 1) * self.w].copy_from_slice(self.key(s));152 vals[t] = self.vals[s];153 }154 self.cap = cap;155 self.keys = keys;156 self.vals = vals;157 }158159 fn add(&mut self, k: &[u64], v: i64) {160 if (self.len + 1) * 10 > self.cap * 7 {161 self.grow();162 }163 let mut s = (hashk(k) as usize) & (self.cap - 1);164 loop {165 if self.empty(s) {166 self.keys[s * self.w..(s + 1) * self.w].copy_from_slice(k);167 self.vals[s] = v;168 self.len += 1;169 return;170 }171 if self.key(s) == k {172 self.vals[s] += v;173 return;174 }175 s = (s + 1) & (self.cap - 1);176 }177 }178}179180// LAYERS181182struct Layer {183 deg: usize,184 pk: Pack,185 w: usize,186 keys: Vec<u64>,187 nus: Vec<i32>,188 starts: Vec<u32>,189 count: u64,190}191192impl Layer {193 fn new(deg: usize) -> Layer {194 let pk = pack_for(deg);195 let w = pk.w;196 let m2 = binom(deg, 2) as usize;197 Layer {198 deg,199 pk,200 w,201 keys: Vec::new(),202 nus: Vec::new(),203 starts: vec![0; (deg + 1) * (m2 + 1) + 1],204 count: 0,205 }206 }207208 fn m2(&self) -> usize {209 binom(self.deg, 2) as usize210 }211212 fn slot(&self, t1: usize, t2: usize) -> usize {213 t1 * (self.m2() + 1) + t2214 }215216 fn bytes(&self) -> usize {217 self.keys.capacity() * 8 + self.nus.capacity() * 4 + self.starts.capacity() * 4218 }219220 fn key(&self, i: usize) -> &[u64] {221 &self.keys[i * self.w..(i + 1) * self.w]222 }223224 fn range(&self, t1: usize, t2: usize) -> (usize, usize) {225 if t1 > self.deg || t2 > self.m2() {226 return (0, 0);227 }228 let s = self.slot(t1, t2);229 (self.starts[s] as usize, self.starts[s + 1] as usize)230 }231232 fn decode(&self, i: usize, out: &mut [u32]) {233 let k = self.key(i);234 for j in 0..=self.deg {235 out[j] = get(k, j, &self.pk);236 }237 }238}239240fn power(l: usize) -> i64 {241 if l == 0 {242 1243 } else if l == 1 {244 -1245 } else {246 0247 }248}249250// LADDER251252struct Stats {253 odd: Vec<u64>,254 cens: Vec<u64>,255 sums: Vec<i64>,256 mx: Vec<u64>,257 coef: Vec<u64>,258 bits: Vec<u32>,259 peak: usize,260 cut: usize,261 rev_bad: u64,262 rev_seen: u64,263}264265fn ladder(top: usize, quiet: bool, budget: usize) -> Stats {266 let t0 = Instant::now();267 let mut layers: Vec<Layer> = Vec::new();268 let mut st = Stats {269 odd: Vec::new(),270 cens: Vec::new(),271 sums: Vec::new(),272 mx: Vec::new(),273 coef: Vec::new(),274 bits: Vec::new(),275 peak: 0,276 cut: top,277 rev_bad: 0,278 rev_seen: 0,279 };280 let mut run: i64 = 0;281 let mut best: u64 = 0;282 let mut total: u64 = 0;283 if !quiet {284 println!(" L odd census M* census M(q^L) max max/2^L log2(max)/L step coef bits heldMB s");285 }286 let mut cut = false;287 for d in 0..=top {288 let mut l0 = Layer::new(d);289 let mut cmax: u64 = 0;290 if d == 0 {291 let mut k = vec![0u64; l0.w];292 put(&mut k, 0, &l0.pk, 1);293 l0.keys.extend_from_slice(&k);294 l0.nus.push(1);295 l0.starts[0] = 0;296 l0.starts[1] = 1;297 l0.count = 1;298 cmax = 1;299 run += 1;300 best = best.max(run.unsigned_abs());301 } else {302 let store = d < top;303 let pk = pack_for(d);304 let w = l0.w;305 let m2 = l0.m2();306 let mut pc = vec![0u32; top + 2];307 let mut shifts: Vec<Vec<u64>> = (0..=top + 1).map(|_| vec![0u64; w]).collect();308 let mut acc = vec![0u64; w];309 if store {310 let guess = layers[d - 1].count as usize * 5 / 2 + 16;311 l0.keys.reserve_exact(guess * w);312 l0.nus.reserve_exact(guess);313 }314 for v1 in 0..=d {315 for v2 in 0..=m2 {316 let slot = v1 * (m2 + 1) + v2;317 l0.starts[slot] = l0.count as u32;318 let mut tab = Table::new(w, 1 << 6);319 for dp in 0..d {320 let j = d - dp;321 for t1 in 0..=dp {322 if v1 < t1 || v1 - t1 > 1 {323 continue;324 }325 let n1 = v1 - t1;326 if j == 1 && n1 != 1 {327 continue;328 }329 for n2 in 0..=1usize {330 if (j == 1 && n2 != 0) || (j == 2 && n2 != 1) {331 continue;332 }333 let off = t1 * n1 + n2;334 if v2 < off {335 continue;336 }337 let (lo, hi) = layers[dp].range(t1, v2 - off);338 for idx in lo..hi {339 let nu = layers[dp].nus[idx] as i64;340 layers[dp].decode(idx, &mut pc);341 for i in 0..=j {342 let s = &mut shifts[i];343 for x in s.iter_mut() {344 *x = 0;345 }346 for k in 0..=dp {347 if pc[k] != 0 {348 put(s, k + i, &pk, pc[k]);349 }350 }351 }352 for x in acc.iter_mut() {353 *x = 0;354 }355 addk(&mut acc, &shifts[0]);356 addk(&mut acc, &shifts[j]);357 if j >= 2 && n1 == 1 {358 addk(&mut acc, &shifts[j - 1]);359 }360 if j >= 3 && n2 == 1 {361 addk(&mut acc, &shifts[j - 2]);362 }363 tab.add(&acc, nu);364 let f = if j >= 3 { j - 3 } else { 0 };365 for g in 1u64..(1u64 << f) {366 let cur = g ^ (g >> 1);367 let prev = (g - 1) ^ ((g - 1) >> 1);368 let dif = cur ^ prev;369 let i = 1 + dif.trailing_zeros() as usize;370 if cur & dif != 0 {371 addk(&mut acc, &shifts[i]);372 } else {373 subk(&mut acc, &shifts[i]);374 }375 tab.add(&acc, nu);376 }377 }378 }379 }380 }381 let mut ord: Vec<u32> = (0..tab.cap as u32)382 .filter(|&s| !tab.empty(s as usize))383 .collect();384 ord.sort_unstable_by(|&x, &y| cmpk(tab.key(x as usize), tab.key(y as usize)));385 let (lo, hi) = layers[d - 1].range(v1, v2);386 let mut sk = vec![0u64; (hi - lo) * w];387 for (n, i) in (lo..hi).enumerate() {388 layers[d - 1].decode(i, &mut pc);389 for c in 0..d {390 if pc[c] != 0 {391 put(&mut sk[n * w..(n + 1) * w], c + 1, &pk, pc[c]);392 }393 }394 }395 let mut a = 0usize;396 let mut c = 0usize;397 while a < ord.len() || c < hi - lo {398 let take_a = if a == ord.len() {399 false400 } else if c == hi - lo {401 true402 } else {403 cmpk(tab.key(ord[a] as usize), &sk[c * w..(c + 1) * w])404 == Ordering::Less405 };406 if take_a {407 let s = ord[a] as usize;408 let nu = -tab.vals[s];409 run += nu;410 best = best.max(run.unsigned_abs());411 assert!(nu >= i32::MIN as i64 && nu <= i32::MAX as i64);412 let k = tab.key(s);413 assert_eq!(get(k, 0, &pk), 1);414 assert_eq!(get(k, d, &pk), 1);415 for i in 0..=d {416 let a = get(k, i, &pk) as u64;417 assert!(a <= binom(d, i));418 cmax = cmax.max(a);419 }420 if store {421 l0.keys.extend_from_slice(k);422 l0.nus.push(nu as i32);423 }424 l0.count += 1;425 a += 1;426 } else {427 run -= layers[d - 1].nus[lo + c] as i64;428 best = best.max(run.unsigned_abs());429 c += 1;430 }431 }432 st.peak = st.peak.max(433 layers.iter().map(|l| l.bytes()).sum::<usize>()434 + l0.bytes()435 + tab.bytes()436 + ord.capacity() * 4437 + sk.len() * 8,438 );439 if st.peak > budget {440 cut = true;441 break;442 }443 }444 if cut {445 break;446 }447 }448 if cut {449 if !quiet {450 println!(451 " CUT inside level {} on the memory budget, {} MB",452 d + 1,453 st.peak / 1_048_576454 );455 }456 break;457 }458 l0.starts[(d + 1) * (m2 + 1)] = l0.count as u32;459 assert!(cmax <= binom(d, d / 2));460 if store {461 let n = l0.count as usize;462 let mut rk = vec![0u64; w];463 for i in 0..n {464 for x in rk.iter_mut() {465 *x = 0;466 }467 for j in 0..=d {468 let a = get(l0.key(i), j, &pk);469 if a != 0 {470 put(&mut rk, d - j, &pk, a);471 }472 }473 let mut lo = 0usize;474 let mut hi = n;475 while lo < hi {476 let mid = (lo + hi) / 2;477 if cmpk(l0.key(mid), &rk) == Ordering::Less {478 lo = mid + 1;479 } else {480 hi = mid;481 }482 }483 if lo >= n484 || cmpk(l0.key(lo), &rk) != Ordering::Equal485 || l0.nus[lo] != l0.nus[i]486 {487 st.rev_bad += 1;488 }489 st.rev_seen += 1;490 }491 }492 }493 total += l0.count;494 let lv = d + 1;495 let sum = run + power(lv);496 let mx = best.max(sum.unsigned_abs());497 st.odd.push(l0.count);498 st.cens.push(total);499 st.sums.push(sum);500 st.mx.push(mx);501 st.coef.push(cmax);502 st.bits.push(l0.pk.bits);503 layers.push(l0);504 if !quiet {505 let step = if d > 0 {506 mx as f64 / st.mx[d - 1] as f64507 } else {508 0.0509 };510 println!(511 " {:>2} {:>11} {:>12} {:>6} {:>9} {:>10.6} {:>12.6} {:>7.4} {:>5} {:>5} {:>6} {:>5.1}",512 lv,513 st.odd[d],514 total,515 sum,516 mx,517 mx as f64 / (1u64 << lv) as f64,518 (mx as f64).log2() / lv as f64,519 step,520 cmax,521 st.bits[d],522 st.peak / 1_048_576,523 t0.elapsed().as_secs_f64()524 );525 }526 st.cut = lv;527 if lv <= 18 {528 assert_eq!(st.mx[d], PIN_MX[d]);529 assert_eq!(st.cens[d], PIN_CENS[d]);530 assert_eq!(st.coef[d], PIN_COEF[d]);531 assert_eq!(st.sums[d], if lv == 1 { 0 } else { -1 });532 if lv == 18 && !quiet {533 println!(" the Python ladder to level 18 is matched term for term");534 }535 }536 }537 st538}539540// REPORT541542fn report(st: &Stats) {543 let n = st.mx.len();544 println!();545 println!(" running maxima {:?}", st.mx);546 println!(" M(q^L) {:?}", st.sums);547 println!(" monoid census {:?}", st.cens);548 println!(" odd census {:?}", st.odd);549 println!(" max coefficient {:?}", st.coef);550 println!();551 for r in 0..2usize {552 let cl: Vec<String> = (1..=n)553 .filter(|l| l % 2 == r)554 .map(|l| format!("{:.6}", st.mx[l - 1] as f64 / (1u64 << l) as f64))555 .collect();556 println!(" L = {} mod 2 max/2^L {}", r, cl.join(", "));557 }558 println!();559 let mut first = 1;560 for l in 2..=n {561 let lo = st.mx[l - 2] as f64 / (1u64 << (l - 1)) as f64;562 let hi = st.mx[l - 1] as f64 / (1u64 << l) as f64;563 if hi <= lo {564 first = l;565 }566 }567 println!(568 " max/2^L falls last at level {} and rises at every step to level {}",569 first, n570 );571 println!();572 let mut lo = f64::MAX;573 let mut hi = f64::MIN;574 for w in [4usize, 6, 8, 10, 12, 14, 16, 18, 20] {575 if n > w {576 let g = (st.mx[n - 1] as f64 / st.mx[n - 1 - w] as f64).powf(1.0 / w as f64);577 lo = lo.min(g);578 hi = hi.max(g);579 println!(580 " geometric mean step over the last {:>2} levels {:.6} log2 {:.6}",581 w,582 g,583 g.log2()584 );585 }586 }587 println!(588 " the window hull of the step at depth {} is [{:.6}, {:.6}], the mass rate is 2",589 n, lo, hi590 );591 println!();592 for w in [4usize, 6, 8] {593 let mut row = Vec::new();594 let first = if n > w + 3 { n - 3 } else { w + 1 };595 for m in first..=n {596 let g = (st.mx[m - 1] as f64 / st.mx[m - 1 - w] as f64).powf(1.0 / w as f64);597 row.push(format!("{:.6}", g));598 }599 println!(600 " the {:>2} level window read at depths {}..{} {}",601 w,602 first,603 n,604 row.join(", ")605 );606 }607 println!();608 println!(609 " the design mass rate is 2, the deepest level is {}, peak {} MB",610 st.cut,611 st.peak / 1_048_576612 );613 println!(614 " every level to {} is complete, the rate below is a reading to that level",615 n616 );617 println!(618 " the reciprocal Q -> x^(deg Q) Q(1/x) holds nu* fixed, {} mismatches over {} elements",619 st.rev_bad, st.rev_seen620 );621}622623// PYTHON PIN624625const PIN_MX: [u64; 18] = [626 1, 1, 2, 3, 4, 7, 15, 23, 45, 86, 162, 331, 741, 1665, 3173, 7508, 17753, 36147,627];628629const PIN_CENS: [u64; 18] = [630 1, 2, 5, 11, 27, 61, 144, 331, 776, 1788, 4147, 9544, 22000, 50420, 115407, 263062, 598540,631 1357535,632];633634const PIN_COEF: [u64; 18] = [635 1, 1, 2, 3, 6, 10, 20, 35, 70, 126, 252, 462, 924, 1716, 3432, 6435, 12870, 24310,636];637638// VERBS639640fn pin() {641 let st = ladder(17, false, usize::MAX);642 report(&st);643 println!();644 let mut bad = 0;645 for l in 1..=18usize {646 if st.mx[l - 1] != PIN_MX[l - 1] {647 bad += 1;648 }649 if st.cens[l - 1] != PIN_CENS[l - 1] {650 bad += 1;651 }652 if st.coef[l - 1] != PIN_COEF[l - 1] {653 bad += 1;654 }655 let want = if l == 1 { 0 } else { -1 };656 if st.sums[l - 1] != want {657 bad += 1;658 }659 }660 println!(" the Python ladder to level 18, mismatches {}", bad);661}662663fn main() {664 let a: Vec<String> = env::args().collect();665 let verb = a.get(1).map(|s| s.as_str()).unwrap_or("pin");666 match verb {667 "pin" => pin(),668 "ladder" => {669 let top: usize = a.get(2).and_then(|s| s.parse().ok()).unwrap_or(21);670 let gb: f64 = a.get(3).and_then(|s| s.parse().ok()).unwrap_or(2.7);671 let st = ladder(top, false, (gb * 1_073_741_824.0) as usize);672 report(&st);673 }674 _ => println!("verbs: pin, ladder <top degree> <GB>"),675 }676}677678// TESTS679680#[cfg(test)]681mod tests {682 use super::*;683684 #[test]685 fn ladder_pins_to_the_python_at_level_eighteen() {686 let st = ladder(17, true, usize::MAX);687 assert_eq!(st.mx, PIN_MX.to_vec());688 assert_eq!(st.cens, PIN_CENS.to_vec());689 assert_eq!(st.coef, PIN_COEF.to_vec());690 assert_eq!(st.sums[0], 0);691 assert!(st.sums[1..].iter().all(|&s| s == -1));692 }693}