sieve.rs
15.9 kB · rust · 424 lines
1use std::f64::consts::FRAC_PI_4;23/// The limit of the plane Wallis sieve's surviving area, pi over four.4pub const PLANE_LIMIT: f64 = FRAC_PI_4;56/// Returns the classical Wallis schedule, the odd sides three, five, seven and on, to the count of levels.7///8/// ```9/// assert_eq!(mrlynum::sieve::odd_word(4), vec![3, 5, 7, 9]);10/// ```11pub fn odd_word(levels: usize) -> Vec<u64> {12 (1..=levels as u64).map(|k| 2 * k + 1).collect()13}1415/// Returns the constant schedule, one odd side repeated to the count of levels, whose limit set is the fixed-ratio carpet.16///17/// ```18/// assert_eq!(mrlynum::sieve::flat_word(3, 3), vec![3, 3, 3]);19/// ```20pub fn flat_word(side: u64, levels: usize) -> Vec<u64> {21 vec![side; levels]22}2324fn letter(side: u64) -> u64 {25 assert!(26 side >= 3 && side % 2 == 1,27 "a letter is an odd side from three"28 );29 side30}3132/// Returns the side of the word, the product of its letters' sides, and panics when that overruns a u128.33///34/// ```35/// assert_eq!(mrlynum::sieve::side(&mrlynum::sieve::odd_word(4)), 945);36/// ```37pub fn side(word: &[u64]) -> u128 {38 word.iter().fold(1u128, |run, &s| {39 run.checked_mul(u128::from(letter(s)))40 .expect("the word's side overruns a u128")41 })42}4344/// Returns the cells the word leaves, the product of its letters' fills, one punctured tile a letter, and panics when that overruns a u128.45///46/// ```47/// assert_eq!(mrlynum::sieve::cells(&mrlynum::sieve::odd_word(3), 2), 9216);48/// ```49pub fn cells(word: &[u64], dimension: u32) -> u128 {50 word.iter().fold(1u128, |run, &s| {51 let fill = u128::from(letter(s)).pow(dimension) - 1;52 run.checked_mul(fill)53 .expect("the word's cells overrun a u128")54 })55}5657/// Returns the punctures the word makes, one per surviving cell at every level, and panics when that overruns a u128.58///59/// ```60/// assert_eq!(mrlynum::sieve::holes(&mrlynum::sieve::odd_word(3), 2), 1 + 8 + 192);61/// ```62pub fn holes(word: &[u64], dimension: u32) -> u128 {63 let mut total = 0u128;64 for place in 0..word.len() {65 total += cells(&word[..place], dimension);66 }67 total68}6970/// Returns the share of the whole the word leaves, the product of one minus the inverse of each letter's site count, exact as a product of the letters' fills.71///72/// ```73/// let flat = mrlynum::sieve::ratio(&mrlynum::sieve::flat_word(3, 2), 2);74/// assert!((flat - 64.0 / 81.0).abs() < 1e-15);75/// ```76pub fn ratio(word: &[u64], dimension: u32) -> f64 {77 word.iter().fold(1.0f64, |run, &s| {78 run * (1.0 - 1.0 / (letter(s) as f64).powi(dimension as i32))79 })80}8182/// Returns the box exponent the word reads at its own scale, the logarithm of its cells over the logarithm of its side, which walks up to the dimension on a schedule of distinct growing letters and stands still on any schedule that reuses its letters.83///84/// Changing the letter is not enough: the alternating word 3, 5, 3, 5 changes at every step and freezes at log 192 / log 15, its ratio falling to nothing like a fixed-ratio schedule. The hypothesis that buys a positive area is strictly increasing odd letters, under which sum s_k^(-d) converges.85///86/// ```87/// let carpet = mrlynum::sieve::exponent(&mrlynum::sieve::flat_word(3, 5), 2);88/// assert!((carpet - 8f64.ln() / 3f64.ln()).abs() < 1e-12);89/// ```90pub fn exponent(word: &[u64], dimension: u32) -> f64 {91 let mut up = 0.0f64;92 let mut down = 0.0f64;93 for &s in word {94 let s = letter(s) as f64;95 up += (s.powi(dimension as i32) - 1.0).ln();96 down += s.ln();97 }98 if down == 0.0 {99 return dimension as f64;100 }101 up / down102}103104/// Builds the plane sieve the word spells as a raster: its side, then one byte a site, row by row, one where the site survives and zero where a level punched it out.105///106/// ```107/// let (side, cells) = mrlynum::sieve::raster(&mrlynum::sieve::odd_word(2));108/// assert_eq!(side, 15);109/// assert_eq!(cells.iter().filter(|&&b| b == 1).count(), 192);110/// ```111pub fn raster(word: &[u64]) -> (usize, Vec<u8>) {112 let mut side = 1usize;113 let mut sites = vec![1u8];114 for &s in word {115 let s = letter(s) as usize;116 let half = s / 2;117 let wide = side * s;118 let mut next = vec![0u8; wide * wide];119 for row in 0..side {120 for col in 0..side {121 if sites[row * side + col] == 0 {122 continue;123 }124 for i in 0..s {125 for j in 0..s {126 if i == half && j == half {127 continue;128 }129 next[(row * s + i) * wide + col * s + j] = 1;130 }131 }132 }133 }134 side = wide;135 sites = next;136 }137 (side, sites)138}139140/// Lists every puncture the word makes in the given dimension: its corner along each axis and then its side, all in units of the word's finest cell, so a level-one hole is the widest block in the list.141///142/// ```143/// let holes = mrlynum::sieve::punctures(&mrlynum::sieve::odd_word(2), 2);144/// assert_eq!(holes.len() / 3, 9);145/// assert_eq!(&holes[..3], &[5, 5, 5]);146/// ```147pub fn punctures(word: &[u64], dimension: u32) -> Vec<u64> {148 let axes = dimension as usize;149 let mut scale = side(word) as u64;150 let mut alive: Vec<u64> = vec![0; axes];151 let mut out: Vec<u64> = Vec::new();152 let mut digits = vec![0u64; axes];153 for (place, &s) in word.iter().enumerate() {154 let s = letter(s);155 let half = s / 2;156 let last = place + 1 == word.len();157 scale /= s;158 let count = s.pow(dimension);159 let mut next: Vec<u64> = Vec::new();160 for origin in alive.chunks(axes) {161 for axis in origin {162 out.push((axis * s + half) * scale);163 }164 out.push(scale);165 if last {166 continue;167 }168 for index in 0..count {169 let mut rest = index;170 let mut centre = true;171 for axis in (0..axes).rev() {172 digits[axis] = rest % s;173 centre &= digits[axis] == half;174 rest /= s;175 }176 if centre {177 continue;178 }179 for axis in 0..axes {180 next.push(origin[axis] * s + digits[axis]);181 }182 }183 }184 if !last {185 alive = next;186 }187 }188 out189}190191// THE SOLID LIMIT192193const SPLIT: u64 = 1001;194195fn odd_tail(start: u64, power: i32) -> f64 {196 let n = start as f64;197 let s = f64::from(power);198 let p = n.powi(-power);199 p * n / (2.0 * (s - 1.0)) + p / 2.0 + s * p / (6.0 * n)200 - s * (s + 1.0) * (s + 2.0) * p / (90.0 * n * n * n)201}202203/// Returns the limit of the solid Wallis sieve's surviving volume, the product of one minus n to the minus three over the odd n from three, in closed form.204///205/// Write m = 2k + 1 and factor m^3 - 1 = (m - 1)(m - w)(m - w^2) at w = exp(2 pi i / 3), the three cube roots of one; dividing by m^3 = 8 (k + 1/2)^3 turns the k-th factor into k (k + (1 - w)/2) (k + (1 - w^2)/2) / (k + 1/2)^3, so the product is a ratio of three rising shifts over one repeated three times.206///207/// The three shifts a1 = 0, a2 = (1 - w)/2 and a3 = (1 - w^2)/2 sum to (2 - w - w^2)/2 = 3/2, exactly three times the shift b = 1/2 below, and prod_{k=1..N} (k + a) = Gamma(N + 1 + a) / Gamma(1 + a) with Gamma(N + 1 + a) / Gamma(N + 1 + b) ~ N^(a - b), so the N-dependent factor tends to N^(a1 + a2 + a3 - 3b) = 1 and the limit is Gamma(1 + b)^3 / (Gamma(1 + a1) Gamma(1 + a2) Gamma(1 + a3)).208///209/// That is Gamma(3/2)^3 / (Gamma(1) Gamma((3 - w)/2) Gamma((3 - w^2)/2)) = pi^(3/2) / (8 |Gamma(7/4 - i sqrt(3)/4)|^2), the two gamma values being conjugates; the same Weierstrass product gives the identity prod_{k >= 1} (1 - a^3/k^3) = 1 / (Gamma(1 - a) Gamma(1 - a w) Gamma(1 - a w^2)) the ratio form of this limit leans on, at a = 1/2 against prod_{n >= 2} (1 - n^-3) = cosh(pi sqrt(3)/2) / (3 pi).210///211/// The value is read from the logarithm rather than from that gamma ratio, which cancels four digits away: log P = sum over odd n from 3 to 999 of log(1 - n^-3), taken from its smallest term, minus sum_{j = 1..3} (1/j) sum_{n odd >= 1001} n^(-3j), each inner sum by Euler-Maclaurin through the fourth Bernoulli term. Every piece carries one sign, so nothing cancels and the exponential lands within a few ulps.212///213/// ```214/// assert!((mrlynum::sieve::solid_limit() - 0.948_815_485_719_679_6).abs() < 4e-16);215/// ```216pub fn solid_limit() -> f64 {217 let mut log = 0.0f64;218 let mut n = SPLIT - 2;219 while n >= 3 {220 log += (-(n as f64).powi(-3)).ln_1p();221 n -= 2;222 }223 for j in 1..=3i32 {224 log -= odd_tail(SPLIT, 3 * j) / f64::from(j);225 }226 log.exp()227}228229/// Returns the limit the word's schedule walks to in the given dimension, when the word names a schedule at all.230///231/// A run of consecutive odd letters from a, of two letters or more, names the schedule a, a + 2, a + 4 and on: its letters are strictly increasing, sum s_k^(-d) converges, and the limit is positive, pi over four in the plane and the closed form above in the cube divided by the head the run skips. Two or more copies of one letter name the fixed-ratio schedule, whose limit is zero; that is what repetition costs, and an alternating word such as 3, 5, 3, 5 changes at every step and still loses the whole measure. Every other word, and every word of one letter, names no schedule and reads None.232///233/// ```234/// assert_eq!(mrlynum::sieve::limit(&mrlynum::sieve::odd_word(3), 2), Some(std::f64::consts::FRAC_PI_4));235/// assert_eq!(mrlynum::sieve::limit(&mrlynum::sieve::flat_word(3, 3), 2), Some(0.0));236/// assert!((mrlynum::sieve::limit(&[5, 7, 9], 2).unwrap() - 0.883_572_933_822_129_3).abs() < 1e-15);237/// assert_eq!(mrlynum::sieve::limit(&[3, 7, 11], 2), None);238/// ```239pub fn limit(word: &[u64], dimension: u32) -> Option<f64> {240 if word.len() < 2 {241 return None;242 }243 let first = letter(word[0]);244 if word.iter().all(|&s| letter(s) == first) {245 return Some(0.0);246 }247 if !word248 .iter()249 .enumerate()250 .all(|(k, &s)| s == first + 2 * k as u64)251 {252 return None;253 }254 let whole = match dimension {255 2 => PLANE_LIMIT,256 3 => solid_limit(),257 _ => return None,258 };259 let head = odd_word(((first - 3) / 2) as usize);260 Some(whole / ratio(&head, dimension))261}262263#[cfg(test)]264mod tests {265 use super::*;266 use crate::series::wallis;267 use std::f64::consts::PI;268269 fn truncated(start: u64, step: u64, stop: u64, scale: f64) -> f64 {270 let count = (stop - start) / step + 1;271 let mut sum = 0.0f64;272 for place in (0..count).rev() {273 let n = (start + place * step) as f64;274 sum += (-scale / n.powi(3)).ln_1p();275 }276 sum.exp()277 }278279 fn brute(word: &[u64]) -> usize {280 raster(word).1.iter().filter(|&&b| b == 1).count()281 }282283 #[test]284 fn the_odd_word_walks_the_sides_of_the_classical_sieve() {285 let word = odd_word(4);286 let sides: Vec<u128> = (1..=4).map(|n| side(&word[..n])).collect();287 assert_eq!(sides, vec![3, 15, 105, 945]);288 let counts: Vec<u128> = (1..=4).map(|n| cells(&word[..n], 2)).collect();289 assert_eq!(counts, vec![8, 192, 9216, 737_280]);290 }291292 #[test]293 fn the_raster_counts_what_the_product_promises() {294 let word = odd_word(3);295 for n in 1..=3 {296 assert_eq!(brute(&word[..n]) as u128, cells(&word[..n], 2), "level {n}");297 }298 for n in 1..=4 {299 let flat = flat_word(3, n);300 assert_eq!(brute(&flat) as u128, cells(&flat, 2), "carpet level {n}");301 }302 let five = flat_word(5, 2);303 assert_eq!(brute(&five) as u128, cells(&five, 2));304 }305306 #[test]307 fn the_raster_is_the_ratio_times_the_area() {308 for n in 1..=3 {309 let word = odd_word(n);310 let (wide, sites) = raster(&word);311 let lit = sites.iter().filter(|&&b| b == 1).count() as f64;312 assert!(313 (lit / (wide * wide) as f64 - ratio(&word, 2)).abs() < 1e-12,314 "level {n}"315 );316 }317 }318319 #[test]320 fn the_punctures_fill_the_gap_the_survivors_leave() {321 for dimension in 2..=3u32 {322 for n in 1..=3 {323 let word = odd_word(n);324 let axes = dimension as usize;325 let list = punctures(&word, dimension);326 assert_eq!(list.len() / (axes + 1), holes(&word, dimension) as usize);327 let volume: u128 = list328 .chunks(axes + 1)329 .map(|hole| u128::from(hole[axes]).pow(dimension))330 .sum();331 let whole = side(&word).pow(dimension);332 assert_eq!(333 volume + cells(&word, dimension),334 whole,335 "{dimension}d level {n}"336 );337 }338 }339 }340341 #[test]342 fn the_plane_ratio_is_the_wallis_product_the_series_walks() {343 for n in 1..=40 {344 assert!(345 (ratio(&odd_word(n), 2) - wallis(n)).abs() < 1e-15,346 "level {n}"347 );348 }349 assert!((ratio(&odd_word(200_000), 2) - PLANE_LIMIT).abs() < 1e-5);350 }351352 #[test]353 fn the_carpet_schedule_freezes_the_ratio_and_the_exponent() {354 for n in 1..=8 {355 let word = flat_word(3, n);356 assert!((ratio(&word, 2) - (8.0f64 / 9.0).powi(n as i32)).abs() < 1e-15);357 assert!((exponent(&word, 2) - 8f64.ln() / 3f64.ln()).abs() < 1e-12);358 }359 let long = odd_word(4_000);360 assert!(exponent(&long, 2) > 1.999_5 && exponent(&long, 2) < 2.0);361 }362363 #[test]364 fn the_solid_limit_holds_the_true_constant_to_four_ulps() {365 assert!((solid_limit() - 0.948_815_485_719_679_6).abs() < 4e-16);366 }367368 #[test]369 fn the_solid_limit_matches_its_truncated_product_to_one_part_in_1e14() {370 let product = truncated(3, 2, 8_000_001, 1.0);371 assert!((solid_limit() - product).abs() < 1e-14, "{product}");372 }373374 #[test]375 fn the_solid_limit_splits_into_the_cosh_and_the_even_product() {376 let all = truncated(2, 1, 2_000_000, 1.0);377 let cosh = (PI * 3f64.sqrt() / 2.0).cosh() / (3.0 * PI);378 assert!((all - cosh).abs() < 1e-12, "{all} {cosh}");379 let even = truncated(1, 1, 1_000_000, 0.125);380 assert!((solid_limit() - cosh / even).abs() < 1e-12);381 }382383 #[test]384 fn the_alternating_word_changes_at_every_step_and_freezes() {385 let pinned = 192f64.ln() / 15f64.ln();386 for pairs in 1..=4 {387 let word: Vec<u64> = [3u64, 5].iter().cycle().take(2 * pairs).copied().collect();388 assert!((exponent(&word, 2) - pinned).abs() < 1e-12, "{pairs}");389 }390 assert_eq!(format!("{pinned:.6}"), "1.941432");391 assert!(ratio(&[3, 5, 3, 5], 2) < ratio(&odd_word(4), 2));392 }393394 #[test]395 fn the_limit_names_a_schedule_or_says_it_cannot() {396 let plane = limit(&[5, 7, 9], 2).unwrap();397 assert!((plane - FRAC_PI_4 / (8.0 / 9.0)).abs() < 1e-15, "{plane}");398 assert_eq!(format!("{plane:.6}"), "0.883573");399 let solid = limit(&[5, 7, 9], 3).unwrap();400 assert!(401 (solid - solid_limit() / (26.0 / 27.0)).abs() < 1e-15,402 "{solid}"403 );404 assert_eq!(limit(&odd_word(4), 2), Some(PLANE_LIMIT));405 assert_eq!(limit(&flat_word(7, 3), 2), Some(0.0));406 assert_eq!(limit(&[3, 7, 11], 2), None);407 assert_eq!(limit(&[5], 2), None);408 assert_eq!(limit(&[3, 5, 3, 5], 2), None);409 assert_eq!(limit(&odd_word(4), 4), None);410 }411412 #[test]413 fn the_plane_ratio_at_two_million_factors_rounds_to_nine_digits() {414 let read = ratio(&odd_word(2_000_000), 2);415 assert_eq!(format!("{read:.9}"), "0.785398262");416 assert!(read > PLANE_LIMIT);417 }418419 #[test]420 #[should_panic(expected = "a letter is an odd side from three")]421 fn an_even_letter_has_no_centre_to_punch() {422 let _ = cells(&[4], 2);423 }424}