main.rs
15.2 kB · rust · 512 lines
1mod census;2mod factor;34use census::{digit_gcd, families, run_family, sweep, Family, Outcome};5use factor::{mobius, mu_sieve, small_primes};6use std::sync::atomic::{AtomicUsize, Ordering};7use std::sync::Mutex;89const SIEVE_LIMIT: usize = 129_140_163;10const KEMPNER_LIMIT: u64 = 100_000_000;1112// READOUT1314fn theta(m: i64, a: u64) -> Option<f64> {15 if m == 0 || a < 2 {16 return None;17 }18 Some((m.unsigned_abs() as f64).ln() / (a as f64).ln())19}2021fn theta_str(m: i64, a: u64) -> String {22 match theta(m, a) {23 Some(t) => format!("{t:.4}"),24 None => "-".to_string(),25 }26}2728fn row(tag: &str, q: u64, label: &str, lev: usize, a: u64, m: i64, mx: u64) {29 println!(30 "{tag} q={q} F={label} l={lev} A={a} M={m} theta={} Mmax={mx} thetamax={}",31 theta_str(m, a),32 theta_str(mx as i64, a)33 );34}3536fn drift(meter: &[u64], counts: &[u64], last: usize) -> String {37 let l = meter.len() - 1;38 let lo = if l > last { l - last + 1 } else { 1 };39 let mut values = Vec::new();40 for lev in lo..=l {41 if let Some(t) = theta(meter[lev] as i64, counts[lev]) {42 values.push(t);43 }44 }45 if values.len() < 2 {46 return "-".to_string();47 }48 let max = values.iter().cloned().fold(f64::MIN, f64::max);49 let min = values.iter().cloned().fold(f64::MAX, f64::min);50 format!("{:.4}", max - min)51}5253// CONTROLS5455struct ControlColumn {56 q: u64,57 levels: usize,58 meter: Vec<i64>,59 mmax: Vec<u64>,60}6162fn controls_and_kempner(63 mu: &[i8],64) -> (65 Vec<ControlColumn>,66 Vec<Vec<i64>>,67 Vec<Vec<u64>>,68 Vec<Vec<u64>>,69) {70 let mut ctrl: Vec<ControlColumn> = [(3u64, 17usize), (4, 13), (5, 11), (10, 8)]71 .iter()72 .map(|&(q, levels)| ControlColumn {73 q,74 levels,75 meter: vec![0i64; levels + 1],76 mmax: vec![0u64; levels + 1],77 })78 .collect();79 let mut marks: Vec<(u64, usize, usize)> = Vec::new();80 for (ci, c) in ctrl.iter().enumerate() {81 let mut p = 1u64;82 for lev in 1..=c.levels {83 p *= c.q;84 marks.push((p, ci, lev));85 }86 }87 marks.sort();88 let mut kem_m = vec![vec![0i64; 9]; 10];89 let mut kem_a = vec![vec![0u64; 9]; 10];90 let mut kem_x = vec![vec![0u64; 9]; 10];91 let mut kem_run = [0i64; 10];92 let mut kem_cnt = [0u64; 10];93 let mut kem_peak = [0u64; 10];94 let mut run = 0i64;95 let mut peak = 0u64;96 let mut next = 0usize;97 for n in 1..=SIEVE_LIMIT as u64 {98 let m = mu[n as usize] as i64;99 run += m;100 peak = peak.max(run.unsigned_abs());101 if n <= KEMPNER_LIMIT {102 let mut mask = 0u16;103 let mut x = n;104 while x > 0 {105 mask |= 1 << (x % 10);106 x /= 10;107 }108 for d in 0..10 {109 if mask >> d & 1 == 0 {110 kem_run[d] += m;111 kem_cnt[d] += 1;112 kem_peak[d] = kem_peak[d].max(kem_run[d].unsigned_abs());113 }114 }115 }116 while next < marks.len() && marks[next].0 == n {117 let (_, ci, lev) = marks[next];118 ctrl[ci].meter[lev] = run;119 ctrl[ci].mmax[lev] = peak;120 if ctrl[ci].q == 10 {121 for d in 0..10 {122 kem_m[d][lev] = kem_run[d];123 kem_a[d][lev] = kem_cnt[d];124 kem_x[d][lev] = kem_peak[d];125 }126 }127 next += 1;128 }129 }130 assert_eq!(next, marks.len());131 (ctrl, kem_m, kem_a, kem_x)132}133134// CROSSCHECK135136fn crosscheck(mu: &[i8], primes: &[u64]) {137 let lmax = 16usize;138 let digits = [1u64, 2];139 let mut by_factor = vec![0i64; lmax + 1];140 let mut by_sieve = vec![0i64; lmax + 1];141 sweep(3, &digits, lmax, &mut |v, len| {142 by_factor[len] += mobius(v, primes) as i64;143 by_sieve[len] += mu[v as usize] as i64;144 });145 let mut cf = 0i64;146 let mut cs = 0i64;147 for lev in 1..=lmax {148 cf += by_factor[lev];149 cs += by_sieve[lev];150 assert_eq!(cf, cs, "method mismatch at level {lev}");151 }152 println!("crosscheck q=3 F=12 L=16 factorization equals sieve at every level, M(3^16)={cf}");153}154155// MAIN156157fn main() {158 println!("mobius-designs generator: CARGO_BUILD_JOBS=4 cargo run --release -p mobius-designs");159 let primes = small_primes(1024);160 let mu = mu_sieve(SIEVE_LIMIT);161 let (ctrl, kem_m, kem_a, kem_x) = controls_and_kempner(&mu);162 for c in &ctrl {163 let mut p = 1u64;164 for lev in 1..=c.levels {165 p *= c.q;166 row("control", c.q, "full", lev, p, c.meter[lev], c.mmax[lev]);167 }168 }169 for d in 0..10 {170 for lev in 1..=8usize {171 row(172 "kempner",173 10,174 &format!("X{d}"),175 lev,176 kem_a[d][lev],177 kem_m[d][lev],178 kem_x[d][lev],179 );180 }181 }182 crosscheck(&mu, &primes);183 drop(mu);184 let jobs = families();185 let results: Vec<Mutex<Option<Outcome>>> = jobs.iter().map(|_| Mutex::new(None)).collect();186 let cursor = AtomicUsize::new(0);187 std::thread::scope(|s| {188 for _ in 0..6 {189 s.spawn(|| loop {190 let i = cursor.fetch_add(1, Ordering::Relaxed);191 if i >= jobs.len() {192 break;193 }194 let out = run_family(&jobs[i], &primes);195 *results[i].lock().unwrap() = Some(out);196 });197 }198 });199 let outcomes: Vec<Outcome> = results200 .into_iter()201 .map(|r| r.into_inner().unwrap().unwrap())202 .collect();203 for (fam, out) in jobs.iter().zip(outcomes.iter()) {204 for lev in 1..=fam.lmax {205 row(206 "row",207 fam.q,208 &fam.label,209 lev,210 out.counts[lev],211 out.meter[lev],212 out.mmax[lev],213 );214 }215 }216 verify_scalings(&jobs, &outcomes);217 for (fam, out) in jobs.iter().zip(outcomes.iter()) {218 let l = fam.lmax;219 println!(220 "slope q={} F={} L={l} theta={} thetamax={} driftmax5={}",221 fam.q,222 fam.label,223 theta_str(out.meter[l], out.counts[l]),224 theta_str(out.mmax[l] as i64, out.counts[l]),225 drift(&out.mmax, &out.counts, 5)226 );227 }228 for q in [3u64, 4, 5] {229 let mut entries: Vec<(String, Option<f64>)> = jobs230 .iter()231 .zip(outcomes.iter())232 .filter(|(f, _)| f.q == q)233 .map(|(f, o)| {234 (235 f.label.clone(),236 theta(o.mmax[f.lmax] as i64, o.counts[f.lmax]),237 )238 })239 .collect();240 entries.sort_by(|a, b| {241 a.1.unwrap_or(f64::MIN)242 .partial_cmp(&b.1.unwrap_or(f64::MIN))243 .unwrap()244 });245 let text: Vec<String> = entries246 .iter()247 .map(|(l, t)| match t {248 Some(t) => format!("F={l} {t:.4}"),249 None => format!("F={l} -"),250 })251 .collect();252 println!(253 "distribution q={q} final-level thetamax ascending: {}",254 text.join(" | ")255 );256 }257 let mut kem: Vec<(usize, f64)> = (0..10)258 .filter_map(|d| theta(kem_x[d][8] as i64, kem_a[d][8]).map(|t| (d, t)))259 .collect();260 kem.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap());261 let text: Vec<String> = kem.iter().map(|(d, t)| format!("X{d} {t:.4}")).collect();262 println!(263 "distribution q=10 kempner thetamax at l=8 ascending: {}",264 text.join(" | ")265 );266 let mut tmax_all: Vec<f64> = Vec::new();267 let mut cut_all: Vec<f64> = Vec::new();268 let mut drift_all: Vec<f64> = Vec::new();269 for (fam, out) in jobs.iter().zip(outcomes.iter()) {270 let l = fam.lmax;271 if let Some(t) = theta(out.mmax[l] as i64, out.counts[l]) {272 tmax_all.push(t);273 if let Some(c) = theta(out.meter[l], out.counts[l]) {274 cut_all.push(c);275 }276 let vals: Vec<f64> = (l - 4..=l)277 .filter_map(|v| theta(out.mmax[v] as i64, out.counts[v]))278 .collect();279 let lo = vals.iter().cloned().fold(f64::MAX, f64::min);280 let hi = vals.iter().cloned().fold(f64::MIN, f64::max);281 drift_all.push(hi - lo);282 }283 }284 for d in 0..10 {285 if let Some(t) = theta(kem_x[d][8] as i64, kem_a[d][8]) {286 tmax_all.push(t);287 }288 if let Some(c) = theta(kem_m[d][8], kem_a[d][8]) {289 cut_all.push(c);290 }291 }292 let mut ctl_all: Vec<f64> = Vec::new();293 for c in &ctrl {294 let mut p = 1u64;295 for _ in 0..c.levels {296 p *= c.q;297 }298 if let Some(t) = theta(c.mmax[c.levels] as i64, p) {299 ctl_all.push(t);300 }301 }302 let band = |v: &[f64]| {303 let lo = v.iter().cloned().fold(f64::MAX, f64::min);304 let hi = v.iter().cloned().fold(f64::MIN, f64::max);305 (lo, hi)306 };307 let (tl, th) = band(&tmax_all);308 let (cl, ch) = band(&cut_all);309 let (dl, dh) = band(&drift_all);310 let (gl, gh) = band(&ctl_all);311 let dev = tmax_all312 .iter()313 .map(|t| (t - 0.5).abs())314 .fold(0.0f64, f64::max);315 println!(316 "band families={} thetamax {tl:.4}..{th:.4} maxdev {dev:.4} driftmax5 {dl:.4}..{dh:.4} cut {cl:.4}..{ch:.4} controls {gl:.4}..{gh:.4}",317 tmax_all.len()318 );319 println!("mobius-designs all checks pass");320}321322// SCALING LAW323324fn verify_scalings(jobs: &[Family], outcomes: &[Outcome]) {325 for (fam, out) in jobs.iter().zip(outcomes.iter()) {326 let g = digit_gcd(&fam.digits);327 if g < 2 {328 continue;329 }330 let base_digits: Vec<u64> = fam.digits.iter().map(|d| d / g).collect();331 let (bi, base) = jobs332 .iter()333 .enumerate()334 .find(|(_, f)| f.q == fam.q && f.digits == base_digits)335 .unwrap();336 let bout = &outcomes[bi];337 let (_, expected) = bout.twisted.iter().find(|(a, _)| *a == g).unwrap();338 let has01 = base.digits.contains(&0) && base.digits.contains(&1);339 for lev in 1..=fam.lmax {340 assert_eq!(341 out.meter[lev], expected[lev],342 "scaling meter q={} F={} l={lev}",343 fam.q, fam.label344 );345 let base_a = bout.counts[lev] - if has01 { 1 } else { 0 };346 assert_eq!(347 out.counts[lev], base_a,348 "scaling count q={} F={} l={lev}",349 fam.q, fam.label350 );351 }352 println!(353 "identity q={} F={} equals the a={g} twist of F={} at every level 1..{}",354 fam.q, fam.label, base.label, fam.lmax355 );356 }357}358359// TESTS360361#[cfg(test)]362mod tests {363 use super::*;364 use factor::{is_prime, isqrt};365366 #[test]367 fn mu_matches_sieve() {368 let primes = small_primes(1024);369 let mu = mu_sieve(50_000);370 for n in 1..=50_000u64 {371 assert_eq!(mobius(n, &primes), mu[n as usize], "mu({n})");372 }373 }374375 #[test]376 fn strong_pseudoprimes_rejected() {377 for n in [378 2047u64,379 1373653,380 25326001,381 3215031751,382 3474749660383,383 341550071728321,384 ] {385 assert!(!is_prime(n), "{n} is composite");386 }387 for n in [388 2u64,389 61,390 1_000_000_007,391 999_999_999_989,392 2_305_843_009_213_693_951,393 ] {394 assert!(is_prime(n), "{n} is prime");395 }396 }397398 #[test]399 fn isqrt_exact() {400 for n in [0u64, 1, 2, 3, 4, 24, 25, 26, 999999999999999999] {401 let r = isqrt(n);402 assert!(r as u128 * r as u128 <= n as u128);403 assert!((r as u128 + 1) * (r as u128 + 1) > n as u128);404 }405 }406407 #[test]408 fn methods_agree_on_family() {409 let primes = small_primes(1024);410 let mu = mu_sieve(6561);411 let mut a = vec![0i64; 9];412 let mut b = vec![0i64; 9];413 sweep(3, &[1, 2], 8, &mut |v, len| {414 a[len] += mobius(v, &primes) as i64;415 b[len] += mu[v as usize] as i64;416 });417 assert_eq!(a, b);418 }419420 #[test]421 fn scaled_family_vanishes() {422 let primes = small_primes(1024);423 let fam = census::make_family_depth(5, vec![0, 4], 8);424 let out = run_family(&fam, &primes);425 for lev in 1..=fam.lmax {426 assert_eq!(out.meter[lev], 0);427 }428 }429430 #[test]431 fn scaling_identity_small() {432 let primes = small_primes(1024);433 for (q, scaled, base, a) in [434 (3u64, vec![0u64, 2], vec![0u64, 1], 2u64),435 (4, vec![0, 3], vec![0, 1], 3),436 ] {437 let fs = census::make_family_depth(q, scaled, 8);438 let fb = census::make_family_depth(q, base, 8);439 let os = run_family(&fs, &primes);440 let ob = run_family(&fb, &primes);441 let (_, expected) = ob.twisted.iter().find(|(c, _)| *c == a).unwrap();442 for lev in 1..=8 {443 assert_eq!(os.meter[lev], expected[lev]);444 }445 }446 }447448 #[test]449 fn counting_closed_form() {450 let primes = small_primes(1024);451 let f01 = census::make_family_depth(3, vec![0, 1], 8);452 let o01 = run_family(&f01, &primes);453 let f12 = census::make_family_depth(3, vec![1, 2], 8);454 let o12 = run_family(&f12, &primes);455 for lev in 1..=8usize {456 assert_eq!(o01.counts[lev], 1 << lev);457 assert_eq!(o12.counts[lev], (1 << (lev + 1)) - 2);458 }459 }460461 #[test]462 fn ordered_enumeration_ascending() {463 for digits in [vec![0u64, 1], vec![0, 2], vec![1, 2]] {464 let mut prev = 0u64;465 for len in 1..=6 {466 census::sweep_length(3, &digits, len, &mut |v| {467 assert!(v > prev);468 prev = v;469 });470 }471 }472 }473474 #[test]475 fn ordered_walk_matches_prefix() {476 let primes = small_primes(1024);477 let mu = mu_sieve(19683);478 let fam = census::make_family_depth(3, vec![0, 1], 8);479 let out = run_family(&fam, &primes);480 let mut cum = vec![0i64; 9];481 let mut cnt = vec![0u64; 9];482 sweep(3, &[0, 1], 8, &mut |v, len| {483 cum[len] += mu[v as usize] as i64;484 cnt[len] += 1;485 });486 let mut run = 0i64;487 let mut tot = 0u64;488 for lev in 1..=8usize {489 run += cum[lev];490 tot += cnt[lev];491 let boundary = if lev == 1 { mu[3] as i64 } else { 0 };492 assert_eq!(out.meter[lev], run + boundary);493 assert_eq!(out.counts[lev], tot + 1);494 assert!(out.mmax[lev] >= out.meter[lev].unsigned_abs());495 assert!(out.mmax[lev] >= out.mmax[lev - 1]);496 }497 }498499 #[test]500 fn mertens_prefix_known() {501 let mu = mu_sieve(100_000);502 let mut run = 0i64;503 let mut got = Vec::new();504 for n in 1..=100_000usize {505 run += mu[n] as i64;506 if n == 10 || n == 100 || n == 1000 || n == 10_000 || n == 100_000 {507 got.push(run);508 }509 }510 assert_eq!(got, vec![-1, 1, 2, -23, -48]);511 }512}