design.rs
10.9 kB · rust · 397 lines
1use mrlynum::factor::mobius_sieve;2use mrlynum::lattice::totients;34pub struct Design {5 pub name: &'static str,6 pub base: u64,7 pub mask: u32,8}910pub const BASE3: Design = Design {11 name: "base 3, digits {0,1}",12 base: 3,13 mask: 0b11,14};1516pub const BASE10: Design = Design {17 name: "base 10, digit 9 missing",18 base: 10,19 mask: 0b01_1111_1111,20};2122pub const LADDER3: [usize; 10] = [9, 27, 81, 243, 729, 2187, 6561, 19683, 59049, 177147];2324pub const LADDER10: [usize; 5] = [10, 100, 1000, 10000, 100000];2526impl Design {27 pub fn alpha(&self) -> f64 {28 (self.mask.count_ones() as f64).ln() / (self.base as f64).ln()29 }3031 pub fn holds(&self, mut n: u64) -> bool {32 if n == 0 {33 return true;34 }35 while n > 0 {36 if self.mask & (1u32 << (n % self.base)) == 0 {37 return false;38 }39 n /= self.base;40 }41 true42 }4344 pub fn flags(&self, q: usize) -> Vec<bool> {45 (0..=q).map(|n| self.holds(n as u64)).collect()46 }47}4849// THE INDEPENDENT COUNT5051pub fn card_both(q: usize, keep: &[bool], mu: &[i8]) -> u64 {52 let mut total = 0i64;53 for d in 1..=q {54 if mu[d] == 0 {55 continue;56 }57 let mut run = 0i64;58 let mut acc = 0i64;59 let mut m = d;60 while m <= q {61 if keep[m] {62 run += 1;63 acc += run;64 }65 m += d;66 }67 total += mu[d] as i64 * acc;68 }69 total as u6470}7172pub fn card_den(q: usize, keep: &[bool], phi: &[u64]) -> u64 {73 (1..=q).filter(|&b| keep[b]).map(|b| phi[b]).sum()74}7576// THE METER7778#[derive(Clone, Copy)]79pub struct Lane {80 total: u64,81 seen: u64,82 s1: f64,83 s2: f64,84 last: f64,85 gap: f64,86 gap_lo: f64,87}8889impl Lane {90 fn new(total: u64) -> Self {91 Lane {92 total,93 seen: 0,94 s1: 0.0,95 s2: 0.0,96 last: 0.0,97 gap: 0.0,98 gap_lo: 0.0,99 }100 }101102 fn push(&mut self, value: f64) {103 self.seen += 1;104 let delta = value - self.seen as f64 / self.total as f64;105 self.s1 += delta.abs();106 self.s2 += delta * delta;107 if value - self.last > self.gap {108 self.gap = value - self.last;109 self.gap_lo = self.last;110 }111 self.last = value;112 }113}114115pub fn sweep(q: usize, keep: &[bool]) -> (Lane, Lane, Lane, u64, u64, u64) {116 let mu = mobius_sieve(q);117 let phi = totients(q);118 let all = vec![true; q + 1];119 let n_both = card_both(q, keep, &mu);120 let n_den = card_den(q, keep, &phi);121 let n_full = card_den(q, &all, &phi);122 let mut both = Lane::new(n_both);123 let mut den = Lane::new(n_den);124 let mut full = Lane::new(n_full);125 let order = q as u64;126 let (mut a, mut b, mut c, mut d) = (0u64, 1u64, 1u64, order);127 while c <= order {128 let k = (order + b) / d;129 (a, b, c, d) = (c, d, k * c - a, k * d - b);130 let value = a as f64 / b as f64;131 full.push(value);132 if keep[b as usize] {133 den.push(value);134 if keep[a as usize] {135 both.push(value);136 }137 }138 }139 (both, den, full, n_both, n_den, n_full)140}141142// THE TABLE143144fn verdict(ok: bool) -> &'static str {145 if ok {146 "PASS"147 } else {148 "FAIL"149 }150}151152fn cell(now: f64, was: Option<f64>, span: f64) -> (String, String) {153 match was {154 Some(old) => (155 format!("{:.3}", now / old),156 format!("{:+.3}", (now / old).ln() / span),157 ),158 None => ("-".to_string(), "-".to_string()),159 }160}161162struct Row {163 q: usize,164 nodes: u64,165 sieve: u64,166 s1: f64,167 s2: f64,168 gap: f64,169 gap_lo: f64,170}171172fn header(title: &str, rule: &str) {173 println!("{title}");174 println!(" {rule}");175}176177fn meter_rows(rows: &[Row]) -> (f64, f64, f64) {178 println!(" Q card sieve chk S2 r S2 e_2 S1 r S1 e_1");179 let mut prev: Option<&Row> = None;180 let (mut em, mut e2, mut e1) = (0.0, 0.0, 0.0);181 for row in rows {182 let span = prev183 .map(|p| (row.q as f64 / p.q as f64).ln())184 .unwrap_or(1.0);185 let (r2, x2) = cell(row.s2, prev.map(|p| p.s2), span);186 let (r1, x1) = cell(row.s1, prev.map(|p| p.s1), span);187 if let Some(p) = prev {188 em = (row.nodes as f64 / p.nodes as f64).ln() / span;189 e2 = (row.s2 / p.s2).ln() / span;190 e1 = (row.s1 / p.s1).ln() / span;191 }192 println!(193 " {:>5} {:>11} {:>11} {:>4} {:>10.3e} {:>5} {:>6} {:>10.3e} {:>5} {:>6}",194 row.q,195 row.nodes,196 row.sieve,197 verdict(row.nodes == row.sieve),198 row.s2,199 r2,200 x2,201 row.s1,202 r1,203 x1204 );205 prev = Some(row);206 }207 (em, e2, e1)208}209210fn scale_rows(rows: &[Row], alpha: f64, e: f64) {211 println!(212 " Q card exp card S2 Q^{:.3} S1/Q^{:.3} S1/card S2/card widest gap from",213 e - alpha,214 alpha / 2.0215 );216 let mut prev: Option<&Row> = None;217 for row in rows {218 let xm = match prev {219 Some(p) => {220 let span = (row.q as f64 / p.q as f64).ln();221 format!("{:+.3}", (row.nodes as f64 / p.nodes as f64).ln() / span)222 }223 None => "-".to_string(),224 };225 let q = row.q as f64;226 println!(227 " {:>5} {:>11} {:>8} {:>10.4e} {:>10.4e} {:>11.4e} {:>11.4e} {:>11.5} {:>7.5}",228 row.q,229 row.nodes,230 xm,231 row.s2 * q.powf(e - alpha),232 row.s1 / q.powf(alpha / 2.0),233 row.s1 / row.nodes as f64,234 row.s2 / row.nodes as f64,235 row.gap,236 row.gap_lo237 );238 prev = Some(row);239 }240}241242pub fn ladder(design: &Design, qs: &[usize]) {243 let alpha = design.alpha();244 let top = *qs.last().unwrap();245 let keep = design.flags(top);246 let mut both: Vec<Row> = Vec::new();247 let mut den: Vec<Row> = Vec::new();248 let mut full: Vec<Row> = Vec::new();249 for &q in qs {250 let (lb, ld, lf, nb, nd, nf) = sweep(q, &keep[..=q]);251 for (lane, sieve, sink) in [(lb, nb, &mut both), (ld, nd, &mut den), (lf, nf, &mut full)] {252 sink.push(Row {253 q,254 nodes: lane.seen,255 sieve,256 s1: lane.s1,257 s2: lane.s2,258 gap: lane.gap,259 gap_lo: lane.gap_lo,260 });261 }262 }263 println!();264 println!("{} alpha = {:.6}", design.name.to_uppercase(), alpha);265 println!();266 header(267 " CONVENTION strict",268 "F_Q(S_F) = { a/b reduced : 0 < a <= b <= Q, a in S_F, b in S_F }",269 );270 let read_both = meter_rows(&both);271 println!();272 scale_rows(&both, alpha, 2.0 * alpha);273 println!();274 header(275 " CONVENTION denominator",276 "F_Q(S_F) = { a/b reduced : 0 < a <= b <= Q, b in S_F }",277 );278 let read_den = meter_rows(&den);279 println!();280 scale_rows(&den, alpha, 1.0 + alpha);281 println!();282 header(283 " CONTROL full set",284 "F_Q = { a/b reduced : 0 < a <= b <= Q }",285 );286 let read_full = meter_rows(&full);287 println!();288 scale_rows(&full, 1.0, 2.0);289 println!();290 println!(" THE TRANSPLANTED SHAPE D_Q = #{{b in S_F, b <= Q}} ~ Q^a and card ~ Q^e.");291 println!(" Square-root cancellation in the denominators is a count error of order sqrt(D_Q),");292 println!(293 " which puts e_2 at a - e, and Cauchy-Schwarz on S1 <= sqrt(card S2) caps e_1 at a/2."294 );295 println!(" Franel-Landau are the case a = 1, e = 2. Both are caps on a limsup, never values:");296 println!(297 " the control's own e_2 and e_1 wander across this ladder, so one rung proves nothing."298 );299 println!(" lane a e a - e e_2 a/2 e_1 exp card");300 shape("strict", alpha, 2.0 * alpha, read_both);301 shape("denominator", alpha, 1.0 + alpha, read_den);302 shape("full set", 1.0, 2.0, read_full);303}304305fn shape(lane: &str, alpha: f64, e: f64, read: (f64, f64, f64)) {306 let (em, e2, e1) = read;307 println!(308 " {:<12} {:>5.3} {:>5.3} {:>+6.3} {:>+6.3} {:>6.3} {:>+6.3} {:>+8.3}",309 lane,310 alpha,311 e,312 alpha - e,313 e2,314 alpha / 2.0,315 e1,316 em317 );318}319320pub fn run() {321 println!("THE RESTRICTED FAREY METER");322 println!(" S_F is the design: every base digit of the integer drawn from the digit set.");323 println!(" rho_1 < ... < rho_m ascending, delta_j = rho_j - j/m, S2 = sum delta^2, S1 = sum |delta|.");324 println!(" nodes counts the enumeration, sieve the restricted totient sum built without it.");325 ladder(&BASE3, &LADDER3);326 ladder(&BASE10, &LADDER10);327 println!();328 println!("CUT base 3 stops at 3^11 = 177147 and base 10 at 10^5. Rungs are powers of the");329 println!(330 " base so the design's set is self-similar at every rung; a rung inside a decade"331 );332 println!(" truncates the top digit and the meter jumps. The full-set control costs one");333 println!(334 " Stern-Brocot step per node of F_Q, 9.6e9 steps at the base 3 top: the reach is"335 );336 println!(" the control's, the restricted lanes riding the same walk for free.");337}338339#[cfg(test)]340mod tests {341 use super::*;342 use mrlynum::factor::gcd;343344 const FULL: Design = Design {345 name: "full set, the control",346 base: 10,347 mask: 0b11_1111_1111,348 };349350 fn brute_both(q: usize, design: &Design) -> u64 {351 let mut count = 0u64;352 for b in 1..=q as u64 {353 if !design.holds(b) {354 continue;355 }356 for a in 1..=b {357 if design.holds(a) && gcd(a as usize, b as usize) == 1 {358 count += 1;359 }360 }361 }362 count363 }364365 fn brute_den(q: usize, design: &Design) -> u64 {366 let mut count = 0u64;367 for b in 1..=q as u64 {368 if !design.holds(b) {369 continue;370 }371 for a in 1..=b {372 if gcd(a as usize, b as usize) == 1 {373 count += 1;374 }375 }376 }377 count378 }379380 #[test]381 fn restricted_count_matches_brute_force() {382 for (design, q) in [(&BASE3, 243usize), (&BASE10, 200), (&FULL, 120)] {383 let keep = design.flags(q);384 let mu = mobius_sieve(q);385 assert_eq!(card_both(q, &keep, &mu), brute_both(q, design));386 }387 }388389 #[test]390 fn restricted_totient_sum_matches_brute_force() {391 for (design, q) in [(&BASE3, 243usize), (&BASE10, 200), (&FULL, 120)] {392 let keep = design.flags(q);393 let phi = totients(q);394 assert_eq!(card_den(q, &keep, &phi), brute_den(q, design));395 }396 }397}