large.rs
16.9 kB · rust · 645 lines
1use num_bigint::BigInt;23use crate::riesz::{alpha, band, families, moment, theta_band, transform_table, Family};4use crate::{bigpow, ratio_f64};56// GRID78pub struct Grid {9 pub l: usize,10 pub n: usize,11 pub k: usize,12 pub fh: Vec<(f64, f64)>,13 pub sq: Vec<f64>,14}1516pub fn grid(fam: &Family, l: usize) -> Grid {17 let q = fam.q as usize;18 let n = q.pow(l as u32);19 let g = transform_table(fam.q, &fam.digits, l);20 let mut fh = vec![(0.0f64, 0.0f64); n];21 let mut sq = vec![0.0f64; n];22 for a in 0..n {23 let mut b = a;24 let (mut re, mut im) = (1.0f64, 0.0f64);25 for _ in 0..l {26 let (gr, gi) = g[b];27 let nr = re * gr - im * gi;28 let ni = re * gi + im * gr;29 re = nr;30 im = ni;31 b = (b * q) % n;32 }33 fh[a] = (re, im);34 sq[a] = re * re + im * im;35 }36 Grid {37 l,38 n,39 k: fam.digits.len(),40 fh,41 sq,42 }43}4445impl Grid {46 fn kl(&self) -> f64 {47 (self.k as f64).powi(self.l as i32)48 }4950 fn l1(&self) -> f64 {51 self.sq.iter().map(|v| v.sqrt()).sum()52 }5354 fn second(&self) -> f64 {55 (1..self.n)56 .map(|a| self.sq[a])57 .fold(0.0f64, f64::max)58 .sqrt()59 / self.kl()60 }6162 pub fn large_set(&self, eta: f64) -> Vec<usize> {63 let t2 = self.kl().powi(2) * (self.n as f64).powf(-2.0 * eta) * (1.0 - 1e-9);64 (0..self.n).filter(|&a| self.sq[a] >= t2).collect()65 }66}6768// ROUNDING6970fn read(v: f64, digits: usize) -> String {71 band(v - 1e-12, v + 1e-12, digits)72}7374fn down(v: f64, digits: usize) -> String {75 band(v - 1e-12, v - 1e-12, digits)76}7778fn up(v: f64, digits: usize) -> String {79 band(v + 1e-12, v + 1e-12, digits)80}8182fn sci(v: f64, upward: bool) -> String {83 if v <= 0.0 {84 return "0".to_string();85 }86 let e = v.log10().floor() as i32;87 let m = v / 10f64.powi(e);88 let scaled = m * 1e5 * if upward { 1.0 + 1e-12 } else { 1.0 - 1e-12 };89 let r = if upward {90 scaled.ceil()91 } else {92 scaled.floor()93 };94 if r >= 1e6 {95 format!("1.00000e{}", e + 1)96 } else {97 format!("{:.5}e{}", r / 1e5, e)98 }99}100101// STRUCTURE102103pub struct Shape {104 pub min_gap: usize,105 pub runs: usize,106 pub max_run: usize,107}108109pub fn shape(set: &[usize], n: usize) -> Shape {110 let m = set.len();111 if m < 2 {112 return Shape {113 min_gap: n,114 runs: m,115 max_run: m,116 };117 }118 let gaps: Vec<usize> = (0..m)119 .map(|i| {120 if i + 1 < m {121 set[i + 1] - set[i]122 } else {123 n - set[i] + set[0]124 }125 })126 .collect();127 let min_gap = *gaps.iter().min().unwrap();128 let breaks = gaps.iter().filter(|&&d| d != 1).count();129 if breaks == 0 {130 return Shape {131 min_gap,132 runs: 1,133 max_run: m,134 };135 }136 let start = (0..m).find(|&i| gaps[(i + m - 1) % m] != 1).unwrap();137 let mut max_run = 0usize;138 let mut run = 0usize;139 for t in 0..m {140 run += 1;141 if gaps[(start + t) % m] != 1 {142 max_run = max_run.max(run);143 run = 0;144 }145 }146 Shape {147 min_gap,148 runs: breaks,149 max_run,150 }151}152153// GRAM154155pub struct Gram {156 pub lo: f64,157 pub hi: f64,158 pub iterations: usize,159}160161pub fn gram_top(g: &Grid, set: &[usize], iterations: usize) -> Gram {162 let m = set.len();163 let entry = |i: usize, j: usize| -> (f64, f64) {164 let d = (set[i] + g.n - set[j]) % g.n;165 g.fh[d]166 };167 let mut gersh = 0.0f64;168 let mut frob = 0.0f64;169 for i in 0..m {170 let mut row = 0.0;171 for j in 0..m {172 let (re, im) = entry(i, j);173 let mag = (re * re + im * im).sqrt();174 row += mag;175 frob += mag * mag;176 }177 gersh = gersh.max(row);178 }179 let q = (g.n as f64).powf(1.0 / g.l as f64).round() as usize;180 let mut level = g.l;181 let mut unit = 1usize;182 while level > 0 && set.iter().all(|&a| a % (unit * q) == 0) {183 unit *= q;184 level -= 1;185 }186 let witness = (g.k as f64).powi((g.l - level) as i32) * (q as f64).powi(level as i32);187 let hi = gersh.min(frob.sqrt()).min(g.n as f64).min(witness);188 let mut v: Vec<(f64, f64)> = (0..m)189 .map(|i| (1.0 + 0.1 * (((i * 7919 + 13) % 97) as f64) / 97.0, 0.0))190 .collect();191 let mut lo = 0.0f64;192 let mut used = 0usize;193 for it in 0..iterations {194 let mut w = vec![(0.0f64, 0.0f64); m];195 for i in 0..m {196 let (mut sr, mut si) = (0.0f64, 0.0f64);197 for j in 0..m {198 let (gr, gi) = entry(i, j);199 let (vr, vi) = v[j];200 sr += gr * vr - gi * vi;201 si += gr * vi + gi * vr;202 }203 w[i] = (sr, si);204 }205 let mut num = 0.0f64;206 let mut den = 0.0f64;207 for i in 0..m {208 let (vr, vi) = v[i];209 let (wr, wi) = w[i];210 num += vr * wr + vi * wi;211 den += vr * vr + vi * vi;212 }213 let ray = num / den;214 used = it + 1;215 let done = (ray - lo).abs() <= 1e-13 * ray.abs();216 lo = lo.max(ray);217 let norm = w.iter().map(|(a, b)| a * a + b * b).sum::<f64>().sqrt();218 v = w.iter().map(|(a, b)| (a / norm, b / norm)).collect();219 if done {220 break;221 }222 }223 Gram {224 lo,225 hi: hi.max(lo),226 iterations: used,227 }228}229230// CELL231232pub struct Cell {233 pub eta: f64,234 pub count: usize,235 pub r9: f64,236 pub parseval: f64,237 pub shape: Shape,238 pub floor: f64,239 pub gram: Option<Gram>,240 pub small: f64,241 pub large: f64,242 pub large_l2: f64,243 pub large9: f64,244 pub c_sup: f64,245 pub c: f64,246}247248pub fn cell(fam: &Family, g: &Grid, s4: f64, eta: f64, cap: usize) -> Cell {249 let x = g.n as f64;250 let lx = x.ln();251 let al = alpha(fam);252 let kl = g.kl();253 let set = g.large_set(eta);254 let count = set.len();255 let r9 = s4 / kl.powi(4) * x.powf(4.0 * eta);256 let parseval = x / kl * x.powf(2.0 * eta);257 assert!(count as f64 <= r9 * (1.0 + 1e-6), "R.9 count fails");258 assert!(259 count as f64 <= parseval * (1.0 + 1e-6),260 "Parseval count fails"261 );262 let floor = set.iter().map(|&a| g.sq[a]).sum::<f64>() / kl;263 let gram = if count <= cap {264 Some(gram_top(g, &set, 300))265 } else {266 None267 };268 if let Some(gr) = &gram {269 assert!(gr.lo >= floor * (1.0 - 1e-9), "Rayleigh below the floor");270 assert!(gr.lo <= x * (1.0 + 1e-9), "Gram above the full grid");271 }272 let small = (al + 0.5 - eta).min((1.0 + al) / 2.0);273 let large = al + (count as f64).ln() / (2.0 * lx);274 let large_l2 = (floor * kl).ln() / (2.0 * lx);275 let large9 = al + r9.ln() / (2.0 * lx);276 let cap = (1.0 + al) / 2.0;277 Cell {278 eta,279 count,280 r9,281 parseval,282 shape: shape(&set, g.n),283 floor,284 gram,285 small,286 large,287 large_l2,288 large9,289 c_sup: al - small.max(large),290 c: al - small.max(large.min(large_l2).min(cap)),291 }292}293294pub fn set_row(fam: &Family, g: &Grid, c: &Cell) -> String {295 let gram = match &c.gram {296 Some(gr) => format!(297 "{}..{} ({})",298 sci(gr.lo, false),299 sci(gr.hi, true),300 gr.iterations301 ),302 None => "skipped".to_string(),303 };304 format!(305 "| {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} |",306 fam.q,307 fam.label,308 g.l,309 read(c.eta, 6),310 c.count,311 sci(c.r9, true),312 sci(c.parseval, true),313 c.shape.min_gap,314 c.shape.runs,315 c.shape.max_run,316 sci(c.floor, false),317 gram,318 (g.k as u128).pow(g.l as u32),319 g.n320 )321}322323pub fn chain_row(fam: &Family, g: &Grid, c: &Cell, c1: f64) -> String {324 format!(325 "| {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} |",326 fam.q,327 fam.label,328 g.l,329 read(c.eta, 6),330 read(alpha(fam), 6),331 read(c1, 6),332 up(c.small, 6),333 up(c.large, 6),334 up(c.large_l2, 6),335 up(c.large9, 6),336 down(c.c_sup, 6),337 down(c.c, 6)338 )339}340341// WITNESS342343fn member(q: u64, digits: &[u64], l: usize, mut v: u64) -> bool {344 for _ in 0..l {345 if !digits.contains(&(v % q)) {346 return false;347 }348 v /= q;349 }350 v == 0351}352353pub struct Witness {354 pub m: u64,355 pub r: u64,356 pub rq: u64,357 pub cntq: u64,358 pub bal: f64,359}360361pub fn witness(fam: &Family, l: usize) -> Witness {362 let q = fam.q;363 let x = q.pow(l as u32);364 let m = ((x as f64).sqrt() / 2.0).floor() as u64;365 let mut r = 0u64;366 let mut rq = 0u64;367 let mut cntq = 0u64;368 for mm in m..2 * m {369 let div = mm % q == 0;370 if div {371 cntq += 1;372 }373 for nn in m..2 * m {374 if member(q, &fam.digits, l, mm * nn) {375 r += 1;376 if div {377 rq += 1;378 }379 }380 }381 }382 let kl = (fam.digits.len() as f64).powi(l as i32);383 let bal = rq as f64 - kl / x as f64 * cntq as f64 * m as f64;384 Witness {385 m,386 r,387 rq,388 cntq,389 bal,390 }391}392393pub fn witness_row(fam: &Family, l: usize, w: &Witness) -> String {394 let x = (fam.q as f64).powi(l as i32);395 let kl = (fam.digits.len() as f64).powi(l as i32);396 let exp = |v: f64| -> String {397 if v > 0.0 {398 read(v.ln() / x.ln(), 6)399 } else {400 "void".to_string()401 }402 };403 format!(404 "| {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} | {} |",405 fam.q,406 fam.label,407 l,408 w.m,409 w.r,410 read(w.r as f64 / kl, 6),411 exp(w.r as f64),412 read(alpha(fam), 6),413 w.rq,414 w.cntq,415 read(w.bal, 1),416 read(w.bal / kl, 6),417 exp(w.bal.abs())418 )419}420421// RUN422423pub fn cells() -> Vec<(u64, &'static str, usize)> {424 vec![425 (3, "01", 8),426 (3, "01", 12),427 (4, "012", 6),428 (4, "012", 9),429 (5, "0123", 6),430 (5, "0123", 8),431 (10, "ex7", 4),432 (10, "ex7", 6),433 (100, "01", 3),434 (100, "0to49", 2),435 (100, "0to49", 3),436 ]437}438439pub fn etas(fam: &Family, l: usize) -> (f64, f64, Vec<f64>) {440 let m = moment(fam, 4, l.max(4));441 let (tlo, thi) = theta_band(fam, &m);442 let al = alpha(fam);443 let eta4 = (1.0 + 3.0 * al - (tlo + thi) / 2.0) / 2.0;444 let s4 = ratio_f64(445 &(bigpow(&BigInt::from(fam.q), l) * &m.emod[l]),446 &BigInt::from(1),447 );448 (s4, eta4, vec![0.0, eta4 / 2.0, eta4, 2.0 * eta4, 0.25, 0.5])449}450451pub fn cell_row(fam: &Family, g: &Grid, s4: f64, eta4: f64, c1: f64) -> String {452 let x = g.n as f64;453 let second = g.second();454 format!(455 "| {} | {} | {} | {} | {} | {} | {} | {} | {} |",456 fam.q,457 fam.label,458 g.l,459 read(alpha(fam), 6),460 read((s4 / g.kl().powi(4)).ln() / x.ln(), 6),461 read(eta4, 6),462 read(c1, 6),463 read(second, 6),464 read(-second.ln() / x.ln(), 6)465 )466}467468pub fn run() {469 let fams = families();470 let cap = 1500usize;471 let mut sets = Vec::new();472 let mut chains = Vec::new();473 let mut wits = Vec::new();474 let mut heads = Vec::new();475 for (q, label, l) in cells() {476 let fam = fams.iter().find(|f| f.q == q && f.label == label).unwrap();477 let g = grid(fam, l);478 let x = g.n as f64;479 let kl = g.kl();480 let (s4, eta4, etas) = etas(fam, l);481 let total2: f64 = g.sq.iter().sum();482 let total4: f64 = g.sq.iter().map(|v| v * v).sum();483 assert!((total2 - x * kl).abs() <= 1e-10 * x * kl, "Parseval fails");484 assert!((total4 - s4).abs() <= 1e-10 * s4, "fourth moment fails");485 let c1 = g.l1().ln() / x.ln();486 heads.push(cell_row(fam, &g, s4, eta4, c1));487 for eta in etas {488 let c = cell(fam, &g, s4, eta, cap);489 sets.push(set_row(fam, &g, &c));490 chains.push(chain_row(fam, &g, &c, c1));491 }492 let w = witness(fam, l);493 wits.push(witness_row(fam, l, &w));494 }495 println!("riesz large values cells");496 println!("| q | F | L | alpha | nu_4(L) | eta_4 | c_1 | second/k^L | eta_1 |");497 for r in heads {498 println!("{r}");499 }500 println!("riesz large values");501 println!("| q | F | L | eta | #A | R.9 | Parseval | min gap | runs | max run | Delta_L floor | Delta_L exact | k^L | q^L |");502 for r in sets {503 println!("{r}");504 }505 println!("riesz large values chain");506 println!(507 "| q | F | L | eta | alpha | c_1 | small | large sup | large l2 | large R.9 | c sup | c |"508 );509 for r in chains {510 println!("{r}");511 }512 println!("riesz large values witness");513 println!("| q | F | L | M | R | R/k^L | log R/log x | alpha | R_q | multiples of q | bal | bal/k^L | log abs bal/log x |");514 for r in wits {515 println!("{r}");516 }517}518519// TESTS520521#[cfg(test)]522mod tests {523 use super::*;524525 fn brute(q: u64, digits: &[u64], l: usize) -> Vec<(f64, f64)> {526 let n = (q as usize).pow(l as u32);527 let mut strings = vec![0u64];528 for _ in 0..l {529 let mut next = Vec::new();530 for &s in &strings {531 for &f in digits {532 next.push(s * q + f);533 }534 }535 strings = next;536 }537 (0..n)538 .map(|a| {539 let y = a as f64 / n as f64;540 let mut re = 0.0;541 let mut im = 0.0;542 for &s in &strings {543 let t = 2.0 * std::f64::consts::PI * (s as f64) * y;544 re += t.cos();545 im += t.sin();546 }547 (re, im)548 })549 .collect()550 }551552 fn fam(q: u64, label: &str) -> Family {553 families()554 .into_iter()555 .find(|f| f.q == q && f.label == label)556 .unwrap()557 }558559 #[test]560 fn grid_matches_brute() {561 for (q, label, l) in [(3, "01", 6), (4, "012", 4), (5, "0123", 3), (10, "ex7", 3)] {562 let f = fam(q, label);563 let g = grid(&f, l);564 let b = brute(q, &f.digits, l);565 for a in 0..g.n {566 assert!((g.fh[a].0 - b[a].0).abs() < 1e-7 && (g.fh[a].1 - b[a].1).abs() < 1e-7);567 }568 }569 }570571 #[test]572 fn large_set_matches_brute() {573 let f = fam(3, "01");574 let g = grid(&f, 6);575 let b = brute(3, &f.digits, 6);576 for eta in [0.0, 0.05, 0.1, 0.2, 0.5] {577 let t = 64.0 * 729f64.powf(-eta);578 let direct: Vec<usize> = (0..729)579 .filter(|&a| (b[a].0 * b[a].0 + b[a].1 * b[a].1).sqrt() >= t * (1.0 - 1e-9))580 .collect();581 assert_eq!(g.large_set(eta), direct);582 }583 }584585 #[test]586 fn gram_pins() {587 let f = fam(3, "01");588 let g = grid(&f, 6);589 let one = gram_top(&g, &[0], 10);590 assert!((one.lo - 64.0).abs() < 1e-9 && (one.hi - 64.0).abs() < 1e-9);591 let sub: Vec<usize> = (0..9).map(|r| r * 81).collect();592 let w = gram_top(&g, &sub, 300);593 assert!((w.lo - 144.0).abs() < 1e-6 && w.hi >= 144.0 - 1e-6 && w.hi <= 144.0 + 1e-6);594 let part: Vec<usize> = vec![0, 81, 162];595 let pw = gram_top(&g, &part, 300);596 assert!(pw.hi <= 144.0 + 1e-6 && pw.lo <= pw.hi);597 assert_eq!(sci(0.38110351, false), "3.81103e-1");598 assert_eq!(sci(0.38110351, true), "3.81104e-1");599 assert_eq!(sci(99999.96, true), "1.00000e5");600 let full: Vec<usize> = (0..729).collect();601 let fg = gram_top(&g, &full, 300);602 assert!((fg.lo - 729.0).abs() < 1e-6 && fg.hi <= 729.0 * (1.0 + 1e-9));603 }604605 #[test]606 fn shape_pins() {607 let s = shape(&[0, 1, 2, 10, 11, 20], 30);608 assert_eq!((s.min_gap, s.runs, s.max_run), (1, 3, 3));609 let w = shape(&[0, 1, 28, 29], 30);610 assert_eq!((w.min_gap, w.runs, w.max_run), (1, 1, 4));611 let e = shape(&[5], 30);612 assert_eq!((e.min_gap, e.runs, e.max_run), (30, 1, 1));613 }614615 #[test]616 fn witness_matches_enumeration() {617 let f = fam(3, "01");618 let l = 8;619 let mut strings = vec![0u64];620 for _ in 0..l {621 let mut next = Vec::new();622 for &s in &strings {623 for &d in &f.digits {624 next.push(s * 3 + d);625 }626 }627 strings = next;628 }629 let w = witness(&f, l);630 let mut r = 0;631 let mut rq = 0;632 for mm in w.m..2 * w.m {633 for nn in w.m..2 * w.m {634 if strings.contains(&(mm * nn)) {635 r += 1;636 if mm % 3 == 0 {637 rq += 1;638 }639 }640 }641 }642 assert_eq!((w.r, w.rq), (r, rq));643 assert_eq!(w.cntq, (w.m..2 * w.m).filter(|m| m % 3 == 0).count() as u64);644 }645}