vaughan.rs
14.0 kB · rust · 425 lines
1use crate::menergy::{bits_of, boxes, column, qpow, Bits};2use crate::signed::{cells, mobius, spf_table, sweep};34// THE PIECES56pub fn iroot(n: u128, k: u32) -> u64 {7 let mut r = (n as f64).powf(1.0 / k as f64) as u64;8 while r > 0 && (r as u128).pow(k) > n {9 r -= 1;10 }11 while ((r + 1) as u128).pow(k) <= n {12 r += 1;13 }14 r15}1617pub fn pieces(mu: &[i8], lim: usize, u: u64) -> (Vec<f64>, Vec<f64>) {18 let mut b = vec![0.0f64; lim + 1];19 let mut c = vec![0.0f64; lim + 1];20 let cap = (u as usize).min(lim);21 for d in 1..=cap {22 if mu[d] == 0 {23 continue;24 }25 let m = mu[d] as f64;26 let ld = (d as f64).ln();27 let mut l = d;28 while l <= lim {29 b[l] += m;30 c[l] += m * ld;31 l += d;32 }33 }34 let blog: Vec<f64> = (0..=lim)35 .map(|l| {36 if l == 0 {37 0.038 } else {39 (l as f64).ln() * b[l] - c[l]40 }41 })42 .collect();43 (b, blog)44}4546pub fn signs(b: &[f64]) -> Vec<f64> {47 b.iter()48 .map(|&v| {49 if v > 0.0 {50 1.051 } else if v < 0.0 {52 -1.053 } else {54 0.055 }56 })57 .collect()58}5960// THE ROW6162pub struct VaughanRow {63 pub r: u128,64 pub form: Vec<f64>,65 pub diag: Vec<f64>,66 pub bmax: Vec<f64>,67 pub nz: Vec<usize>,68}6970pub fn vaughan_row(71 bits: &Bits,72 mm: u64,73 nn: u64,74 g: u128,75 x: u128,76 coefs: &[Vec<f64>],77) -> VaughanRow {78 let nb = coefs.len();79 let delta = g as f64 / x as f64;80 let bsum: Vec<f64> = coefs81 .iter()82 .map(|b| (nn..2 * nn).map(|l| b[l as usize]).sum())83 .collect();84 let mut cprime = vec![0u32; nn as usize];85 let mut form = vec![0.0f64; nb];86 let mut u = vec![0.0f64; nb];87 let mut r: u128 = 0;88 for m in mm..2 * mm {89 u.iter_mut().for_each(|v| *v = 0.0);90 for l in nn..2 * nn {91 if bits.get((m * l) as usize) {92 r += 1;93 cprime[(l - nn) as usize] += 1;94 for (j, b) in coefs.iter().enumerate() {95 u[j] += b[l as usize];96 }97 }98 }99 for j in 0..nb {100 let v = u[j] - delta * bsum[j];101 form[j] += v * v;102 }103 }104 let mut diag = vec![0.0f64; nb];105 let mut bmax = vec![0.0f64; nb];106 let mut nz = vec![0usize; nb];107 for (j, b) in coefs.iter().enumerate() {108 let mut d = 0.0f64;109 for l in nn..2 * nn {110 let v = b[l as usize];111 if v != 0.0 {112 nz[j] += 1;113 }114 if v.abs() > bmax[j] {115 bmax[j] = v.abs();116 }117 let c = cprime[(l - nn) as usize] as f64;118 d += v * v * ((1.0 - delta) * (1.0 - delta) * c + delta * delta * (mm as f64 - c));119 }120 diag[j] = d;121 }122 VaughanRow {123 r,124 form,125 diag,126 bmax,127 nz,128 }129}130131// THE STUDY132133const BAND: (f64, f64) = (0.2343, 2.4727);134135pub fn run() {136 println!("menergy signed vaughan");137 println!("| q | F | L | box | live/seen | M | N | R | U3 | U5 | 2N/U5 | nz3 | nz5 | bmax3 | bmax5 | fd3 | fd5 | fdlog3 | fdlog5 | fdsgn | w5 | bf5 | raw5 |");138 let mut count = 0usize;139 let mut inband = 0usize;140 let mut under = 0usize;141 let mut inrange = 0usize;142 let mut boxrows = 0usize;143 let mut fdlo = f64::INFINITY;144 let mut fdhi = f64::NEG_INFINITY;145 let mut bflo = f64::INFINITY;146 let mut bfhi = f64::NEG_INFINITY;147 let mut rawlo = f64::INFINITY;148 let mut rawhi = f64::NEG_INFINITY;149 let mut nzlo = f64::INFINITY;150 let mut nzhi = f64::NEG_INFINITY;151 let mut fd5 = (f64::INFINITY, f64::NEG_INFINITY);152 let mut fd5l5 = (f64::INFINITY, f64::NEG_INFINITY);153 let mut bfl5 = (f64::INFINITY, f64::NEG_INFINITY);154 let mut pml5 = 0usize;155 let mut piece = [0usize; 5];156 let mut edge = f64::INFINITY;157 let mut msmall = 0usize;158 let mut nzall = (f64::INFINITY, f64::NEG_INFINITY);159 let mut bmaxall = 0.0f64;160 let mut bmaxhi = 0.0f64;161 let mut plusminus = 0usize;162 let mut wlo = f64::INFINITY;163 let mut whi = f64::NEG_INFINITY;164 for c in cells().iter() {165 for &l in c.depths.iter() {166 let x = qpow(c.q, l);167 let mut lim = 16usize;168 for &(_, nn) in boxes(x).iter() {169 lim = lim.max(2 * nn as usize - 1);170 }171 let vals = column(c.q, &c.digits, l);172 let g = vals.len() as u128;173 let bits = bits_of(&vals, x);174 let ones = vec![vec![1i8; lim + 1]];175 let sw = match sweep(&bits, x, g, &ones) {176 Some(s) => s,177 None => continue,178 };179 let spf = spf_table(lim);180 let mu = mobius(lim, &spf);181 let u3 = iroot(x, 3);182 let u5 = iroot(x * x, 5);183 let (b3, b3l) = pieces(&mu, lim, u3);184 let (b5, b5l) = pieces(&mu, lim, u5);185 let sg = signs(&b5);186 let coefs = vec![b3, b5, b3l, b5l, sg];187 for (tag, row) in [Some(("top", &sw.top)), sw.l5.as_ref().map(|r| ("l5", r))]188 .into_iter()189 .flatten()190 {191 let v = vaughan_row(&bits, row.m, row.n, g, x, &coefs);192 assert_eq!(v.r, row.r, "R disagrees with the signed sweep");193 let delta = g as f64 / x as f64;194 let floor = (1.0 - delta) * (row.m as f64 * v.r as f64).sqrt();195 let fd: Vec<f64> = (0..coefs.len())196 .map(|j| {197 if v.diag[j] > 0.0 {198 v.form[j] / v.diag[j]199 } else {200 0.0201 }202 })203 .collect();204 let raw = (row.m as f64 * v.form[1]).sqrt() / floor;205 let bf = if v.bmax[1] > 0.0 {206 raw / v.bmax[1]207 } else {208 0.0209 };210 let w = v.diag[1]211 / (v.bmax[1] * v.bmax[1] * (1.0 - delta) * (1.0 - delta) * v.r as f64);212 assert!(213 (bf - (fd[1] * w).sqrt()).abs() < 1e-9 * (1.0 + bf),214 "bound/floor is not (form/diag times the diagonal share)^(1/2)"215 );216 let share = |j: usize| v.nz[j] as f64 / row.n as f64;217 boxrows += 1;218 if row.n > u5 {219 inrange += 1;220 }221 if row.m > u5 {222 msmall += 1;223 }224 for j in 0..coefs.len() {225 nzall.0 = nzall.0.min(share(j));226 nzall.1 = nzall.1.max(share(j));227 bmaxall = bmaxall.max(v.bmax[j]);228 }229 for j in 0..coefs.len() {230 count += 1;231 fdlo = fdlo.min(fd[j]);232 fdhi = fdhi.max(fd[j]);233 if fd[j] >= BAND.0 && fd[j] <= BAND.1 {234 inband += 1;235 piece[j] += 1;236 }237 edge = edge.min((fd[j] - BAND.0).abs().min((fd[j] - BAND.1).abs()));238 if fd[j] < 0.1 {239 under += 1;240 }241 }242 bflo = bflo.min(bf);243 bfhi = bfhi.max(bf);244 rawlo = rawlo.min(raw);245 rawhi = rawhi.max(raw);246 nzlo = nzlo.min(share(0).min(share(1)));247 nzhi = nzhi.max(share(0).max(share(1)));248 fd5.0 = fd5.0.min(fd[1]);249 fd5.1 = fd5.1.max(fd[1]);250 if tag == "l5" {251 fd5l5.0 = fd5l5.0.min(fd[1]);252 fd5l5.1 = fd5l5.1.max(fd[1]);253 bfl5.0 = bfl5.0.min(bf);254 bfl5.1 = bfl5.1.max(bf);255 if v.bmax[1] == 1.0 {256 pml5 += 1;257 }258 }259 bmaxhi = bmaxhi.max(v.bmax[0].max(v.bmax[1]));260 wlo = wlo.min(w);261 whi = whi.max(w);262 if v.bmax[1] == 1.0 {263 plusminus += 1;264 }265 println!(266 "| {} | {} | {} | {} | {}/{} | {} | {} | {} | {} | {} | {:.2} | {:.4} | {:.4} | {} | {} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} |",267 c.q,268 c.label,269 l,270 tag,271 sw.live,272 sw.seen,273 row.m,274 row.n,275 v.r,276 u3,277 u5,278 2.0 * row.n as f64 / u5 as f64,279 share(0),280 share(1),281 v.bmax[0],282 v.bmax[1],283 fd[0],284 fd[1],285 fd[2],286 fd[3],287 fd[4],288 w,289 bf,290 raw291 );292 }293 }294 }295 println!("menergy signed vaughan summary");296 println!("| box rows | N > U5 | M > U5 | coefficient cells | fd lo | fd hi | inside band | under 0.1 | bf lo | bf hi | raw lo | raw hi | nz lo | nz hi | bmax hi | bmax5 = 1 | at l5 | w lo | w hi | fd5 lo | fd5 hi | fd5 l5 lo | fd5 l5 hi | bf l5 lo | bf l5 hi | in band by piece | band edge | nz all | bmax all |");297 println!(298 "| {} | {} | {} | {} | {:.4} | {:.4} | {} | {} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {} | {} | {} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {:.4} | {} {} {} {} {} | {:.4} | {:.4} {:.4} | {:.4} |",299 boxrows, inrange, msmall, count, fdlo, fdhi, inband, under, bflo, bfhi, rawlo, rawhi, nzlo, nzhi,300 bmaxhi, plusminus, pml5, wlo, whi, fd5.0, fd5.1, fd5l5.0, fd5l5.1, bfl5.0, bfl5.1, piece[0], piece[1], piece[2],301 piece[3], piece[4], edge, nzall.0, nzall.1, bmaxall302 );303}304305// TESTS306307#[cfg(test)]308mod tests {309 use super::*;310311 #[test]312 fn the_vaughan_pieces_match_a_direct_divisor_sum() {313 let lim = 1200usize;314 let spf = spf_table(lim);315 let mu = mobius(lim, &spf);316 for u in [7u64, 40, 300] {317 let (b, bl) = pieces(&mu, lim, u);318 for l in 1..=lim {319 let mut s = 0.0f64;320 let mut sl = 0.0f64;321 for d in 1..=l {322 if l % d == 0 && (d as u64) <= u {323 s += mu[d] as f64;324 sl += mu[d] as f64 * ((l / d) as f64).ln();325 }326 }327 assert!((b[l] - s).abs() < 1e-9, "b at l={l} U={u}");328 assert!(329 (bl[l] - sl).abs() < 1e-9 * (1.0 + sl.abs()),330 "blog at l={l} U={u}"331 );332 }333 }334 }335336 #[test]337 fn the_vaughan_form_matches_the_pair_definition() {338 let (q, digits, l, mm, nn) = (3u64, vec![0u64, 1], 8usize, 8u64, 16u64);339 let x = qpow(q, l);340 let vals = column(q, &digits, l);341 let g = vals.len() as u128;342 let bits = bits_of(&vals, x);343 let lim = 2 * nn as usize;344 let spf = spf_table(lim);345 let mu = mobius(lim, &spf);346 let u = iroot(x, 3);347 let (b, blog) = pieces(&mu, lim, u);348 let sg = signs(&b);349 let coefs = vec![b, blog, sg];350 let row = vaughan_row(&bits, mm, nn, g, x, &coefs);351 let delta = g as f64 / x as f64;352 let psi = |m: u64, l: u64| {353 if bits.get((m * l) as usize) {354 1.0 - delta355 } else {356 -delta357 }358 };359 for (j, c) in coefs.iter().enumerate() {360 let mut full = 0.0f64;361 let mut dg = 0.0f64;362 for l1 in nn..2 * nn {363 for l2 in nn..2 * nn {364 let mut w = 0.0f64;365 for m in mm..2 * mm {366 w += psi(m, l1) * psi(m, l2);367 }368 full += c[l1 as usize] * c[l2 as usize] * w;369 if l1 == l2 {370 dg += c[l1 as usize] * c[l1 as usize] * w;371 }372 }373 }374 assert!(375 (full - row.form[j]).abs() < 1e-9 * (1.0 + full.abs()),376 "form {j}: {full} against {}",377 row.form[j]378 );379 assert!(380 (dg - row.diag[j]).abs() < 1e-9 * (1.0 + dg.abs()),381 "diag {j}: {dg} against {}",382 row.diag[j]383 );384 }385 }386387 #[test]388 fn the_mu_piece_form_is_the_exact_integer() {389 let (q, digits, l, mm, nn) = (5u64, vec![0u64, 1, 2, 3], 5usize, 16u64, 32u64);390 let x = qpow(q, l);391 let vals = column(q, &digits, l);392 let g = vals.len() as u128;393 let bits = bits_of(&vals, x);394 let lim = 2 * nn as usize;395 let spf = spf_table(lim);396 let mu = mobius(lim, &spf);397 let u = iroot(x * x, 5);398 let (b, _) = pieces(&mu, lim, u);399 let row = vaughan_row(&bits, mm, nn, g, x, std::slice::from_ref(&b));400 let bi: Vec<i128> = (0..=lim).map(|i| b[i].round() as i128).collect();401 for i in 1..=lim {402 assert_eq!(bi[i] as f64, b[i], "piece is not integral at {i}");403 }404 let xi = x as i128;405 let gi = g as i128;406 let bs: i128 = (nn..2 * nn).map(|l| bi[l as usize]).sum();407 let mut exact = 0i128;408 for m in mm..2 * mm {409 let mut um = 0i128;410 for l in nn..2 * nn {411 if bits.get((m * l) as usize) {412 um += bi[l as usize];413 }414 }415 let v = xi * um - gi * bs;416 exact += v * v;417 }418 let want = exact as f64 / (xi as f64 * xi as f64);419 assert!(420 (want - row.form[0]).abs() < 1e-9 * (1.0 + want.abs()),421 "{want} against {}",422 row.form[0]423 );424 }425}