carry.rs
21.1 kB · rust · 787 lines
1use num_bigint::BigInt;2use std::collections::BTreeMap;34use crate::{bigpow, pow_checked, ratio_f64};56// MODULAR ARITHMETIC78const MR_LIMIT: u128 = 3_317_044_064_679_887_385_961_981;9const MR_BASES: [u128; 13] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41];1011fn addmod(a: u128, b: u128, m: u128) -> u128 {12 let s = a + b;13 if s >= m {14 s - m15 } else {16 s17 }18}1920fn mulmod(a: u128, b: u128, m: u128) -> u128 {21 if m < 1u128 << 64 {22 return (a % m) * (b % m) % m;23 }24 let mut r = 0u128;25 let mut x = a % m;26 let mut y = b % m;27 while y > 0 {28 if y & 1 == 1 {29 r = addmod(r, x, m);30 }31 x = addmod(x, x, m);32 y >>= 1;33 }34 r35}3637fn powmod(a: u128, mut e: u128, m: u128) -> u128 {38 let mut r = 1u128 % m;39 let mut b = a % m;40 while e > 0 {41 if e & 1 == 1 {42 r = mulmod(r, b, m);43 }44 b = mulmod(b, b, m);45 e >>= 1;46 }47 r48}4950fn gcd_u128(mut a: u128, mut b: u128) -> u128 {51 while b != 0 {52 let t = a % b;53 a = b;54 b = t;55 }56 a57}5859// PRIMALITY6061fn is_prime(n: u128) -> bool {62 if n < 2 {63 return false;64 }65 for &p in MR_BASES.iter() {66 if n % p == 0 {67 return n == p;68 }69 }70 let mut d = n - 1;71 let mut r = 0;72 while d % 2 == 0 {73 d /= 2;74 r += 1;75 }76 'outer: for &a in MR_BASES.iter() {77 let mut x = powmod(a, d, n);78 if x == 1 || x == n - 1 {79 continue;80 }81 for _ in 1..r {82 x = mulmod(x, x, n);83 if x == n - 1 {84 continue 'outer;85 }86 }87 return false;88 }89 true90}9192// FACTORISATION9394fn absdiff(a: u128, b: u128) -> u128 {95 if a > b {96 a - b97 } else {98 b - a99 }100}101102fn brent(n: u128, c: u128) -> u128 {103 let step = |x: u128| addmod(mulmod(x, x, n), c % n, n);104 let mut y = 2u128;105 let mut g = 1u128;106 let mut r = 1u128;107 let mut q = 1u128;108 let mut x = 0u128;109 let mut ys = 0u128;110 while g == 1 {111 x = y;112 for _ in 0..r {113 y = step(y);114 }115 let mut k = 0u128;116 while k < r && g == 1 {117 ys = y;118 let lim = if r - k < 128 { r - k } else { 128 };119 for _ in 0..lim {120 y = step(y);121 q = mulmod(q, absdiff(x, y), n);122 }123 g = gcd_u128(q, n);124 k += lim;125 }126 r *= 2;127 }128 if g == n {129 g = 1;130 while g == 1 {131 ys = step(ys);132 let d = absdiff(x, ys);133 if d == 0 {134 return n;135 }136 g = gcd_u128(d, n);137 }138 }139 g140}141142fn split(n: u128) -> u128 {143 let mut c = 1u128;144 loop {145 let d = brent(n, c);146 if d > 1 && d < n {147 return d;148 }149 c += 1;150 }151}152153fn factor_into(n: u128, out: &mut BTreeMap<u128, u32>, unknown: &mut bool) {154 if n == 1 {155 return;156 }157 if is_prime(n) {158 if n >= MR_LIMIT {159 *unknown = true;160 }161 *out.entry(n).or_insert(0) += 1;162 return;163 }164 let d = split(n);165 factor_into(d, out, unknown);166 factor_into(n / d, out, unknown);167}168169fn factorise(mut n: u128) -> (BTreeMap<u128, u32>, bool) {170 assert!(n < 1u128 << 126, "modulus outside the mulmod envelope");171 let mut out = BTreeMap::new();172 let mut unknown = false;173 let mut p = 2u128;174 while p <= 100_000 && p * p <= n {175 while n % p == 0 {176 *out.entry(p).or_insert(0) += 1;177 n /= p;178 }179 p += 1;180 }181 factor_into(n, &mut out, &mut unknown);182 (out, unknown)183}184185fn quot_map(m: &BTreeMap<u128, u32>, g: u64) -> BTreeMap<u128, u32> {186 let mut out = m.clone();187 let mut r = g as u128;188 let mut p = 2u128;189 while r > 1 {190 while r % p == 0 {191 let c = out.get_mut(&p).expect("g fails to divide q^t - 1");192 *c -= 1;193 if *c == 0 {194 out.remove(&p);195 }196 r /= p;197 }198 p += 1;199 }200 out201}202203fn mu_map(m: &BTreeMap<u128, u32>) -> i8 {204 let mut s = 1i8;205 for e in m.values() {206 if *e >= 2 {207 return 0;208 }209 s = -s;210 }211 s212}213214fn fac_string(m: &BTreeMap<u128, u32>) -> String {215 m.iter()216 .map(|(p, e)| {217 if *e == 1 {218 p.to_string()219 } else {220 format!("{p}^{e}")221 }222 })223 .collect::<Vec<String>>()224 .join(" * ")225}226227// CYCLOTOMIC SPLIT228229fn cyclotomic(q: u64, tmax: usize) -> Vec<BigInt> {230 let mut phi: Vec<BigInt> = vec![BigInt::from(1); tmax + 1];231 for d in 1..=tmax {232 let mut v = bigpow(&BigInt::from(q), d) - BigInt::from(1);233 for e in 1..d {234 if d % e == 0 {235 v /= &phi[e];236 }237 }238 phi[d] = v;239 }240 phi241}242243fn to_u128(x: &BigInt) -> u128 {244 x.to_string()245 .parse::<u128>()246 .expect("cyclotomic value overflows u128")247}248249// CARRY DP250251fn poly_powers(digits: &[u64], nmax: usize) -> Vec<Vec<u128>> {252 let top = *digits.iter().max().unwrap() as usize;253 let mut out: Vec<Vec<u128>> = vec![vec![1u128]];254 for n in 1..=nmax {255 let prev = &out[n - 1];256 let mut cur = vec![0u128; prev.len() + top];257 for (i, &c) in prev.iter().enumerate() {258 if c == 0 {259 continue;260 }261 for &f in digits {262 cur[i + f as usize] += c;263 }264 }265 out.push(cur);266 }267 out268}269270fn count_div(q: u64, polys: &[Vec<u128>], l: usize, t: usize, g: u64) -> u128 {271 let s = l / t;272 let u = l - s * t;273 let dig = (q - 1) / g;274 let cmax = s + 1;275 let mut state = vec![0u128; cmax + 1];276 let mut next = vec![0u128; cmax + 1];277 let mut td = vec![0u64; t];278 let mut total: u128 = 0;279 for j in 0..=(s as u64 + 1) * g {280 let mut carry = 0u64;281 for slot in td.iter_mut() {282 let p = dig * j + carry;283 *slot = p % q;284 carry = p / q;285 }286 let high = carry as usize;287 if high > cmax {288 continue;289 }290 state.iter_mut().for_each(|x| *x = 0);291 state[0] = 1;292 for c in 0..t {293 let sc = if c < u { s + 1 } else { s };294 let poly = &polys[sc];295 next.iter_mut().for_each(|x| *x = 0);296 let want = td[c];297 for cr in 0..=cmax {298 let v = state[cr];299 if v == 0 {300 continue;301 }302 let mut m = ((want + q - (cr as u64) % q) % q) as usize;303 while m < poly.len() {304 let w = poly[m];305 if w != 0 {306 let nc = (((cr + m) as u64) - want) / q;307 next[nc as usize] += v * w;308 }309 m += q as usize;310 }311 }312 std::mem::swap(&mut state, &mut next);313 }314 total += state[high];315 }316 total317}318319// STUDY320321struct Ctx {322 q: u64,323 digits: Vec<u64>,324 label: &'static str,325 gs: Vec<u64>,326 polys: Vec<Vec<u128>>,327 qt: Vec<BigInt>,328 fac: Vec<(BTreeMap<u128, u32>, bool)>,329}330331impl Ctx {332 fn new(q: u64, digits: Vec<u64>, label: &'static str, lmax: usize) -> Ctx {333 let gs: Vec<u64> = (1..=q - 1).filter(|g| (q - 1) % g == 0).collect();334 let polys = poly_powers(&digits, lmax);335 let qt: Vec<BigInt> = (0..=lmax)336 .map(|t| bigpow(&BigInt::from(q), t) - BigInt::from(1))337 .collect();338 let phi = cyclotomic(q, lmax);339 let split: Vec<(BTreeMap<u128, u32>, bool)> = (0..=lmax)340 .map(|d| {341 if d == 0 {342 (BTreeMap::new(), false)343 } else {344 factorise(to_u128(&phi[d]))345 }346 })347 .collect();348 let mut fac: Vec<(BTreeMap<u128, u32>, bool)> = Vec::with_capacity(lmax + 1);349 for t in 0..=lmax {350 let mut m: BTreeMap<u128, u32> = BTreeMap::new();351 let mut unk = false;352 for d in 1..=t {353 if t % d == 0 {354 unk |= split[d].1;355 for (p, e) in &split[d].0 {356 *m.entry(*p).or_insert(0) += e;357 }358 }359 }360 fac.push((m, unk));361 }362 Ctx {363 q,364 digits,365 label,366 gs,367 polys,368 qt,369 fac,370 }371 }372373 fn k(&self) -> u64 {374 self.digits.len() as u64375 }376}377378struct Term {379 t: usize,380 e: BigInt,381 mu: i8,382 n: u128,383 et: BigInt,384 aet: BigInt,385}386387struct Res {388 l: usize,389 ratio: f64,390 absnorm: f64,391 live: usize,392 nterms: usize,393 unknown: usize,394 top: Vec<Term>,395}396397fn build(c: &Ctx, l: usize) -> (Vec<Term>, usize) {398 let kl = BigInt::from(pow_checked(c.k(), l));399 let two = BigInt::from(2);400 let zero = BigInt::from(0);401 let mut out: Vec<Term> = Vec::new();402 let mut unknown = 0usize;403 for t in 1..=l {404 for &g in &c.gs {405 let e = &c.qt[t] / g;406 if e < two {407 continue;408 }409 if c.fac[t].1 {410 unknown += 1;411 continue;412 }413 let mu = mu_map("_map(&c.fac[t].0, g));414 if mu == 0 {415 continue;416 }417 let n = count_div(c.q, &c.polys, l, t, g);418 let et = &e * BigInt::from(n) - &kl;419 let aet = if et < zero { -&et } else { et.clone() };420 out.push(Term {421 t,422 e,423 mu,424 n,425 et,426 aet,427 });428 }429 }430 (out, unknown)431}432433fn study(c: &Ctx, l: usize) -> Res {434 let (mut terms, unknown) = build(c, l);435 let kl = BigInt::from(pow_checked(c.k(), l));436 let zero = BigInt::from(0);437 let mut den = BigInt::from(1);438 for t in terms.iter() {439 den *= &t.e;440 }441 let mut sig = BigInt::from(0);442 let mut abs = BigInt::from(0);443 for t in terms.iter() {444 let w = &den / &t.e;445 sig += BigInt::from(t.mu) * &t.et * &w;446 abs += &t.aet * &w;447 }448 let ratio = if abs == zero {449 0.0450 } else {451 let neg = sig < zero;452 let mag = if neg { -&sig } else { sig.clone() };453 let v = ratio_f64(&mag, &abs);454 if neg {455 -v456 } else {457 v458 }459 };460 let absnorm = ratio_f64(&abs, &(&den * &kl));461 terms.sort_by(|a, b| (&b.aet * &a.e).cmp(&(&a.aet * &b.e)));462 let nterms = terms.len();463 let live = if terms.is_empty() || terms[0].aet == zero {464 0465 } else {466 let (bm, be) = (terms[0].aet.clone(), terms[0].e.clone());467 terms468 .iter()469 .filter(|t| BigInt::from(10) * &t.aet * &be >= &bm * &t.e)470 .count()471 };472 let top: Vec<Term> = terms473 .into_iter()474 .take(4)475 .map(|t| Term {476 t: t.t,477 e: t.e,478 mu: t.mu,479 n: t.n,480 et: t.et,481 aet: t.aet,482 })483 .collect();484 Res {485 l,486 ratio,487 absnorm,488 live,489 nterms,490 unknown,491 top,492 }493}494495// RENDER496497fn depth_row(c: &Ctx, r: &Res) -> String {498 format!(499 "| {} | {} | {} | {:+.3} | {:.1e} | {} | {} | {} |",500 c.q, c.label, r.l, r.ratio, r.absnorm, r.live, r.nterms, r.unknown501 )502}503504fn top_row(c: &Ctx, t: &Term) -> String {505 format!(506 "| {} | {} | {} | {} | {:+} | {} | {} |",507 c.q, c.label, t.t, t.e, t.mu, t.n, t.et508 )509}510511fn mu_row(c: &Ctx, g: u64, lmax: usize) -> String {512 let two = BigInt::from(2);513 let mut plus: Vec<String> = Vec::new();514 let mut minus: Vec<String> = Vec::new();515 let mut zeros = 0usize;516 for t in 1..=lmax {517 if &c.qt[t] / g < two {518 continue;519 }520 match mu_map("_map(&c.fac[t].0, g)) {521 1 => plus.push(t.to_string()),522 -1 => minus.push(t.to_string()),523 _ => zeros += 1,524 }525 }526 format!(527 "| {} | Q_t/{} | {} | {} | {} |",528 c.q,529 g,530 plus.join(","),531 minus.join(","),532 zeros533 )534}535536fn fac_row(name: &str, m: &BTreeMap<u128, u32>) -> String {537 format!("| {} | {} |", name, fac_string(m))538}539540fn seq_row(c: &Ctx, res: &[Res], lo: usize, hi: usize) -> String {541 let vals: Vec<String> = res542 .iter()543 .filter(|r| r.l >= lo && r.l <= hi)544 .map(|r| format!("{:+.2}", r.ratio))545 .collect();546 format!(547 "| {} | {} | {}..{} | {} |",548 c.q,549 c.label,550 lo,551 hi,552 vals.join(" ")553 )554}555556fn ex7_digits() -> Vec<u64> {557 (0..10).filter(|&f| f != 7).collect()558}559560fn contexts(lmax: usize) -> Vec<Ctx> {561 vec![562 Ctx::new(3, vec![0, 1], "01", lmax),563 Ctx::new(10, ex7_digits(), "ex7", lmax),564 ]565}566567fn fac_rows(cs: &[Ctx]) -> Vec<String> {568 let mut out = Vec::new();569 for &t in &[7usize, 37, 39] {570 out.push(fac_row(&format!("3^{t} - 1"), &cs[0].fac[t].0));571 }572 for &t in &[19usize, 23, 31, 37] {573 out.push(fac_row(&format!("R_{t}"), "_map(&cs[1].fac[t].0, 9)));574 }575 out576}577578pub fn run() {579 let lmax = 40;580 let cs = contexts(lmax);581 let res: Vec<Vec<Res>> = cs582 .iter()583 .map(|c| (3..=lmax).map(|l| study(c, l)).collect())584 .collect();585 println!("carry depth");586 println!("| q | F | L | ratio | abs/k^L | live | terms | unknown |");587 for (i, c) in cs.iter().enumerate() {588 for &l in &[10usize, 20, 30, 40] {589 println!("{}", depth_row(c, &res[i][l - 3]));590 }591 }592 println!("carry top L=40");593 println!("| q | F | t | e | mu | N | eT |");594 for (i, c) in cs.iter().enumerate() {595 for t in res[i][lmax - 3].top.iter() {596 println!("{}", top_row(c, t));597 }598 }599 println!("carry mobius");600 println!("| q | e | mu=+1 | mu=-1 | zeros |");601 for c in cs.iter() {602 for &g in &c.gs {603 println!("{}", mu_row(c, g, lmax));604 }605 }606 println!("carry factors");607 println!("| n | factorisation |");608 for row in fac_rows(&cs) {609 println!("{row}");610 }611 println!("carry ratio sequence");612 println!("| q | F | L | Sigma_L/Abs_L |");613 for (i, c) in cs.iter().enumerate() {614 for &(lo, hi) in &[(3usize, 12usize), (13, 22), (23, 32), (33, 40)] {615 println!("{}", seq_row(c, &res[i], lo, hi));616 }617 }618}619620// TESTS621622#[cfg(test)]623mod tests {624 use super::*;625 use crate::n_div;626627 fn brute(q: u64, digits: &[u64], l: usize, e: u128) -> u128 {628 fn rec(v: u128, len: usize, q: u64, digits: &[u64], l: usize, e: u128, hits: &mut u128) {629 if len == l {630 if v % e == 0 {631 *hits += 1;632 }633 return;634 }635 for &f in digits {636 rec(v * q as u128 + f as u128, len + 1, q, digits, l, e, hits);637 }638 }639 let mut hits = 0;640 rec(0, 0, q, digits, l, e, &mut hits);641 hits642 }643644 #[test]645 fn carry_matches_brute() {646 for (q, digits, l) in [647 (3u64, vec![0u64, 1], 8usize),648 (10, ex7_digits(), 5),649 (5, vec![0, 2, 4], 6),650 ] {651 let polys = poly_powers(&digits, l);652 let gs: Vec<u64> = (1..=q - 1).filter(|g| (q - 1) % g == 0).collect();653 for t in 1..=l {654 let qt = (q as u128).pow(t as u32) - 1;655 for &g in &gs {656 let e = qt / g as u128;657 if e < 2 {658 continue;659 }660 assert_eq!(661 count_div(q, &polys, l, t, g),662 brute(q, &digits, l, e),663 "q={q} l={l} t={t} g={g}"664 );665 }666 }667 }668 }669670 #[test]671 fn carry_matches_residue_dp() {672 let mut checks = 0;673 for (q, digits) in [(3u64, vec![0u64, 1]), (10, ex7_digits())] {674 let polys = poly_powers(&digits, 12);675 let gs: Vec<u64> = (1..=q - 1).filter(|g| (q - 1) % g == 0).collect();676 for &l in &[8usize, 12] {677 for t in 1..=l {678 let qt = (q as u128).pow(t as u32) - 1;679 for &g in &gs {680 let e = qt / g as u128;681 if e < 2 || e > 30000 {682 continue;683 }684 assert_eq!(685 count_div(q, &polys, l, t, g),686 n_div(q, &digits, l, e as u64),687 "q={q} l={l} t={t} g={g}"688 );689 checks += 1;690 }691 }692 }693 }694 assert!(695 checks >= 40,696 "residue cross-check too thin at {checks} cells"697 );698 }699700 #[test]701 fn mobius_is_certified() {702 for c in contexts(40).iter() {703 for t in 1..=40 {704 assert!(!c.fac[t].1, "uncertified cofactor at q={} t={t}", c.q);705 }706 }707 }708709 #[test]710 fn mobius_rows_pinned() {711 let cs = contexts(40);712 let mut got: Vec<String> = Vec::new();713 for c in cs.iter() {714 for &g in &c.gs {715 got.push(mu_row(c, g, 40));716 }717 }718 got.extend(fac_rows(&cs));719 let want = [720 "| 3 | Q_t/1 | 3,7,13,21,27,29,31 | 1,9,11,17,19,23,33,37 | 25 |",721 "| 3 | Q_t/2 | 9,11,17,19,23,33,37 | 3,7,13,21,27,29,31 | 25 |",722 "| 10 | Q_t/1 | | | 40 |",723 "| 10 | Q_t/3 | 2,13,19,20,23,25,29,31,32,35,37,38,40 | 1,4,5,7,8,10,11,14,16,17,26,28,34 | 14 |",724 "| 10 | Q_t/9 | 3,4,5,7,8,10,11,14,15,16,17,24,26,28,33,34,39 | 2,6,12,13,19,20,21,23,25,29,30,31,32,35,37,38,40 | 5 |",725 "| 3^7 - 1 | 2 * 1093 |",726 "| 3^37 - 1 | 2 * 13097927 * 17189128703 |",727 "| 3^39 - 1 | 2 * 13^2 * 313 * 6553 * 7333 * 797161 |",728 "| R_19 | 1111111111111111111 |",729 "| R_23 | 11111111111111111111111 |",730 "| R_31 | 2791 * 6943319 * 57336415063790604359 |",731 "| R_37 | 2028119 * 247629013 * 2212394296770203368013 |",732 ];733 assert_eq!(got, want);734 }735736 #[test]737 fn carry_rows_pinned() {738 let cs = contexts(40);739 let res: Vec<Vec<Res>> = cs740 .iter()741 .map(|c| (3..=40).map(|l| study(c, l)).collect())742 .collect();743 let mut got: Vec<String> = Vec::new();744 for (i, c) in cs.iter().enumerate() {745 for &l in &[10usize, 20, 30, 40] {746 got.push(depth_row(c, &res[i][l - 3]));747 }748 }749 for (i, c) in cs.iter().enumerate() {750 for t in res[i][37].top.iter() {751 got.push(top_row(c, t));752 }753 }754 for (i, c) in cs.iter().enumerate() {755 for &(lo, hi) in &[(3usize, 12usize), (13, 22), (23, 32), (33, 40)] {756 got.push(seq_row(c, &res[i], lo, hi));757 }758 }759 let want = [760 "| 3 | 01 | 10 | -0.211 | 4.7e-2 | 4 | 7 | 0 |",761 "| 3 | 01 | 20 | -0.123 | 4.1e-3 | 6 | 15 | 0 |",762 "| 3 | 01 | 30 | +0.069 | 1.0e-3 | 6 | 23 | 0 |",763 "| 3 | 01 | 40 | -0.498 | 2.1e-4 | 6 | 29 | 0 |",764 "| 10 | ex7 | 10 | +0.812 | 1.1e-5 | 5 | 15 | 0 |",765 "| 10 | ex7 | 20 | -0.495 | 3.5e-8 | 2 | 31 | 0 |",766 "| 10 | ex7 | 30 | -0.127 | 4.7e-10 | 2 | 44 | 0 |",767 "| 10 | ex7 | 40 | -0.192 | 3.9e-12 | 4 | 60 | 0 |",768 "| 3 | 01 | 7 | 1093 | -1 | 1059181242 | 58173469730 |",769 "| 3 | 01 | 9 | 19682 | -1 | 107075926 | 1007956747756 |",770 "| 3 | 01 | 7 | 2186 | +1 | 460789966 | -92224762100 |",771 "| 3 | 01 | 9 | 9841 | +1 | 148454776 | 361431822840 |",772 "| 10 | ex7 | 5 | 33333 | -1 | 4434309484921670553282063418321580 | 8646548121236467809716529928539 |",773 "| 10 | ex7 | 8 | 33333333 | -1 | 4434111044140794648321659572234 | -5129421023116418959440572469821679 |",774 "| 10 | ex7 | 10 | 3333333333 | -1 | 44419676099006184777353657556 | 256757567534800575426920054498816547 |",775 "| 10 | ex7 | 7 | 3333333 | -1 | 44342590913069608659343530190385 | -207818310865474807662686276694396 |",776 "| 3 | 01 | 3..12 | -0.33 -0.64 -1.00 -1.00 -0.81 -0.83 -0.33 -0.21 -0.09 -0.09 |",777 "| 3 | 01 | 13..22 | -0.13 -0.24 -0.36 -0.41 -0.40 -0.33 -0.26 -0.12 +0.02 +0.14 |",778 "| 3 | 01 | 23..32 | +0.19 +0.18 +0.19 +0.17 +0.12 +0.10 +0.09 +0.07 +0.07 +0.06 |",779 "| 3 | 01 | 33..40 | +0.02 -0.05 -0.17 -0.23 -0.31 -0.38 -0.44 -0.50 |",780 "| 10 | ex7 | 3..12 | +0.76 +0.12 +0.02 +0.15 -0.31 +0.57 +0.53 +0.81 +0.92 +0.13 |",781 "| 10 | ex7 | 13..22 | +0.57 +0.82 +0.20 -0.82 -0.79 -0.73 -0.24 -0.49 -0.86 -0.91 |",782 "| 10 | ex7 | 23..32 | -0.95 -0.31 +0.68 +0.63 +0.58 -0.96 -0.99 -0.13 -0.10 -0.27 |",783 "| 10 | ex7 | 33..40 | -0.31 -0.45 -0.60 -0.64 -0.65 -0.59 -0.23 -0.19 |",784 ];785 assert_eq!(got, want);786 }787}