riesz.rs
48.5 kB · rust · 1707 lines
1use num_bigint::{BigInt, Sign};23use crate::{bigpow, gcd, ratio_f64};45// DIGIT SUMS67fn sum_dist(digits: &[u64], r: usize) -> Vec<u128> {8 let top = *digits.iter().max().unwrap() as usize;9 let mut out = vec![1u128];10 for _ in 0..r {11 let mut cur = vec![0u128; out.len() + top];12 for (i, &c) in out.iter().enumerate() {13 if c == 0 {14 continue;15 }16 for &f in digits {17 cur[i + f as usize] += c;18 }19 }20 out = cur;21 }22 out23}2425fn digit_energy(digits: &[u64]) -> u128 {26 sum_dist(digits, 2).iter().map(|c| c * c).sum()27}2829// TRANSFER3031pub struct Transfer {32 states: Vec<(usize, usize)>,33 m: Vec<Vec<u128>>,34}3536fn ordered(a: &[u128], q: usize, c: usize, cp: usize, d: usize, dp: usize) -> u128 {37 let mut s = 0u128;38 for rho in 0..q {39 let i = d * q + rho;40 let j = dp * q + rho;41 if i < c || j < cp {42 continue;43 }44 let (i, j) = (i - c, j - cp);45 if i >= a.len() || j >= a.len() {46 continue;47 }48 s += a[i] * a[j];49 }50 s51}5253pub fn transfer(q: u64, digits: &[u64], p: usize) -> Transfer {54 assert!(p >= 2 && p % 2 == 0);55 let r = p / 2;56 let a = sum_dist(digits, r);57 let mut states = Vec::new();58 for c in 0..r {59 for cp in c..r {60 states.push((c, cp));61 }62 }63 let n = states.len();64 let mut m = vec![vec![0u128; n]; n];65 for (i, &(c, cp)) in states.iter().enumerate() {66 for (j, &(d, dp)) in states.iter().enumerate() {67 let mut v = ordered(&a, q as usize, c, cp, d, dp);68 if c != cp {69 v += ordered(&a, q as usize, cp, c, d, dp);70 }71 m[i][j] = v;72 }73 }74 Transfer { states, m }75}7677impl Transfer {78 fn n(&self) -> usize {79 self.states.len()80 }8182 fn step(&self, w: &[BigInt]) -> Vec<BigInt> {83 let n = self.n();84 let mut out = vec![BigInt::from(0); n];85 for i in 0..n {86 if w[i] == BigInt::from(0) {87 continue;88 }89 for j in 0..n {90 if self.m[i][j] != 0 {91 out[j] += &w[i] * BigInt::from(self.m[i][j]);92 }93 }94 }95 out96 }9798 pub fn energies(&self, lmax: usize) -> (Vec<BigInt>, Vec<BigInt>) {99 let n = self.n();100 let mut w = vec![BigInt::from(0); n];101 w[0] = BigInt::from(1);102 let mut emod = Vec::new();103 let mut eint = Vec::new();104 for _ in 0..=lmax {105 let mut sm = BigInt::from(0);106 let mut si = BigInt::from(0);107 for (j, &(d, dp)) in self.states.iter().enumerate() {108 if d == dp {109 si += &w[j];110 sm += &w[j];111 } else {112 sm += &w[j] * 2;113 }114 }115 emod.push(sm);116 eint.push(si);117 w = self.step(&w);118 }119 (emod, eint)120 }121122 fn max_row_sum(&self) -> u128 {123 self.m.iter().map(|r| r.iter().sum::<u128>()).max().unwrap()124 }125126 fn charpoly(&self) -> Vec<BigInt> {127 let n = self.n();128 let a: Vec<Vec<BigInt>> = self129 .m130 .iter()131 .map(|r| r.iter().map(|&v| BigInt::from(v)).collect())132 .collect();133 let ident = |n: usize| -> Vec<Vec<BigInt>> {134 (0..n)135 .map(|i| (0..n).map(|j| BigInt::from((i == j) as u8)).collect())136 .collect()137 };138 let mul = |x: &Vec<Vec<BigInt>>, y: &Vec<Vec<BigInt>>| -> Vec<Vec<BigInt>> {139 let mut z = vec![vec![BigInt::from(0); n]; n];140 for i in 0..n {141 for l in 0..n {142 if x[i][l] == BigInt::from(0) {143 continue;144 }145 for j in 0..n {146 z[i][j] += &x[i][l] * &y[l][j];147 }148 }149 }150 z151 };152 let mut c = vec![BigInt::from(0); n + 1];153 c[n] = BigInt::from(1);154 let mut mk = ident(n);155 for step in 1..=n {156 if step > 1 {157 let am = mul(&a, &mk);158 mk = am;159 for i in 0..n {160 mk[i][i] += &c[n - step + 1];161 }162 }163 let am = mul(&a, &mk);164 let mut tr = BigInt::from(0);165 for i in 0..n {166 tr += &am[i][i];167 }168 let kk = BigInt::from(step as u64);169 assert!(&tr % &kk == BigInt::from(0), "leverrier division not exact");170 c[n - step] = -(&tr / &kk);171 }172 c173 }174175 fn reachable(&self) -> Vec<usize> {176 let n = self.n();177 let mut seen = vec![false; n];178 let mut stack = vec![0usize];179 seen[0] = true;180 while let Some(i) = stack.pop() {181 for j in 0..n {182 if self.m[i][j] != 0 && !seen[j] {183 seen[j] = true;184 stack.push(j);185 }186 }187 }188 (0..n).filter(|&i| seen[i]).collect()189 }190191 fn power(&self, idx: &[usize]) -> Vec<f64> {192 let n = idx.len();193 let mut v = vec![1.0f64; n];194 for _ in 0..4000 {195 let mut nv = vec![0.0f64; n];196 for a in 0..n {197 let mut s = v[a];198 for b in 0..n {199 s += self.m[idx[a]][idx[b]] as f64 * v[b];200 }201 nv[a] = s;202 }203 let mx = nv.iter().cloned().fold(0.0f64, f64::max);204 for x in nv.iter_mut() {205 *x /= mx;206 }207 v = nv;208 }209 v210 }211212 pub fn bracket(&self) -> (Frac, Frac) {213 let reach = self.reachable();214 let v0 = self.power(&reach);215 let idx: Vec<usize> = reach216 .iter()217 .zip(v0.iter())218 .filter(|(_, &x)| x > 1e-9)219 .map(|(&i, _)| i)220 .collect();221 let v = self.power(&idx);222 let n = idx.len();223 let vv: Vec<BigInt> = v224 .iter()225 .map(|&x| {226 let s = (x * 2f64.powi(60)).round() as u128;227 BigInt::from(s.max(1))228 })229 .collect();230 let mut lo: Option<Frac> = None;231 let mut hi: Option<Frac> = None;232 for a in 0..n {233 let mut mv = BigInt::from(0);234 for b in 0..n {235 mv += BigInt::from(self.m[idx[a]][idx[b]]) * &vv[b];236 }237 let f = Frac {238 num: mv,239 den: vv[a].clone(),240 };241 lo = Some(match lo {242 None => f.clone(),243 Some(l) => {244 if f.lt(&l) {245 f.clone()246 } else {247 l248 }249 }250 });251 hi = Some(match hi {252 None => f,253 Some(h) => {254 if h.lt(&f) {255 f256 } else {257 h258 }259 }260 });261 }262 (lo.unwrap(), hi.unwrap())263 }264}265266// FRACTIONS267268#[derive(Clone)]269pub struct Frac {270 pub num: BigInt,271 pub den: BigInt,272}273274impl Frac {275 fn int(v: i64) -> Frac {276 Frac {277 num: BigInt::from(v),278 den: BigInt::from(1),279 }280 }281282 fn lt(&self, o: &Frac) -> bool {283 &self.num * &o.den < &o.num * &self.den284 }285286 fn mul(&self, o: &Frac) -> Frac {287 Frac {288 num: &self.num * &o.num,289 den: &self.den * &o.den,290 }291 }292293 fn add(&self, o: &Frac) -> Frac {294 Frac {295 num: &self.num * &o.den + &o.num * &self.den,296 den: &self.den * &o.den,297 }298 }299300 fn sub(&self, o: &Frac) -> Frac {301 Frac {302 num: &self.num * &o.den - &o.num * &self.den,303 den: &self.den * &o.den,304 }305 }306307 fn abs(&self) -> Frac {308 Frac {309 num: if self.num.sign() == Sign::Minus {310 -&self.num311 } else {312 self.num.clone()313 },314 den: self.den.clone(),315 }316 }317318 fn sign(&self) -> Sign {319 self.num.sign()320 }321322 fn to_f64(&self) -> f64 {323 ratio_f64(&self.num, &self.den)324 }325326 fn floor_digits(&self, digits: usize) -> String {327 let scaled = &self.num * bigpow(&BigInt::from(10), digits) / &self.den;328 place(&scaled, digits)329 }330331 fn ceil_digits(&self, digits: usize) -> String {332 let ten = bigpow(&BigInt::from(10), digits);333 let scaled = (&self.num * &ten + &self.den - BigInt::from(1)) / &self.den;334 place(&scaled, digits)335 }336}337338fn eval_frac(c: &[BigInt], x: &Frac) -> Frac {339 let mut acc = Frac::int(0);340 for coef in c.iter().rev() {341 acc = acc.mul(x).add(&Frac {342 num: coef.clone(),343 den: BigInt::from(1),344 });345 }346 acc347}348349fn place(scaled: &BigInt, digits: usize) -> String {350 let s = scaled.to_string();351 if digits == 0 {352 return s;353 }354 let s = if s.len() <= digits {355 format!("{}{}", "0".repeat(digits + 1 - s.len()), s)356 } else {357 s358 };359 let (a, b) = s.split_at(s.len() - digits);360 format!("{a}.{b}")361}362363fn frac_band(lo: &Frac, hi: &Frac, digits: usize) -> String {364 let l = lo.floor_digits(digits);365 let h = hi.ceil_digits(digits);366 if l == h {367 l368 } else {369 format!("{l}..{h}")370 }371}372373pub(crate) fn band(lo: f64, hi: f64, digits: usize) -> String {374 let scale = 10f64.powi(digits as i32);375 let l = (lo * scale).floor() / scale;376 let h = (hi * scale).ceil() / scale;377 if (l - h).abs() < 0.5 / scale {378 format!("{l:.digits$}")379 } else {380 format!("{l:.digits$}..{h:.digits$}")381 }382}383384// FACTORS385386fn eval_big(c: &[BigInt], x: &BigInt) -> BigInt {387 let mut acc = BigInt::from(0);388 for coef in c.iter().rev() {389 acc = acc * x + coef;390 }391 acc392}393394fn eval_i128(c: &[i128], x: i128) -> Option<i128> {395 let mut acc: i128 = 0;396 for &coef in c.iter().rev() {397 acc = acc.checked_mul(x)?.checked_add(coef)?;398 }399 Some(acc)400}401402fn to_i128(c: &[BigInt]) -> Option<Vec<i128>> {403 let mut out = Vec::new();404 for x in c {405 let s = x.to_string();406 out.push(s.parse::<i128>().ok()?);407 }408 Some(out)409}410411fn divide_root(c: &[BigInt], r: &BigInt) -> Vec<BigInt> {412 let n = c.len() - 1;413 let mut qv = vec![BigInt::from(0); n];414 let mut carry = BigInt::from(0);415 for i in (0..n).rev() {416 carry = &carry * r + &c[i + 1];417 qv[i] = carry.clone();418 }419 assert!(&carry * r + &c[0] == BigInt::from(0));420 qv421}422423fn quad_rem_i128(c: &[i128], b: i128, cc: i128) -> Option<(i128, i128)> {424 let n = c.len() - 1;425 let mut q = vec![0i128; n - 1];426 for i in (0..=n - 2).rev() {427 let q1 = if i + 1 < q.len() { q[i + 1] } else { 0 };428 let q2 = if i + 2 < q.len() { q[i + 2] } else { 0 };429 q[i] = c[i + 2]430 .checked_sub(b.checked_mul(q1)?)?431 .checked_sub(cc.checked_mul(q2)?)?;432 }433 let q1 = if 1 < q.len() { q[1] } else { 0 };434 let r1 = c[1]435 .checked_sub(b.checked_mul(q[0])?)?436 .checked_sub(cc.checked_mul(q1)?)?;437 let r0 = c[0].checked_sub(cc.checked_mul(q[0])?)?;438 Some((r1, r0))439}440441fn divide_quad(c: &[BigInt], b: &BigInt, cc: &BigInt) -> Vec<BigInt> {442 let n = c.len() - 1;443 let mut q = vec![BigInt::from(0); n - 1];444 for i in (0..=n - 2).rev() {445 let q1 = if i + 1 < q.len() {446 q[i + 1].clone()447 } else {448 BigInt::from(0)449 };450 let q2 = if i + 2 < q.len() {451 q[i + 2].clone()452 } else {453 BigInt::from(0)454 };455 q[i] = &c[i + 2] - b * q1 - cc * q2;456 }457 let q1 = if 1 < q.len() {458 q[1].clone()459 } else {460 BigInt::from(0)461 };462 assert!(&c[1] - b * &q[0] - cc * q1 == BigInt::from(0));463 assert!(&c[0] - cc * &q[0] == BigInt::from(0));464 q465}466467fn divisors(n: u128) -> Vec<u128> {468 let mut out = vec![1u128];469 let mut m = n;470 let mut f = 2u128;471 while f * f <= m {472 if m % f == 0 {473 let mut e = 0;474 while m % f == 0 {475 m /= f;476 e += 1;477 }478 let base = out.clone();479 let mut pw = 1u128;480 for _ in 0..e {481 pw *= f;482 out.extend(base.iter().map(|d| d * pw));483 }484 }485 f += 1;486 }487 if m > 1 {488 let base = out.clone();489 out.extend(base.iter().map(|d| d * m));490 }491 out.sort();492 out493}494495pub struct Factors {496 pub linear: Vec<BigInt>,497 pub quads: Vec<Vec<BigInt>>,498 pub rest: Vec<BigInt>,499 pub scanned: bool,500 pub witness: Option<u64>,501}502503pub fn factor_bounded(c: &[BigInt], bound: u128) -> Factors {504 let mut cur: Vec<BigInt> = c.to_vec();505 let mut linear = Vec::new();506 let mut quads = Vec::new();507 let small = to_i128(c);508 if small.is_none() && bound > 100000 {509 return Factors {510 linear,511 quads,512 rest: cur,513 scanned: false,514 witness: None,515 };516 }517 let bound_i = bound as i128;518 let mut r: i128 = -bound_i;519 while r <= bound_i {520 let rb = BigInt::from(r);521 let is_root = match small.as_ref().and_then(|s| eval_i128(s, r)) {522 Some(v) => v == 0,523 None => eval_big(c, &rb) == BigInt::from(0),524 };525 if is_root {526 while cur.len() > 1 && eval_big(&cur, &rb) == BigInt::from(0) {527 cur = divide_root(&cur, &rb);528 linear.push(rb.clone());529 }530 }531 r += 1;532 }533 let mut again = cur.len() >= 4;534 while again {535 again = false;536 let ci = match to_i128(&cur) {537 Some(v) => v,538 None => break,539 };540 let a0 = ci[0].unsigned_abs();541 if a0 == 0 || a0 > 1_000_000_000_000_000 {542 break;543 }544 let divs = divisors(a0);545 let bb = bound_i.checked_mul(bound_i).unwrap_or(i128::MAX);546 'search: for &d in &divs {547 if d as i128 > bb {548 break;549 }550 for cc in [d as i128, -(d as i128)] {551 let mut b = -2 * bound_i;552 while b <= 2 * bound_i {553 let hit = match quad_rem_i128(&ci, b, cc) {554 Some((r1, r0)) => r1 == 0 && r0 == 0,555 None => false,556 };557 if hit {558 let bq = BigInt::from(b);559 let cq = BigInt::from(cc);560 cur = divide_quad(&cur, &bq, &cq);561 quads.push(vec![cq, bq, BigInt::from(1)]);562 again = cur.len() >= 4;563 break 'search;564 }565 b += 1;566 }567 }568 }569 }570 let witness = if cur.len() >= 5 {571 irreducibility_witness(&cur)572 } else {573 None574 };575 Factors {576 linear,577 quads,578 rest: cur,579 scanned: true,580 witness,581 }582}583584impl Factors {585 fn all(&self) -> Vec<(Vec<BigInt>, bool)> {586 let mut out: Vec<(Vec<BigInt>, bool)> = Vec::new();587 for r in &self.linear {588 out.push((vec![-r, BigInt::from(1)], true));589 }590 for q in &self.quads {591 out.push((q.clone(), true));592 }593 if self.rest.len() >= 2 {594 let certified = self.scanned && (self.rest.len() <= 4 || self.witness.is_some());595 out.push((self.rest.clone(), certified));596 }597 out598 }599}600601fn nonzero_on(g: &[BigInt], lo: &Frac, hi: &Frac) -> bool {602 let at_lo = eval_frac(g, lo).abs();603 let mut slope = Frac::int(0);604 let mut pw = Frac::int(1);605 for (i, coef) in g.iter().enumerate().skip(1) {606 let term = Frac {607 num: BigInt::from(i as u64) * coef,608 den: BigInt::from(1),609 }610 .abs()611 .mul(&pw);612 slope = slope.add(&term);613 pw = pw.mul(hi);614 }615 let width = hi.sub(lo);616 width.mul(&slope).lt(&at_lo)617}618619pub fn locate(factors: &[(Vec<BigInt>, bool)], lo: &Frac, hi: &Frac) -> Option<usize> {620 let mut found = None;621 for (i, (g, _)) in factors.iter().enumerate() {622 let sl = eval_frac(g, lo).sign();623 let sh = eval_frac(g, hi).sign();624 let changes = sl == Sign::NoSign || sh == Sign::NoSign || sl != sh;625 if changes {626 if found.is_some() {627 return None;628 }629 found = Some(i);630 } else if !nonzero_on(g, lo, hi) {631 return None;632 }633 }634 found635}636637// IRREDUCIBILITY MOD P638639fn pm_trim(a: &mut Vec<u64>) {640 while a.len() > 1 && *a.last().unwrap() == 0 {641 a.pop();642 }643}644645fn pm_rem(a: &[u64], m: &[u64], p: u64) -> Vec<u64> {646 let mut a = a.to_vec();647 pm_trim(&mut a);648 let dm = m.len() - 1;649 if dm == 0 {650 return vec![0];651 }652 let inv = modinv(m[dm], p);653 while a.len() > dm {654 let lead = a[a.len() - 1];655 if lead != 0 {656 let f = lead * inv % p;657 let shift = a.len() - 1 - dm;658 for i in 0..=dm {659 a[shift + i] = (a[shift + i] + p - f * m[i] % p) % p;660 }661 }662 a.pop();663 if a.is_empty() {664 a.push(0);665 }666 }667 pm_trim(&mut a);668 a669}670671fn modinv(a: u64, p: u64) -> u64 {672 let mut r = 1u64;673 let mut b = a % p;674 let mut e = p - 2;675 while e > 0 {676 if e & 1 == 1 {677 r = r * b % p;678 }679 b = b * b % p;680 e >>= 1;681 }682 r683}684685fn pm_mulmod(a: &[u64], b: &[u64], m: &[u64], p: u64) -> Vec<u64> {686 let mut c = vec![0u64; a.len() + b.len() - 1];687 for (i, &x) in a.iter().enumerate() {688 if x == 0 {689 continue;690 }691 for (j, &y) in b.iter().enumerate() {692 c[i + j] = (c[i + j] + x * y) % p;693 }694 }695 pm_rem(&c, m, p)696}697698fn pm_powmod(base: &[u64], mut e: u64, m: &[u64], p: u64) -> Vec<u64> {699 let mut r = vec![1u64];700 let mut b = pm_rem(base, m, p);701 while e > 0 {702 if e & 1 == 1 {703 r = pm_mulmod(&r, &b, m, p);704 }705 b = pm_mulmod(&b, &b, m, p);706 e >>= 1;707 }708 r709}710711fn pm_gcd(a: &[u64], b: &[u64], p: u64) -> Vec<u64> {712 let mut a = a.to_vec();713 let mut b = b.to_vec();714 pm_trim(&mut a);715 pm_trim(&mut b);716 while !(b.len() == 1 && b[0] == 0) {717 let r = pm_rem(&a, &b, p);718 a = b;719 b = r;720 }721 a722}723724fn pm_sub(a: &[u64], b: &[u64], p: u64) -> Vec<u64> {725 let n = a.len().max(b.len());726 let mut c = vec![0u64; n];727 for i in 0..n {728 let x = if i < a.len() { a[i] } else { 0 };729 let y = if i < b.len() { b[i] } else { 0 };730 c[i] = (x + p - y) % p;731 }732 pm_trim(&mut c);733 c734}735736fn reduce_mod(c: &[BigInt], p: u64) -> Vec<u64> {737 let pb = BigInt::from(p);738 c.iter()739 .map(|x| {740 let m = ((x % &pb) + &pb) % &pb;741 m.to_string().parse::<u64>().unwrap()742 })743 .collect()744}745746fn prime_factors(mut d: u64) -> Vec<u64> {747 let mut out = Vec::new();748 let mut f = 2;749 while f * f <= d {750 if d % f == 0 {751 out.push(f);752 while d % f == 0 {753 d /= f;754 }755 }756 f += 1;757 }758 if d > 1 {759 out.push(d);760 }761 out762}763764pub fn irreducible_mod(c: &[BigInt], p: u64) -> bool {765 let m = reduce_mod(c, p);766 let d = m.len() - 1;767 if d < 1 || m[d] == 0 {768 return false;769 }770 if d == 1 {771 return true;772 }773 let deriv: Vec<u64> = (1..=d).map(|i| m[i] * (i as u64 % p) % p).collect();774 let g = pm_gcd(&m, &deriv, p);775 if g.len() > 1 {776 return false;777 }778 let x = vec![0u64, 1];779 let mut powers = vec![x.clone()];780 for _ in 0..d {781 let last = powers.last().unwrap().clone();782 powers.push(pm_powmod(&last, p, &m, p));783 }784 let top = pm_sub(&powers[d], &x, p);785 if !(top.len() == 1 && top[0] == 0) {786 return false;787 }788 for l in prime_factors(d as u64) {789 let i = d / l as usize;790 let diff = pm_sub(&powers[i], &x, p);791 let g = pm_gcd(&m, &diff, p);792 if g.len() > 1 {793 return false;794 }795 }796 true797}798799const PRIMES: [u64; 30] = [800 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,801 101, 103, 107, 109, 113,802];803804pub fn irreducibility_witness(c: &[BigInt]) -> Option<u64> {805 if c.len() <= 2 {806 return Some(0);807 }808 PRIMES.iter().cloned().find(|&p| irreducible_mod(c, p))809}810811// RENDER812813fn poly_string(c: &[BigInt]) -> String {814 let mut parts = Vec::new();815 for i in (0..c.len()).rev() {816 if c[i] == BigInt::from(0) {817 continue;818 }819 let sign = if c[i].sign() == Sign::Minus { "-" } else { "+" };820 let mag = if c[i].sign() == Sign::Minus {821 -&c[i]822 } else {823 c[i].clone()824 };825 let term = match i {826 0 => format!("{mag}"),827 1 => {828 if mag == BigInt::from(1) {829 "x".to_string()830 } else {831 format!("{mag} x")832 }833 }834 _ => {835 if mag == BigInt::from(1) {836 format!("x^{i}")837 } else {838 format!("{mag} x^{i}")839 }840 }841 };842 if parts.is_empty() {843 parts.push(if sign == "-" {844 format!("-{term}")845 } else {846 term847 });848 } else {849 parts.push(format!("{sign} {term}"));850 }851 }852 parts.join(" ")853}854855// STUDY856857pub struct Family {858 pub q: u64,859 pub digits: Vec<u64>,860 pub label: &'static str,861}862863pub struct Moment {864 pub p: usize,865 pub states: usize,866 pub charpoly: Vec<BigInt>,867 pub factors: Factors,868 pub lo: Frac,869 pub hi: Frac,870 pub located: Option<usize>,871 pub emod: Vec<BigInt>,872 pub eint: Vec<BigInt>,873}874875pub fn moment(fam: &Family, p: usize, lmax: usize) -> Moment {876 let t = transfer(fam.q, &fam.digits, p);877 let (emod, eint) = t.energies(lmax);878 let charpoly = t.charpoly();879 let n = t.n();880 for l in 0..=(lmax - n) {881 let mut sm = BigInt::from(0);882 let mut si = BigInt::from(0);883 for i in 0..=n {884 sm += &charpoly[i] * &emod[l + i];885 si += &charpoly[i] * &eint[l + i];886 }887 assert!(sm == BigInt::from(0), "recurrence fails for E_mod");888 assert!(si == BigInt::from(0), "recurrence fails for E_int");889 }890 let factors = factor_bounded(&charpoly, t.max_row_sum());891 let (lo, hi) = t.bracket();892 let located = if factors.scanned {893 locate(&factors.all(), &lo, &hi)894 } else {895 None896 };897 Moment {898 p,899 states: n,900 charpoly,901 factors,902 lo,903 hi,904 located,905 emod,906 eint,907 }908}909910fn assert_bounds(fam: &Family, m: &Moment) {911 let q = BigInt::from(fam.q);912 let k = BigInt::from(fam.digits.len() as u64);913 let r = m.p / 2;914 for l in 0..m.emod.len() {915 let ql = bigpow(&q, l);916 let kl = bigpow(&k, l);917 let s = &ql * &m.emod[l];918 let hoelder = &ql * bigpow(&kl, r);919 let zero = bigpow(&kl, m.p);920 let upper = &ql * bigpow(&kl, m.p - 1);921 assert!(s >= hoelder, "hoelder floor fails");922 assert!(s >= zero, "zero frequency floor fails");923 assert!(s <= upper, "trivial ceiling fails");924 assert!(m.eint[l] <= m.emod[l], "integer energy exceeds modular");925 }926}927928pub(crate) fn alpha(fam: &Family) -> f64 {929 (fam.digits.len() as f64).ln() / (fam.q as f64).ln()930}931932pub(crate) fn theta_band(fam: &Family, m: &Moment) -> (f64, f64) {933 let lq = (fam.q as f64).ln();934 let lo = 1.0 + m.lo.to_f64().ln() / lq - 1e-12;935 let hi = 1.0 + m.hi.to_f64().ln() / lq + 1e-12;936 (lo, hi)937}938939fn minpoly_string(m: &Moment) -> String {940 if !m.factors.scanned {941 return "divides charpoly, scan skipped".to_string();942 }943 let all = m.factors.all();944 match m.located {945 Some(i) => {946 let (g, certified) = &all[i];947 if *certified {948 match (g.len(), m.factors.witness) {949 (l, Some(p)) if l >= 5 => format!("{} [irreducible mod {p}]", poly_string(g)),950 _ => poly_string(g),951 }952 } else {953 format!("divides {} [not certified]", poly_string(g))954 }955 }956 None => "not located".to_string(),957 }958}959960fn factors_string(m: &Moment) -> String {961 let mut parts = Vec::new();962 for r in &m.factors.linear {963 parts.push(format!("(x - {r})").replace("- -", "+ "));964 }965 for q in &m.factors.quads {966 parts.push(format!("({})", poly_string(q)));967 }968 if m.factors.rest.len() >= 2 {969 parts.push(format!("({})", poly_string(&m.factors.rest)));970 }971 parts.join(" ")972}973974pub fn algebra_row(fam: &Family, m: &Moment) -> String {975 let carries = 2 * *fam.digits.iter().max().unwrap() >= fam.q;976 format!(977 "| {} | {} | {} | {} | {} | {} | {} | {} | {} |",978 fam.q,979 fam.label,980 fam.digits.len(),981 m.p,982 if carries { "yes" } else { "no" },983 m.states,984 poly_string(&m.charpoly),985 factors_string(m),986 minpoly_string(m)987 )988}989990pub fn growth_row(fam: &Family, m: &Moment) -> String {991 let q = fam.q as f64;992 let k = fam.digits.len() as f64;993 let (tlo, thi) = theta_band(fam, m);994 let al = alpha(fam);995 let lam_lo = Frac {996 num: &m.lo.num * BigInt::from(fam.q),997 den: m.lo.den.clone(),998 };999 let lam_hi = Frac {1000 num: &m.hi.num * BigInt::from(fam.q),1001 den: m.hi.den.clone(),1002 };1003 let r = (m.p / 2) as i32;1004 let floor = (k.powi(m.p as i32)).max(q * k.powi(r));1005 let ceil = q * k.powi(m.p as i32 - 1);1006 let root_lo = (lam_lo.to_f64()).powf(1.0 / m.p as f64) / k;1007 let root_hi = (lam_hi.to_f64()).powf(1.0 / m.p as f64) / k;1008 format!(1009 "| {} | {} | {} | {} | {} | {} | {} | {} | {} |",1010 fam.q,1011 fam.label,1012 m.p,1013 frac_band(&lam_lo, &lam_hi, 9),1014 floor,1015 ceil,1016 band(tlo, thi, 9),1017 band(al - 1e-12, al + 1e-12, 6),1018 band(root_lo - 1e-12, root_hi + 1e-12, 6)1019 )1020}10211022pub fn typeii_row(fam: &Family, m: &Moment) -> String {1023 let (tlo, thi) = theta_band(fam, m);1024 let al = alpha(fam);1025 let nu_lo = tlo - 4.0 * al - 1e-12;1026 let nu_hi = thi - 4.0 * al + 1e-12;1027 let t2_lo = tlo / 4.0 + 0.25;1028 let t2_hi = thi / 4.0 + 0.25;1029 let miss_lo = t2_lo - al - 1e-12;1030 let miss_hi = t2_hi - al + 1e-12;1031 let need_lo = 1.0 - al / 2.0 - 1e-12;1032 let need_hi = 1.0 - al / 2.0 + 1e-12;1033 let eta_lo = (1.0 + 3.0 * al - thi) / 2.0 - 1e-12;1034 let eta_hi = (1.0 + 3.0 * al - tlo) / 2.0 + 1e-12;1035 format!(1036 "| {} | {} | {} | {} | {} | {} | {} | {} | {} |",1037 fam.q,1038 fam.label,1039 band(al - 1e-12, al + 1e-12, 6),1040 band(tlo, thi, 6),1041 band(nu_lo, nu_hi, 6),1042 band(t2_lo, t2_hi, 6),1043 band(miss_lo, miss_hi, 6),1044 band(need_lo, need_hi, 6),1045 band(eta_lo, eta_hi, 6)1046 )1047}10481049pub fn values_row(fam: &Family, m: &Moment, l: usize) -> String {1050 let s = bigpow(&BigInt::from(fam.q), l) * &m.emod[l];1051 format!(1052 "| {} | {} | {} | {} | {} | {} |",1053 fam.q, fam.label, l, m.emod[l], m.eint[l], s1054 )1055}10561057pub fn ratio_row(fam: &Family, m: &Moment) -> String {1058 let lmax = m.emod.len() - 1;1059 let q = BigInt::from(fam.q);1060 let k4 = bigpow(&BigInt::from(fam.digits.len() as u64), m.p);1061 let s = |l: usize| bigpow(&q, l) * &m.emod[l];1062 let prim = |l: usize| s(l) - &k4 * s(l - 1);1063 let growth = ratio_f64(&m.emod[lmax], &m.emod[lmax - 1]);1064 let prim_growth = {1065 let a = prim(lmax);1066 let b = prim(lmax - 1);1067 if b == BigInt::from(0) {1068 0.01069 } else {1070 ratio_f64(&a, &b)1071 }1072 };1073 let mi = ratio_f64(&m.emod[lmax], &m.eint[lmax]);1074 let rho = m.lo.to_f64();1075 let zero_share = {1076 let z = bigpow(&BigInt::from(fam.digits.len() as u64), m.p * lmax);1077 ratio_f64(&z, &s(lmax))1078 };1079 format!(1080 "| {} | {} | {} | {} | {:.3e} | {:.6} | {:.6} | {:.3e} |",1081 fam.q,1082 fam.label,1083 m.p,1084 lmax,1085 (growth - rho).abs(),1086 mi,1087 prim_growth / (fam.q as f64 * rho),1088 zero_share1089 )1090}10911092// ARCS10931094pub(crate) fn transform_table(q: u64, digits: &[u64], l: usize) -> Vec<(f64, f64)> {1095 let n = (q as usize).pow(l as u32);1096 let mut g = vec![(0.0f64, 0.0f64); n];1097 for b in 0..n {1098 let y = b as f64 / n as f64;1099 let mut re = 0.0;1100 let mut im = 0.0;1101 for &f in digits {1102 let t = 2.0 * std::f64::consts::PI * f as f64 * y;1103 re += t.cos();1104 im += t.sin();1105 }1106 g[b] = (re, im);1107 }1108 g1109}11101111pub struct Arcs {1112 pub l: usize,1113 pub dmax: u64,1114 pub points: usize,1115 pub major: usize,1116 pub parseval_rel: f64,1117 pub fourth_rel: f64,1118 pub major_l2: f64,1119 pub major_l4: f64,1120 pub zero_l4: f64,1121}11221123pub fn arcs(fam: &Family, l: usize, s4_exact: &BigInt) -> Arcs {1124 let q = fam.q;1125 let n = (q as usize).pow(l as u32);1126 let g = transform_table(q, &fam.digits, l);1127 let mut sq = vec![0.0f64; n];1128 for a in 0..n {1129 let mut b = a;1130 let mut v = 1.0f64;1131 for _ in 0..l {1132 let (re, im) = g[b];1133 v *= re * re + im * im;1134 b = (b * q as usize) % n;1135 }1136 sq[a] = v;1137 }1138 let x = n as f64;1139 let dmax = x.powf(0.4).floor() as u64;1140 let qq = x.powf(0.6);1141 let mut major = vec![false; n];1142 for d in 1..=dmax {1143 for lnum in 0..d {1144 if gcd(lnum, d) != 1 {1145 continue;1146 }1147 let c = lnum as f64 * x / d as f64;1148 let h = x / (d as f64 * qq);1149 let lo = (c - h).ceil() as i64;1150 let hi = (c + h).floor() as i64;1151 for a in lo..=hi {1152 let idx = a.rem_euclid(n as i64) as usize;1153 major[idx] = true;1154 }1155 }1156 }1157 let total2: f64 = sq.iter().sum();1158 let total4: f64 = sq.iter().map(|v| v * v).sum();1159 let m2: f64 = (0..n).filter(|&a| major[a]).map(|a| sq[a]).sum();1160 let m4: f64 = (0..n).filter(|&a| major[a]).map(|a| sq[a] * sq[a]).sum();1161 let k = fam.digits.len() as f64;1162 let parseval = x * k.powi(l as i32);1163 let s4 = ratio_f64(s4_exact, &BigInt::from(1));1164 Arcs {1165 l,1166 dmax,1167 points: n,1168 major: major.iter().filter(|&&b| b).count(),1169 parseval_rel: (total2 - parseval).abs() / parseval,1170 fourth_rel: (total4 - s4).abs() / s4,1171 major_l2: m2 / total2,1172 major_l4: m4 / total4,1173 zero_l4: sq[0] * sq[0] / total4,1174 }1175}11761177pub fn arcs_row(fam: &Family, a: &Arcs) -> String {1178 let share = (fam.digits.len() as f64 / fam.q as f64).powi(a.l as i32);1179 format!(1180 "| {} | {} | {} | {} | {} | {} | {:.1e} | {:.1e} | {:.4} | {:.4} | {:.4} | {:.4} |",1181 fam.q,1182 fam.label,1183 a.l,1184 a.dmax,1185 a.points,1186 a.major,1187 a.parseval_rel,1188 a.fourth_rel,1189 1.0 - a.major_l2,1190 share,1191 a.major_l4,1192 a.zero_l41193 )1194}11951196// FAMILIES11971198fn ex(q: u64, e: u64) -> Vec<u64> {1199 (0..q).filter(|&f| f != e).collect()1200}12011202pub fn families() -> Vec<Family> {1203 vec![1204 Family {1205 q: 3,1206 digits: vec![0, 1],1207 label: "01",1208 },1209 Family {1210 q: 3,1211 digits: vec![0, 2],1212 label: "02",1213 },1214 Family {1215 q: 3,1216 digits: vec![1, 2],1217 label: "12",1218 },1219 Family {1220 q: 4,1221 digits: vec![0, 1, 2],1222 label: "012",1223 },1224 Family {1225 q: 5,1226 digits: vec![0, 1, 2, 3],1227 label: "0123",1228 },1229 Family {1230 q: 5,1231 digits: vec![0, 2, 4],1232 label: "024",1233 },1234 Family {1235 q: 10,1236 digits: ex(10, 0),1237 label: "ex0",1238 },1239 Family {1240 q: 10,1241 digits: ex(10, 1),1242 label: "ex1",1243 },1244 Family {1245 q: 10,1246 digits: ex(10, 3),1247 label: "ex3",1248 },1249 Family {1250 q: 10,1251 digits: ex(10, 5),1252 label: "ex5",1253 },1254 Family {1255 q: 10,1256 digits: ex(10, 7),1257 label: "ex7",1258 },1259 Family {1260 q: 10,1261 digits: ex(10, 9),1262 label: "ex9",1263 },1264 Family {1265 q: 100,1266 digits: ex(100, 37),1267 label: "ex37",1268 },1269 Family {1270 q: 100,1271 digits: (0..50).collect(),1272 label: "0to49",1273 },1274 Family {1275 q: 100,1276 digits: vec![0, 1],1277 label: "01",1278 },1279 ]1280}12811282fn find<'a>(fams: &'a [Family], q: u64, label: &str) -> &'a Family {1283 fams.iter().find(|f| f.q == q && f.label == label).unwrap()1284}12851286pub fn run() {1287 let lmax = 60;1288 let fams = families();1289 let fourth: Vec<Moment> = fams.iter().map(|f| moment(f, 4, lmax)).collect();1290 for (f, m) in fams.iter().zip(fourth.iter()) {1291 assert_bounds(f, m);1292 }1293 for (a, b) in [((3, "02"), (3, "01")), ((3, "12"), (3, "01"))] {1294 let ia = fams1295 .iter()1296 .position(|f| f.q == a.0 && f.label == a.1)1297 .unwrap();1298 let ib = fams1299 .iter()1300 .position(|f| f.q == b.0 && f.label == b.1)1301 .unwrap();1302 assert!(fourth[ia].emod == fourth[ib].emod && fourth[ia].eint == fourth[ib].eint);1303 }1304 {1305 let aux = Family {1306 q: 5,1307 digits: vec![0, 1, 2],1308 label: "012",1309 };1310 let ma = moment(&aux, 4, lmax);1311 let i = fams1312 .iter()1313 .position(|f| f.q == 5 && f.label == "024")1314 .unwrap();1315 assert!(fourth[i].emod == ma.emod && fourth[i].eint == ma.eint);1316 }1317 for (f, m) in fams.iter().zip(fourth.iter()) {1318 if 2 * *f.digits.iter().max().unwrap() < f.q {1319 let e = BigInt::from(digit_energy(&f.digits));1320 for l in 0..=lmax {1321 assert!(m.emod[l] == bigpow(&e, l) && m.eint[l] == bigpow(&e, l));1322 }1323 }1324 }1325 println!("riesz algebra");1326 println!(1327 "| q | F | k | p | carries | states | charpoly | factors | minimal polynomial of rho |"1328 );1329 for (f, m) in fams.iter().zip(fourth.iter()) {1330 println!("{}", algebra_row(f, m));1331 }1332 println!("riesz growth");1333 println!("| q | F | p | Lambda(p) | floor | ceiling | theta_p | alpha | Lambda^(1/p)/k |");1334 for (f, m) in fams.iter().zip(fourth.iter()) {1335 println!("{}", growth_row(f, m));1336 }1337 println!("riesz values");1338 println!("| q | F | L | E_mod | E_int | S_4 |");1339 for (f, m) in fams.iter().zip(fourth.iter()) {1340 for l in [1usize, 2, 3] {1341 println!("{}", values_row(f, m, l));1342 }1343 }1344 println!("riesz ratios");1345 println!("| q | F | p | L | growth-rho | E_mod/E_int | prim growth/Lambda | zero share |");1346 for (f, m) in fams.iter().zip(fourth.iter()) {1347 println!("{}", ratio_row(f, m));1348 }1349 println!("riesz type II");1350 println!("| q | F | alpha | theta_4 | nu_4 | L4 exponent | miss | delta needed | eta_4 |");1351 for (f, m) in fams.iter().zip(fourth.iter()) {1352 println!("{}", typeii_row(f, m));1353 }1354 println!("riesz higher moments");1355 println!(1356 "| q | F | k | p | carries | states | charpoly | factors | minimal polynomial of rho |"1357 );1358 let higher: Vec<(&Family, Moment)> = [1359 (3, "01", 6),1360 (3, "01", 8),1361 (3, "01", 10),1362 (3, "12", 6),1363 (4, "012", 6),1364 (5, "0123", 6),1365 (10, "ex7", 6),1366 (100, "01", 6),1367 ]1368 .iter()1369 .map(|&(q, label, p)| {1370 let f = find(&fams, q, label);1371 let m = moment(f, p, lmax);1372 assert_bounds(f, &m);1373 (f, m)1374 })1375 .collect();1376 for (f, m) in higher.iter() {1377 println!("{}", algebra_row(f, m));1378 }1379 println!("| q | F | p | Lambda(p) | floor | ceiling | theta_p | alpha | Lambda^(1/p)/k |");1380 for (f, m) in higher.iter() {1381 println!("{}", growth_row(f, m));1382 }1383 println!("| q | F | p | L | growth-rho | E_mod/E_int | prim growth/Lambda | zero share |");1384 for (f, m) in higher.iter() {1385 println!("{}", ratio_row(f, m));1386 }1387 println!("riesz arcs");1388 println!("| q | F | L | D | points | major | parseval rel | fourth rel | minor l2 | (k/q)^L | major l4 | zero l4 |");1389 for (q, label, l) in [1390 (3, "01", 8),1391 (3, "01", 12),1392 (4, "012", 6),1393 (4, "012", 9),1394 (5, "0123", 6),1395 (5, "0123", 8),1396 (10, "ex7", 4),1397 (10, "ex7", 6),1398 (100, "01", 3),1399 (100, "0to49", 2),1400 (100, "0to49", 3),1401 ] {1402 let i = fams1403 .iter()1404 .position(|f| f.q == q && f.label == label)1405 .unwrap();1406 let s4 = bigpow(&BigInt::from(q), l) * &fourth[i].emod[l];1407 let a = arcs(&fams[i], l, &s4);1408 println!("{}", arcs_row(&fams[i], &a));1409 }1410}14111412// TESTS14131414#[cfg(test)]1415mod tests {1416 use super::*;14171418 fn roots_string(r: &[BigInt]) -> String {1419 r.iter()1420 .map(|x| x.to_string())1421 .collect::<Vec<_>>()1422 .join(",")1423 }14241425 fn strings(q: u64, digits: &[u64], l: usize) -> Vec<u128> {1426 let mut out = vec![0u128];1427 for _ in 0..l {1428 let mut next = Vec::with_capacity(out.len() * digits.len());1429 for &v in &out {1430 for &f in digits {1431 next.push(v * q as u128 + f as u128);1432 }1433 }1434 out = next;1435 }1436 out1437 }14381439 fn histogram_energy(q: u64, digits: &[u64], l: usize, r: usize) -> (u128, u128) {1440 let n = (q as u128).pow(l as u32);1441 let vals = strings(q, digits, l);1442 let mut hist = vec![1u128];1443 for _ in 0..r {1444 let mut next = vec![0u128; hist.len() + n as usize];1445 for (s, &c) in hist.iter().enumerate() {1446 if c == 0 {1447 continue;1448 }1449 for &v in &vals {1450 next[s + v as usize] += c;1451 }1452 }1453 hist = next;1454 }1455 let eint: u128 = hist.iter().map(|c| c * c).sum();1456 let mut modhist = vec![0u128; n as usize];1457 for (s, &c) in hist.iter().enumerate() {1458 modhist[s % n as usize] += c;1459 }1460 let emod: u128 = modhist.iter().map(|c| c * c).sum();1461 (emod, eint)1462 }14631464 fn direct_fourth(q: u64, digits: &[u64], l: usize, grid: usize) -> f64 {1465 let n = (q as usize).pow(l as u32);1466 let mut total = 0.0;1467 for a in 0..grid {1468 let mut v = 1.0f64;1469 for j in 0..l {1470 let y = a as f64 * (q as f64).powi(j as i32) / grid as f64;1471 let mut re = 0.0;1472 let mut im = 0.0;1473 for &f in digits {1474 let t = 2.0 * std::f64::consts::PI * f as f64 * y;1475 re += t.cos();1476 im += t.sin();1477 }1478 v *= re * re + im * im;1479 }1480 total += v * v;1481 }1482 let _ = n;1483 total1484 }14851486 #[test]1487 fn energy_matches_histogram() {1488 for (q, digits, l) in [1489 (3u64, vec![0u64, 1], 7usize),1490 (3, vec![1, 2], 7),1491 (3, vec![0, 2], 6),1492 (4, vec![0, 1, 2], 5),1493 (5, vec![0, 1, 2, 3], 4),1494 (5, vec![0, 2, 4], 4),1495 (10, ex(10, 7), 3),1496 (10, ex(10, 0), 3),1497 (100, ex(100, 37), 1),1498 (100, (0..50).collect(), 1),1499 (100, vec![0, 1], 2),1500 ] {1501 let t = transfer(q, &digits, 4);1502 let (emod, eint) = t.energies(l);1503 for ll in 0..=l {1504 let (hm, hi) = histogram_energy(q, &digits, ll, 2);1505 assert_eq!(emod[ll], BigInt::from(hm), "E_mod q={q} L={ll}");1506 assert_eq!(eint[ll], BigInt::from(hi), "E_int q={q} L={ll}");1507 }1508 }1509 }15101511 #[test]1512 fn higher_energy_matches_histogram() {1513 for (q, digits, p, l) in [1514 (3u64, vec![0u64, 1], 6usize, 5usize),1515 (3, vec![0, 1], 8, 4),1516 (3, vec![1, 2], 6, 4),1517 (4, vec![0, 1, 2], 6, 3),1518 (5, vec![0, 1, 2, 3], 6, 3),1519 ] {1520 let t = transfer(q, &digits, p);1521 let (emod, eint) = t.energies(l);1522 for ll in 0..=l {1523 let (hm, hi) = histogram_energy(q, &digits, ll, p / 2);1524 assert_eq!(emod[ll], BigInt::from(hm));1525 assert_eq!(eint[ll], BigInt::from(hi));1526 }1527 }1528 }15291530 #[test]1531 fn fourth_moment_matches_direct_sum() {1532 for (q, digits, l) in [1533 (3u64, vec![0u64, 1], 6usize),1534 (3, vec![1, 2], 6),1535 (4, vec![0, 1, 2], 5),1536 (5, vec![0, 1, 2, 3], 4),1537 (10, ex(10, 7), 4),1538 (100, vec![0, 1], 2),1539 ] {1540 let t = transfer(q, &digits, 4);1541 let (emod, eint) = t.energies(l);1542 let n = (q as usize).pow(l as u32);1543 let s4 = ratio_f64(&(bigpow(&BigInt::from(q), l) * &emod[l]), &BigInt::from(1));1544 let grid = direct_fourth(q, &digits, l, n);1545 assert!((grid - s4).abs() <= 1e-9 * s4, "S_4 q={q} L={l}");1546 let e_int = ratio_f64(&eint[l], &BigInt::from(1));1547 let integral = direct_fourth(q, &digits, l, 2 * n) / (2 * n) as f64;1548 assert!(1549 (integral - e_int).abs() <= 1e-9 * e_int,1550 "E_int q={q} L={l}"1551 );1552 }1553 }15541555 #[test]1556 fn charpoly_and_scaling_pins() {1557 let f = Family {1558 q: 3,1559 digits: vec![0, 1],1560 label: "01",1561 };1562 let m = moment(&f, 4, 20);1563 assert_eq!(poly_string(&m.charpoly), "x^3 - 8 x^2 + 13 x - 6");1564 assert_eq!(roots_string(&m.factors.linear), "1,1,6");1565 assert_eq!(minpoly_string(&m), "x - 6");1566 assert_eq!(digit_energy(&[0, 1]), 6);1567 assert_eq!(digit_energy(&(0..50).collect::<Vec<u64>>()), 83350);1568 let g = Family {1569 q: 3,1570 digits: vec![1, 2],1571 label: "12",1572 };1573 let mg = moment(&g, 4, 20);1574 assert_eq!(mg.emod, m.emod);1575 assert_eq!(minpoly_string(&mg), "x - 6");1576 let h = Family {1577 q: 3,1578 digits: vec![0, 2],1579 label: "02",1580 };1581 let mh = moment(&h, 4, 20);1582 assert_eq!(poly_string(&mh.charpoly), "x^3 - 15 x^2 + 74 x - 120");1583 assert_eq!(minpoly_string(&mh), "x - 6");1584 }15851586 #[test]1587 fn irreducibility_tester_pins() {1588 let poly = |c: &[i64]| -> Vec<BigInt> { c.iter().map(|&x| BigInt::from(x)).collect() };1589 assert!(irreducible_mod(&poly(&[1, 0, 1]), 3));1590 assert!(!irreducible_mod(&poly(&[1, 0, 1]), 5));1591 assert!(irreducible_mod(&poly(&[-2, 0, 1]), 3));1592 assert!(!irreducible_mod(&poly(&[-2, 0, 1]), 7));1593 assert!(irreducible_mod(&poly(&[-2, 0, 0, 1]), 7));1594 assert!(!irreducible_mod(&poly(&[-1, 0, 0, 1]), 7));1595 for p in [2u64, 3, 5, 7, 11, 13] {1596 assert!(!irreducible_mod(&poly(&[1, 0, 0, 0, 1]), p));1597 }1598 assert!(irreducible_mod(&poly(&[1, 1, 1, 1, 1]), 3));1599 assert_eq!(irreducibility_witness(&poly(&[1, 1, 1, 1, 1])), Some(2));1600 let f = factor_bounded(&poly(&[-6, 11, -6, 1]), 10);1601 assert!(f.scanned);1602 assert_eq!(roots_string(&f.linear), "1,2,3");1603 assert_eq!(f.rest.len(), 1);1604 let g = factor_bounded(&poly(&[2, 0, 3, 0, 1]), 5);1605 assert!(g.linear.is_empty());1606 assert_eq!(g.quads.len(), 1);1607 assert_eq!(poly_string(&g.rest), "x^2 + 2");1608 assert_eq!(divisors(12), vec![1, 2, 3, 4, 6, 12]);1609 }16101611 #[test]1612 fn rows_pinned() {1613 let fams = families();1614 let row = |q: u64, label: &str, p: usize| {1615 let f = find(&fams, q, label);1616 let m = moment(f, p, 30);1617 (1618 algebra_row(f, &m),1619 values_row(f, &m, 3),1620 typeii_row(f, &m),1621 m,1622 )1623 };1624 let (a, v, t, m) = row(3, "01", 4);1625 assert_eq!(a, "| 3 | 01 | 2 | 4 | no | 3 | x^3 - 8 x^2 + 13 x - 6 | (x - 1) (x - 1) (x - 6) | x - 6 |");1626 assert_eq!(v, "| 3 | 01 | 3 | 216 | 216 | 5832 |");1627 assert_eq!(t, "| 3 | 01 | 0.630929..0.630930 | 2.630929..2.630930 | 0.107210..0.107211 | 0.907732..0.907733 | 0.276802..0.276803 | 0.684535..0.684536 | 0.130929..0.130930 |");1628 assert_eq!(frac_band(&m.lo, &m.hi, 9), "6.000000000");1629 let (a, _, _, _) = row(3, "02", 4);1630 assert_eq!(a, "| 3 | 02 | 2 | 4 | yes | 3 | x^3 - 15 x^2 + 74 x - 120 | (x - 4) (x - 5) (x - 6) | x - 6 |");1631 let (a, _, _, _) = row(5, "024", 4);1632 assert_eq!(a, "| 5 | 024 | 3 | 4 | yes | 3 | x^3 - 42 x^2 + 563 x - 2394 | (x - 9) (x - 14) (x - 19) | x - 19 |");1633 let (a, v, _, m) = row(4, "012", 4);1634 assert_eq!(a, "| 4 | 012 | 3 | 4 | yes | 3 | x^3 - 27 x^2 + 136 x - 176 | (x - 4) (x^2 - 23 x + 44) | x^2 - 23 x + 44 |");1635 assert_eq!(v, "| 4 | 012 | 3 | 9173 | 8203 | 587072 |");1636 let rho = (23.0 + 353f64.sqrt()) / 2.0;1637 assert!(m.lo.to_f64() <= rho + 1e-9 && rho - 1e-9 <= m.hi.to_f64());1638 let (a, v, _, m) = row(5, "0123", 4);1639 assert_eq!(a, "| 5 | 0123 | 4 | 4 | yes | 3 | x^3 - 64 x^2 + 659 x - 1476 | (x - 9) (x^2 - 55 x + 164) | x^2 - 55 x + 164 |");1640 assert_eq!(v, "| 5 | 0123 | 3 | 139752 | 116864 | 17469000 |");1641 let rho = (55.0 + 2369f64.sqrt()) / 2.0;1642 assert!(m.lo.to_f64() <= rho + 1e-9 && rho - 1e-9 <= m.hi.to_f64());1643 let (a, v, _, m) = row(10, "ex7", 4);1644 assert_eq!(a, "| 10 | ex7 | 9 | 4 | yes | 3 | x^3 - 731 x^2 + 49474 x - 434304 | (x - 64) (x^2 - 667 x + 6786) | x^2 - 667 x + 6786 |");1645 assert_eq!(v, "| 10 | ex7 | 3 | 283307409 | 188677161 | 283307409000 |");1646 let rho = (667.0 + 417745f64.sqrt()) / 2.0;1647 assert!(m.lo.to_f64() <= rho + 1e-8 && rho - 1e-8 <= m.hi.to_f64());1648 let (a, v, _, m) = row(100, "ex37", 4);1649 assert_eq!(a, "| 100 | ex37 | 99 | 4 | yes | 3 | x^3 - 970301 x^2 + 9322916414 x - 925656819304 | (x - 9604) (x^2 - 960697 x + 96382426) | x^2 - 960697 x + 96382426 |");1650 assert_eq!(v, "| 100 | ex37 | 3 | 886386992219168729 | 588561774108546671 | 886386992219168729000000 |");1651 let rho = (960697.0 + (960697f64 * 960697.0 - 4.0 * 96382426.0).sqrt()) / 2.0;1652 assert!(m.lo.to_f64() <= rho + 1e-4 && rho - 1e-4 <= m.hi.to_f64());1653 let (a, _, _, _) = row(100, "0to49", 4);1654 assert_eq!(a, "| 100 | 0to49 | 50 | 4 | no | 3 | x^3 - 83350 x^2 | (x - 0) (x - 0) (x - 83350) | x - 83350 |");1655 for (q, label, p, want) in [1656 (3, "01", 6, "x^2 - 26 x + 90"),1657 (3, "01", 8, "x^2 - 99 x + 1134"),1658 (3, "01", 10, "x^3 - 392 x^2 + 17469 x - 96228"),1659 (4, "012", 6, "x^2 - 197 x + 2604"),1660 (5, "0123", 6, "x^3 - 858 x^2 + 31475 x - 62700"),1661 (10, "ex7", 6, "x^3 - 53697 x^2 + 29368547 x - 460410771"),1662 (100, "01", 6, "x - 20"),1663 ] {1664 let (_, _, _, m) = row(q, label, p);1665 assert_eq!(minpoly_string(&m), want, "{q} {label} {p}");1666 }1667 let (_, _, _, m) = row(3, "01", 6);1668 let rho = 13.0 + 79f64.sqrt();1669 assert!(m.lo.to_f64() <= rho + 1e-9 && rho - 1e-9 <= m.hi.to_f64());1670 }16711672 #[test]1673 fn arcs_pinned() {1674 let fams = families();1675 let f = find(&fams, 3, "01");1676 let m = moment(f, 4, 12);1677 let s4 = bigpow(&BigInt::from(3), 12) * &m.emod[12];1678 let a = arcs(f, 12, &s4);1679 assert_eq!((a.dmax, a.points, a.major), (195, 531441, 46525));1680 assert!(a.parseval_rel < 1e-10 && a.fourth_rel < 1e-10);1681 assert!((a.major_l2 - 0.2005).abs() < 5e-4);1682 assert!((a.major_l4 - 0.6462).abs() < 5e-4);1683 assert!((a.zero_l4 - (16f64 / 18.0).powi(12)).abs() < 1e-9);1684 }16851686 #[test]1687 fn bracket_contains_growth() {1688 for f in families() {1689 let m = moment(&f, 4, 40);1690 let g = ratio_f64(&m.emod[40], &m.emod[39]);1691 let rho = m.lo.to_f64();1692 assert!(1693 (g - rho).abs() <= 1e-3 * rho,1694 "growth off at {} {}",1695 f.q,1696 f.label1697 );1698 let width = m.hi.to_f64() - rho;1699 assert!(1700 width <= 1e-9 * rho,1701 "bracket too wide at {} {}",1702 f.q,1703 f.label1704 );1705 }1706 }1707}