menergy.rs
25.4 kB · rust · 910 lines
1// RANDOM23pub struct Rng(u64);45impl Rng {6 pub fn new(seed: u64) -> Rng {7 Rng(seed ^ 0x9e3779b97f4a7c15)8 }910 fn step(&mut self) -> u64 {11 self.0 = self.0.wrapping_add(0x9e3779b97f4a7c15);12 let mut z = self.0;13 z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);14 z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);15 z ^ (z >> 31)16 }1718 fn wide(&mut self) -> u128 {19 ((self.step() as u128) << 64) | self.step() as u12820 }2122 pub fn sign(&mut self) -> f64 {23 if self.step() & 1 == 0 {24 1.025 } else {26 -1.027 }28 }29}3031// COLUMNS3233pub fn column(q: u64, digits: &[u64], l: usize) -> Vec<u128> {34 let mut out = vec![0u128];35 for _ in 0..l {36 let mut next = Vec::with_capacity(out.len() * digits.len());37 for &v in &out {38 for &f in digits {39 next.push(v * q as u128 + f as u128);40 }41 }42 out = next;43 }44 out.retain(|&v| v != 0);45 out.sort_unstable();46 out47}4849pub fn random_column(x: u128, count: usize, rng: &mut Rng) -> Vec<u128> {50 let mut out: Vec<u128> = Vec::with_capacity(count);51 while out.len() < count {52 for _ in 0..count - out.len() {53 out.push(1 + rng.wide() % (x - 1));54 }55 out.sort_unstable();56 out.dedup();57 }58 out59}6061// GLOBAL ENERGY6263pub trait Prod: Copy + Ord {64 fn times(self, other: Self) -> Self;65 fn scramble(self) -> u64;66}6768impl Prod for u64 {69 fn times(self, other: u64) -> u64 {70 self * other71 }7273 fn scramble(self) -> u64 {74 (self ^ (self >> 33)).wrapping_mul(0xff51afd7ed558ccd) >> 2475 }76}7778impl Prod for u128 {79 fn times(self, other: u128) -> u128 {80 self * other81 }8283 fn scramble(self) -> u64 {84 let v = (self as u64) ^ ((self >> 64) as u64);85 (v ^ (v >> 33)).wrapping_mul(0xff51afd7ed558ccd) >> 2486 }87}8889pub fn energy<T: Prod>(vals: &[T], cap: usize) -> (u128, u32) {90 let n = vals.len() as u128;91 let total = n * n;92 let mut parts = 1usize;93 while total / parts as u128 > cap as u128 {94 parts <<= 1;95 }96 let mask = (parts - 1) as u64;97 let mut e: u128 = 0;98 let mut top: u32 = 0;99 let mut buf: Vec<T> = Vec::with_capacity((total / parts as u128) as usize * 5 / 4 + 64);100 for p in 0..parts as u64 {101 buf.clear();102 for &a in vals {103 for &b in vals {104 let m = a.times(b);105 if parts == 1 || m.scramble() & mask == p {106 buf.push(m);107 }108 }109 }110 buf.sort_unstable();111 let mut i = 0;112 while i < buf.len() {113 let mut j = i + 1;114 while j < buf.len() && buf[j] == buf[i] {115 j += 1;116 }117 let r = (j - i) as u128;118 e += r * r;119 top = top.max((j - i) as u32);120 i = j;121 }122 }123 (e, top)124}125126pub fn energy_of(vals: &[u128], cap: usize) -> (u128, u32) {127 if *vals.last().unwrap() < 1u128 << 32 {128 let small: Vec<u64> = vals.iter().map(|&v| v as u64).collect();129 energy(&small, cap)130 } else {131 energy(vals, cap)132 }133}134135// SHIFT SOLUTIONS136137fn live(k: u128, m: i64) -> u128 {138 if m >= 1 {139 let mut p = k - 1;140 for _ in 1..m {141 p *= k;142 }143 p144 } else {145 0146 }147}148149pub fn shift_excess(k: u128, l: usize) -> u128 {150 let mut t = 0u128;151 for s in 0..2 * l {152 for i in 0..=s {153 for ip in 0..=s {154 if i == ip {155 continue;156 }157 let a = live(k, l as i64 - i.max(ip) as i64);158 let b = live(k, l as i64 - (s - i).max(s - ip) as i64);159 if a == 0 || b == 0 {160 continue;161 }162 t += a * b - a.min(b);163 }164 }165 }166 t167}168169// MEMBERSHIP170171pub struct Bits {172 w: Vec<u64>,173}174175impl Bits {176 fn new(x: usize) -> Bits {177 Bits {178 w: vec![0u64; (x >> 6) + 1],179 }180 }181182 fn set(&mut self, n: usize) {183 self.w[n >> 6] |= 1u64 << (n & 63);184 }185186 pub fn get(&self, n: usize) -> bool {187 self.w[n >> 6] >> (n & 63) & 1 == 1188 }189}190191pub fn bits_of(vals: &[u128], x: u128) -> Bits {192 let mut b = Bits::new(x as usize);193 for &v in vals {194 b.set(v as usize);195 }196 b197}198199// BOXES200201pub struct BoxRow {202 pub m: u64,203 pub n: u64,204 pub r: u128,205 pub e: u128,206 pub ebal: f64,207 pub diag: f64,208 pub triv: f64,209 pub bound: f64,210 pub check: (f64, f64),211}212213pub fn box_row(bits: &Bits, mm: u64, nn: u64, delta: f64, cap: usize) -> Option<BoxRow> {214 let mut cs = vec![0u32; nn as usize];215 let mut r: u128 = 0;216 let mut e: u128 = 0;217 let mut spread = 0.0f64;218 for m in mm..2 * mm {219 let mut c = 0u32;220 for l in nn..2 * nn {221 if bits.get((m * l) as usize) {222 c += 1;223 cs[(l - nn) as usize] += 1;224 }225 }226 r += c as u128;227 e += (c as u128) * (c as u128);228 let d = c as f64 - nn as f64 * delta;229 spread += d * d;230 }231 if e > cap as u128 {232 return None;233 }234 let mut keys: Vec<u64> = Vec::with_capacity(e as usize + 8);235 let mut hits: Vec<u32> = Vec::with_capacity(nn as usize);236 for m in mm..2 * mm {237 hits.clear();238 for l in nn..2 * nn {239 if bits.get((m * l) as usize) {240 hits.push((l - nn) as u32);241 }242 }243 for &a in &hits {244 for &b in &hits {245 keys.push(a as u64 * nn + b as u64);246 }247 }248 }249 keys.sort_unstable();250 let base = mm as f64 * delta * delta;251 let mut sorted = cs.clone();252 sorted.sort_unstable();253 let mut dv: Vec<(f64, f64)> = Vec::new();254 let mut i = 0;255 while i < sorted.len() {256 let mut j = i + 1;257 while j < sorted.len() && sorted[j] == sorted[i] {258 j += 1;259 }260 dv.push((sorted[i] as f64, (j - i) as f64));261 i = j;262 }263 let mut ebal = 0.0f64;264 let mut signed = 0.0f64;265 for &(v1, n1) in &dv {266 for &(v2, n2) in &dv {267 ebal += n1 * n2 * (base - delta * (v1 + v2)).abs();268 signed += n1 * n2 * (base - delta * (v1 + v2));269 }270 }271 signed += e as f64;272 let mut diag = 0.0f64;273 for &c in &cs {274 diag += (c as f64 * (1.0 - 2.0 * delta) + base).abs();275 }276 let mut i = 0;277 while i < keys.len() {278 let mut j = i + 1;279 while j < keys.len() && keys[j] == keys[i] {280 j += 1;281 }282 let t = (j - i) as f64;283 let l1 = (keys[i] / nn) as usize;284 let l2 = (keys[i] % nn) as usize;285 let z = base - delta * (cs[l1] as f64 + cs[l2] as f64);286 ebal += (t + z).abs() - z.abs();287 i = j;288 }289 let cells = mm as f64 * nn as f64;290 Some(BoxRow {291 m: mm,292 n: nn,293 r,294 e,295 ebal,296 diag,297 triv: r as f64 * (1.0 - delta) + (cells - r as f64) * delta,298 bound: (mm as f64 * ebal).sqrt(),299 check: (signed, spread),300 })301}302303pub fn boxes(x: u128) -> Vec<(u64, u64)> {304 let mut out = Vec::new();305 let mut mm = 8u64;306 while (mm as u128) * 32 <= x {307 let mut nn = 8u64;308 while (nn as u128) * 4 * mm as u128 <= x {309 if (nn as u128) * 32 * mm as u128 >= x {310 out.push((mm, nn));311 }312 nn <<= 1;313 }314 mm <<= 1;315 }316 out317}318319pub struct Sweep {320 pub worst: BoxRow,321 pub best: BoxRow,322 pub l5: Option<BoxRow>,323 pub seen: usize,324 pub live: usize,325}326327fn sweep(bits: &Bits, x: u128, delta: f64, cap: usize) -> Option<Sweep> {328 let mut worst: Option<BoxRow> = None;329 let mut best: Option<BoxRow> = None;330 let all = boxes(x);331 let mut live = 0;332 let mut l5: Option<BoxRow> = None;333 let edge = (x as f64).powf(0.4);334 for &(mm, nn) in all.iter() {335 let row = match box_row(bits, mm, nn, delta, cap) {336 Some(r) => r,337 None => continue,338 };339 if row.r == 0 {340 continue;341 }342 live += 1;343 if worst.as_ref().map(|b| row.bound > b.bound).unwrap_or(true) {344 worst = Some(row_copy(&row));345 }346 if mm as f64 >= edge347 && nn as f64 >= edge348 && l5.as_ref().map(|b| row.bound > b.bound).unwrap_or(true)349 {350 l5 = Some(row_copy(&row));351 }352 if best.as_ref().map(|b| row.bound < b.bound).unwrap_or(true) {353 best = Some(row);354 }355 }356 match (worst, best) {357 (Some(w), Some(b)) => Some(Sweep {358 worst: w,359 best: b,360 l5,361 seen: all.len(),362 live,363 }),364 _ => None,365 }366}367368fn row_copy(r: &BoxRow) -> BoxRow {369 BoxRow {370 m: r.m,371 n: r.n,372 r: r.r,373 e: r.e,374 ebal: r.ebal,375 diag: r.diag,376 triv: r.triv,377 bound: r.bound,378 check: r.check,379 }380}381382// SIGMA383384fn sigma_max(bits: &Bits, mm: u64, nn: u64, delta: f64, rounds: usize, rng: &mut Rng) -> f64 {385 let mut worst = 0.0f64;386 let mut asign = vec![0.0f64; mm as usize];387 let mut bsign = vec![0.0f64; nn as usize];388 for _ in 0..rounds {389 for a in asign.iter_mut() {390 *a = rng.sign();391 }392 for b in bsign.iter_mut() {393 *b = rng.sign();394 }395 let mut hit = 0.0f64;396 for m in mm..2 * mm {397 let am = asign[(m - mm) as usize];398 for l in nn..2 * nn {399 if bits.get((m * l) as usize) {400 hit += am * bsign[(l - nn) as usize];401 }402 }403 }404 let asum: f64 = asign.iter().sum();405 let bsum: f64 = bsign.iter().sum();406 worst = worst.max((hit - delta * asum * bsum).abs());407 }408 worst409}410411// STUDY412413struct Cell {414 q: u64,415 digits: Vec<u64>,416 label: &'static str,417 lmax: usize,418 lrand: usize,419 typeii: Vec<usize>,420}421422fn ex(q: u64, e: u64) -> Vec<u64> {423 (0..q).filter(|&f| f != e).collect()424}425426fn cells() -> Vec<Cell> {427 vec![428 Cell {429 q: 3,430 digits: vec![0, 1],431 label: "01",432 lmax: 12,433 lrand: 12,434 typeii: vec![10, 12, 14],435 },436 Cell {437 q: 3,438 digits: vec![0, 2],439 label: "02",440 lmax: 10,441 lrand: 99,442 typeii: vec![],443 },444 Cell {445 q: 3,446 digits: vec![1, 2],447 label: "12",448 lmax: 12,449 lrand: 12,450 typeii: vec![12],451 },452 Cell {453 q: 4,454 digits: vec![0, 1, 2],455 label: "012",456 lmax: 8,457 lrand: 8,458 typeii: vec![9],459 },460 Cell {461 q: 5,462 digits: vec![0, 1, 2, 3],463 label: "0123",464 lmax: 6,465 lrand: 6,466 typeii: vec![8],467 },468 Cell {469 q: 5,470 digits: vec![0, 2, 4],471 label: "024",472 lmax: 7,473 lrand: 7,474 typeii: vec![7],475 },476 Cell {477 q: 10,478 digits: ex(10, 7),479 label: "ex7",480 lmax: 4,481 lrand: 4,482 typeii: vec![6],483 },484 Cell {485 q: 10,486 digits: ex(10, 0),487 label: "ex0",488 lmax: 3,489 lrand: 3,490 typeii: vec![4],491 },492 Cell {493 q: 100,494 digits: (0..50).collect(),495 label: "0to49",496 lmax: 2,497 lrand: 2,498 typeii: vec![3],499 },500 Cell {501 q: 100,502 digits: vec![0, 1],503 label: "01",504 lmax: 9,505 lrand: 9,506 typeii: vec![3],507 },508 Cell {509 q: 100,510 digits: ex(100, 37),511 label: "ex37",512 lmax: 1,513 lrand: 1,514 typeii: vec![2],515 },516 ]517}518519pub fn qpow(q: u64, l: usize) -> u128 {520 let mut p = 1u128;521 for _ in 0..l {522 p *= q as u128;523 }524 p525}526527pub fn exps(v: f64, x: u128) -> f64 {528 if v <= 0.0 {529 0.0530 } else {531 v.ln() / (x as f64).ln()532 }533}534535const CAP: usize = 1 << 24;536const PAIRCAP: usize = 8_000_000;537538pub fn run() {539 let cs = cells();540 println!("menergy census");541 println!("| q | F | L | K | E_x | theta_x | 2 alpha | E_x/(2K^2 - K) | max r | shift T | shift share | excess share | E_rand | E_rand/(2K^2 - K) | E_x/E_rand |");542 for c in cs.iter() {543 let alpha = (c.digits.len() as f64).ln() / (c.q as f64).ln();544 for l in 1..=c.lmax {545 let x = qpow(c.q, l);546 assert!(x <= u128::MAX / x, "q^(2L) overflows u128");547 let vals = column(c.q, &c.digits, l);548 let k = vals.len() as u128;549 let (e, top) = energy_of(&vals, CAP);550 assert!(e >= 2 * k * k - k, "diagonal floor fails");551 assert!(e <= k * k * k, "trivial ceiling fails");552 assert!(553 e as f64 <= (top as f64) * (k * k) as f64 + 0.5,554 "divisor ceiling fails"555 );556 let diag = (2 * k * k - k) as f64;557 let shift = if c.digits.contains(&0) {558 shift_excess(c.digits.len() as u128, l)559 } else {560 0561 };562 assert!(e >= 2 * k * k - k + shift, "shift floor fails");563 let (rand, randratio, overrand) = if l >= c.lrand {564 let mut rng = Rng::new(0x5eed + l as u64 + c.q * 977);565 let rv = random_column(x, k as usize, &mut rng);566 let (re, _) = energy_of(&rv, CAP);567 (568 format!("{re}"),569 format!("{:.4}", re as f64 / diag),570 format!("{:.4}", e as f64 / re as f64),571 )572 } else {573 ("-".to_string(), "-".to_string(), "-".to_string())574 };575 let excess = e as f64 - diag;576 let exshare = if excess > 0.0 {577 format!("{:.4}", shift as f64 / excess)578 } else {579 "-".to_string()580 };581 println!(582 "| {} | {} | {} | {} | {} | {:.6} | {:.6} | {:.4} | {} | {} | {:.4} | {} | {} | {} | {} |",583 c.q,584 c.label,585 l,586 k,587 e,588 exps(e as f64, x),589 2.0 * alpha,590 e as f64 / diag,591 top,592 shift,593 (diag + shift as f64) / e as f64,594 exshare,595 rand,596 randratio,597 overrand598 );599 }600 }601 println!("menergy type II");602 println!("| q | F | L | x | alpha | kind | box | boxes | M | N | R | E_x(M,N) | E_bal | diag share | bound | bound/triv | bound exp | alpha - exp |");603 for c in cs.iter() {604 let alpha = (c.digits.len() as f64).ln() / (c.q as f64).ln();605 for &l in c.typeii.iter() {606 let x = qpow(c.q, l);607 let vals = column(c.q, &c.digits, l);608 let delta = vals.len() as f64 / x as f64;609 let mut rng = Rng::new(0xb0a7 + l as u64 + c.q * 131);610 let rv = random_column(x, vals.len(), &mut rng);611 for (kind, set) in [("digit", &vals), ("random", &rv)] {612 let bits = bits_of(set, x);613 let sw = match sweep(&bits, x, delta, PAIRCAP) {614 Some(s) => s,615 None => {616 println!(617 "| {} | {} | {} | {} | {:.6} | {} | none | 0/{} | - | - | - | - | - | - | - | - | - | - |",618 c.q, c.label, l, x, alpha, kind, boxes(x).len()619 );620 continue;621 }622 };623 let mut shown: Vec<(&str, &BoxRow)> = vec![("worst", &sw.worst)];624 if let Some(r) = sw.l5.as_ref() {625 shown.push(("l5", r));626 }627 shown.push(("best", &sw.best));628 for (tag, row) in shown {629 let expo = exps(row.bound, x);630 println!(631 "| {} | {} | {} | {} | {:.6} | {} | {} | {}/{} | {} | {} | {} | {} | {:.4e} | {:.4} | {:.4e} | {:.4} | {:.6} | {:.6} |",632 c.q,633 c.label,634 l,635 x,636 alpha,637 kind,638 tag,639 sw.live,640 sw.seen,641 row.m,642 row.n,643 row.r,644 row.e,645 row.ebal,646 row.diag / row.ebal,647 row.bound,648 row.bound / row.triv,649 expo,650 alpha - expo651 );652 }653 }654 }655 }656 println!("menergy sigma");657 println!("| q | F | L | M | N | rounds | max Sigma | bound | ratio | signed check |");658 for (q, digits, label, l, mm, nn) in [659 (3u64, vec![0u64, 1], "01", 6usize, 8u64, 8u64),660 (3, vec![0, 1], "01", 8, 8, 32),661 (3, vec![0, 1], "01", 8, 16, 16),662 (3, vec![1, 2], "12", 8, 16, 16),663 (4, vec![0, 1, 2], "012", 6, 16, 32),664 (5, vec![0, 1, 2, 3], "0123", 5, 16, 32),665 ] {666 let x = qpow(q, l);667 let vals = column(q, &digits, l);668 let delta = vals.len() as f64 / x as f64;669 let bits = bits_of(&vals, x);670 let row = box_row(&bits, mm, nn, delta, PAIRCAP).unwrap();671 let mut rng = Rng::new(0x5169 + q * 7 + l as u64);672 let worst = sigma_max(&bits, mm, nn, delta, 40, &mut rng);673 assert!(674 worst <= row.bound * (1.0 + 1e-9),675 "sigma exceeds the energy bound at q={q} L={l}"676 );677 println!(678 "| {} | {} | {} | {} | {} | {} | {:.4e} | {:.4e} | {:.4} | {:.3e} |",679 q,680 label,681 l,682 mm,683 nn,684 40,685 worst,686 row.bound,687 worst / row.bound,688 (row.check.0 - row.check.1).abs() / row.check.1.max(1.0)689 );690 }691}692693// TESTS694695#[cfg(test)]696mod tests {697 use super::*;698699 fn brute_energy(vals: &[u128]) -> u128 {700 let mut e = 0u128;701 for &a in vals {702 for &b in vals {703 for &c in vals {704 for &d in vals {705 if a * b == c * d {706 e += 1;707 }708 }709 }710 }711 }712 e713 }714715 #[test]716 fn energy_matches_brute() {717 for (q, digits, l) in [718 (3u64, vec![0u64, 1], 4usize),719 (3, vec![1, 2], 4),720 (4, vec![0, 1, 2], 3),721 (5, vec![0, 1, 2, 3], 2),722 (10, ex(10, 7), 2),723 (100, vec![0, 1], 3),724 ] {725 let vals = column(q, &digits, l);726 let (e, _) = energy_of(&vals, CAP);727 assert_eq!(e, brute_energy(&vals), "energy q={q} L={l}");728 }729 }730731 #[test]732 fn partition_invariance() {733 let vals = column(3, &[0, 1], 8);734 let (a, ta) = energy_of(&vals, CAP);735 let (b, tb) = energy_of(&vals, 64);736 let (c, tc) = energy(&vals, 7);737 assert_eq!((a, ta), (b, tb));738 assert_eq!((a, ta), (c, tc));739 }740741 #[test]742 fn diagonal_and_ceiling() {743 for (q, digits, l) in [744 (3u64, vec![0u64, 1], 9usize),745 (4, vec![0, 1, 2], 6),746 (10, ex(10, 0), 3),747 (100, (0..50).collect::<Vec<u64>>(), 2),748 ] {749 let vals = column(q, &digits, l);750 let k = vals.len() as u128;751 let (e, top) = energy_of(&vals, CAP);752 assert!(e >= 2 * k * k - k);753 assert!(e <= k * k * k);754 assert!(e <= top as u128 * k * k);755 }756 }757758 #[test]759 fn scaling_invariance() {760 for l in 1..=9 {761 let a = energy_of(&column(3, &[0, 1], l), CAP);762 let b = energy_of(&column(3, &[0, 2], l), CAP);763 assert_eq!(a, b, "scaling q=3 L={l}");764 }765 for l in 1..=5 {766 let a = energy_of(&column(5, &[0, 1, 2], l), CAP);767 let b = energy_of(&column(5, &[0, 2, 4], l), CAP);768 assert_eq!(a, b, "scaling q=5 L={l}");769 }770 let a = energy_of(&column(3, &[0, 1], 6), CAP);771 let b = energy_of(&column(3, &[1, 2], 6), CAP);772 assert!(a != b, "translation is not a multiplicative invariance");773 }774775 #[test]776 fn restricted_energy_matches_brute() {777 let (q, digits, l, mm, nn) = (3u64, vec![0u64, 1], 8usize, 8u64, 16u64);778 let x = qpow(q, l);779 let vals = column(q, &digits, l);780 let delta = vals.len() as f64 / x as f64;781 let bits = bits_of(&vals, x);782 let row = box_row(&bits, mm, nn, delta, PAIRCAP).unwrap();783 let mut r = 0u128;784 let mut e = 0u128;785 for m in mm..2 * mm {786 for l1 in nn..2 * nn {787 if bits.get((m * l1) as usize) {788 r += 1;789 }790 for l2 in nn..2 * nn {791 if bits.get((m * l1) as usize) && bits.get((m * l2) as usize) {792 e += 1;793 }794 }795 }796 }797 assert_eq!(row.r, r);798 assert_eq!(row.e, e);799 let mut ebal = 0.0f64;800 for l1 in nn..2 * nn {801 for l2 in nn..2 * nn {802 let mut w = 0.0f64;803 for m in mm..2 * mm {804 let p1 = if bits.get((m * l1) as usize) {805 1.0806 } else {807 0.0808 } - delta;809 let p2 = if bits.get((m * l2) as usize) {810 1.0811 } else {812 0.0813 } - delta;814 w += p1 * p2;815 }816 ebal += w.abs();817 }818 }819 assert!(820 (row.ebal - ebal).abs() < 1e-6 * ebal,821 "E_bal {} against {ebal}",822 row.ebal823 );824 assert!((row.check.0 - row.check.1).abs() < 1e-6 * row.check.1);825 }826827 #[test]828 fn sigma_stays_under_the_bound() {829 for (q, digits, l, mm, nn) in [830 (3u64, vec![0u64, 1], 8usize, 8u64, 16u64),831 (3, vec![1, 2], 8, 16, 16),832 (4, vec![0, 1, 2], 6, 8, 32),833 ] {834 let x = qpow(q, l);835 let vals = column(q, &digits, l);836 let delta = vals.len() as f64 / x as f64;837 let bits = bits_of(&vals, x);838 let row = box_row(&bits, mm, nn, delta, PAIRCAP).unwrap();839 let mut rng = Rng::new(11 + q + l as u64);840 let worst = sigma_max(&bits, mm, nn, delta, 60, &mut rng);841 assert!(842 worst <= row.bound,843 "sigma {worst} over bound {} at q={q}",844 row.bound845 );846 }847 }848849 #[test]850 fn shift_solutions_are_solutions() {851 for (q, digits, l) in [852 (3u64, vec![0u64, 1], 6usize),853 (4, vec![0, 1, 2], 4),854 (5, vec![0, 2, 4], 4),855 (10, ex(10, 7), 3),856 ] {857 let k = digits.len() as u128;858 let vals = column(q, &digits, l);859 let (e, _) = energy_of(&vals, CAP);860 let kk = vals.len() as u128;861 let t = shift_excess(k, l);862 assert!(t > 0);863 assert!(e >= 2 * kk * kk - kk + t, "shift floor q={q} L={l}");864 }865 let mut seen = std::collections::HashSet::new();866 let (q, digits, l) = (3u64, vec![0u64, 1], 5usize);867 let vals = column(q, &digits, l);868 let set: std::collections::HashSet<u128> = vals.iter().copied().collect();869 let mut count = 0u128;870 for s in 0..2 * l {871 for i in 0..=s {872 for ip in 0..=s {873 if i == ip {874 continue;875 }876 for &u in vals.iter() {877 for &v in vals.iter() {878 if u == v || u % q as u128 == 0 || v % q as u128 == 0 {879 continue;880 }881 let p = |e: usize| (q as u128).pow(e as u32);882 let quad = (p(i) * u, p(s - i) * v, p(ip) * u, p(s - ip) * v);883 if !set.contains(&quad.0)884 || !set.contains(&quad.1)885 || !set.contains(&quad.2)886 || !set.contains(&quad.3)887 {888 continue;889 }890 assert_eq!(quad.0 * quad.1, quad.2 * quad.3);891 if seen.insert(quad) {892 count += 1;893 }894 }895 }896 }897 }898 }899 assert_eq!(count, shift_excess(2, l), "shift count L={l}");900 }901902 #[test]903 fn random_column_is_a_set() {904 let mut rng = Rng::new(4);905 let v = random_column(1000, 400, &mut rng);906 assert_eq!(v.len(), 400);907 assert!(v.windows(2).all(|w| w[0] < w[1]));908 assert!(v[0] >= 1 && *v.last().unwrap() < 1000);909 }910}