dissection.rs

1use crate::num::design;2use crate::num::factor::{gcd, mobius_sieve, totient};3use crate::num::prime::flags;4use std::f64::consts::PI;56/// The bar region A sets on the certificate exponent, `alpha_1 < 1/5`: the whole `l^1` mass against the minor-arc `x^(4/5)`.7pub const BAR_A: f64 = 0.2;8/// The bar region B sets, `alpha_1 < 1/4`: the hybrid `l^1` mass against the `d^(-1/2)` decay.9pub const BAR_B: f64 = 0.25;10/// The constant `gamma'` of the digit-uniform chain, `(2/pi)(gamma + log(8/pi))` rounded up at seven decimals.11pub const GAMMA: f64 = 0.9625229;12/// The base the chain is scanned to; from here the closed-form cap `1 + sqrt(2 (2/pi) log base + 0.97)` keeps it below the bar.13pub const CAP_BASE: u64 = 1272;14/// The least base whose every one-missing-digit set carries a window certificate below `1/5`, `lab/py/prime-dissection`, verb `window`.15pub const WINDOW_WALL: u64 = 301;16/// The least base whose every one-missing-digit set carries a per-digit shifted-grid certificate below `1/5`, `lab/py/digit-transform-norms`, verbs `fifth` and `fifthbelow`; base 114 missing 56 is certified above.17pub const DIGIT_WALL: u64 = 115;18/// The least base with some missing digit certified below `1/5`, `lab/py/digit-transform-norms`, verb `fifth`.19pub const FIRST_BELOW: u64 = 65;20/// The sets below [`DIGIT_WALL`] certified below `1/5`, as `(base, missing digit)`, printed by `lab/py/digit-transform-norms`, verb `fifth`.21pub const CERTIFIED: [(u64, u64); 1] = [(65, 0)];2223/// The four regions the dissection cuts the frequencies `a/y` into by their Dirichlet fraction `l/d` at `Q = y^(3/5)` and height `h = |a d - l y|`.24#[derive(Clone, Copy, Debug, PartialEq, Eq)]25pub enum Region {26    /// The minor arcs, `d >= y^(2/5)`.27    A,28    /// The middle, `d < y^(2/5)` and `max(d, h) >= Z`.29    B,30    /// Near a fraction whose denominator has a prime outside the base, `d < Z` and `h < Z`.31    C1,32    /// Near a fraction whose denominator divides a power of the base, `d < Z` and `h < Z`.33    C2,34}3536/// Returns the last continued-fraction convergent `l/d` of `a/y` whose denominator is at most the cap, the Dirichlet fraction of the dissection.37///38/// ```39/// assert_eq!(mrlyrs::num::dissection::fraction(333, 1000, 63), (1, 3));40/// ```41pub fn fraction(a: u64, y: u64, cap: u64) -> (u64, u64) {42    let (mut p0, mut q0, mut p1, mut q1) = (0u64, 1u64, 1u64, 0u64);43    let (mut n, mut d) = (a, y);44    let mut best = (0, 1);45    while d > 0 {46        let c = n / d;47        let (p2, q2) = (c * p1 + p0, c * q1 + q0);48        if q2 > cap {49            break;50        }51        (p0, q0, p1, q1) = (p1, q1, p2, q2);52        best = (p2, q2);53        (n, d) = (d, n - c * d);54    }55    best56}5758/// Returns `Q = floor(y^(3/5))`, the largest denominator the dissection admits, exact in integers.59///60/// ```61/// assert_eq!(mrlyrs::num::dissection::cap(1000), 63);62/// ```63pub fn cap(y: u64) -> u64 {64    let cube = u128::from(y).pow(3);65    let mut c = (y as f64).powf(0.6) as u64;66    while u128::from(c + 1).pow(5) <= cube {67        c += 1;68    }69    while c > 0 && u128::from(c).pow(5) > cube {70        c -= 1;71    }72    c73}7475fn smooth(base: u64, mut d: u64) -> bool {76    loop {77        let g = gcd(u128::from(d), u128::from(base)) as u64;78        if g == 1 {79            return d == 1;80        }81        while d.is_multiple_of(g) {82            d /= g;83        }84    }85}8687/// Returns the region of every frequency `a/y`, `a < y = base^level`, at the cut `Z`, the fraction taken by [`fraction`] at [`cap`].88pub fn regions(base: u64, level: u32, z: u64) -> Vec<Region> {89    let y = base.pow(level);90    let top = cap(y);91    let square = u128::from(y).pow(2);92    (0..y)93        .map(|a| {94            let (l, d) = fraction(a, y, top);95            let h = (i128::from(a) * i128::from(d) - i128::from(l) * i128::from(y)).unsigned_abs();96            if u128::from(d).pow(5) >= square {97                Region::A98            } else if d < z && h < u128::from(z) {99                if smooth(base, d) {100                    Region::C2101                } else {102                    Region::C1103                }104            } else {105                Region::B106            }107        })108        .collect()109}110111fn modulus(base: u64, missing: &[u64], digits: &[u64], a: u64, y: u64) -> f64 {112    let turn =113        |k: u64| PI * ((u128::from(k) * u128::from(a) % (2 * u128::from(y))) as f64) / y as f64;114    let term = |k: u64| (2.0 * turn(k)).sin_cos();115    let (im, re) = if missing.len() < digits.len() {116        let kernel = if a == 0 {117            base as f64118        } else {119            turn(base).sin() / turn(1).sin()120        };121        let (sin, cos) = turn(base - 1).sin_cos();122        missing123            .iter()124            .fold((kernel * sin, kernel * cos), |(im, re), &e| {125                let (s, c) = term(e);126                (im - s, re - c)127            })128    } else {129        digits.iter().fold((0.0, 0.0), |(im, re), &d| {130            let (s, c) = term(d);131            (im + s, re + c)132        })133    };134    re.hypot(im)135}136137fn rise(base: u64, digits: &[u64], row: &[f64]) -> Vec<f64> {138    let missing: Vec<u64> = (0..base).filter(|d| !digits.contains(d)).collect();139    let below = row.len() as u64;140    let y = below * base;141    (0..y)142        .map(|a| modulus(base, &missing, digits, a, y) * row[(a % below) as usize])143        .collect()144}145146/// Returns `|hat F_level(a/base^level)|` at every `a < base^level`, built one digit at a time from `hat F_j(t) = hat F(t) hat F_(j-1)(base t)`.147///148/// ```149/// let row = mrlyrs::num::dissection::weights(3, &[0, 1], 1);150/// assert_eq!(row.iter().map(|w| format!("{w:.3}")).collect::<Vec<_>>(), ["2.000", "1.000", "1.000"]);151/// ```152pub fn weights(base: u64, digits: &[u64], level: u32) -> Vec<f64> {153    (0..level).fold(vec![1.0], |row, _| rise(base, digits, &row))154}155156/// Returns the unshifted masses `c_j = sum_(a < base^j) |hat F_j(a/base^j)|` for `j = 0..=level`, the `l^1` mass region A pays, `c_0 = 1`.157pub fn masses(base: u64, digits: &[u64], level: u32) -> Vec<f64> {158    let mut row = vec![1.0];159    let mut out = vec![1.0];160    for _ in 0..level {161        row = rise(base, digits, &row);162        out.push(row.iter().sum());163    }164    out165}166167/// Returns the `l^1` exponent the top two masses read, `log_base(c_j/(fill c_(j-1)))`: a reading of the growth region A pays, never a certificate.168pub fn reading(base: u64, fill: usize, masses: &[f64]) -> f64 {169    let n = masses.len();170    if n < 2 {171        return f64::NAN;172    }173    (masses[n - 1] / (fill as f64 * masses[n - 2])).ln() / (base as f64).ln()174}175176fn cubic(base: u64, z: f64) -> f64 {177    let slope = 2.0 / PI;178    (z - 1.0).powi(3)179        - (slope * (base as f64).ln() * z180            + GAMMA * (z - 1.0)181            + slope * (z - 1.0).powi(2) / (base as f64 * z - 1.0))182}183184/// Returns the root `z > 1` of the digit-uniform chain at one missing digit, `(z-1)^3 = (2/pi)(log base) z + gamma'(z-1) + (2/pi)(z-1)^2/(base z - 1)`, by bisection.185pub fn chain_root(base: u64) -> f64 {186    let (mut low, mut high) = (1.0, 11.0);187    while cubic(base, high) < 0.0 {188        high *= 2.0;189    }190    for _ in 0..200 {191        let mid = 0.5 * (low + high);192        if cubic(base, mid) < 0.0 {193            low = mid;194        } else {195            high = mid;196        }197    }198    high199}200201/// Returns the chain's certificate exponent at one missing digit, `alpha_1 = log_base(z base/(base - 1))`, the same at every missing digit.202pub fn chain_exponent(base: u64) -> f64 {203    let q = base as f64;204    (chain_root(base) * q / (q - 1.0)).ln() / q.ln()205}206207/// Returns the chain's margin at the bar, the cubic cleared of denominators at `w = base^(1/5)(1 - 1/base)`, positive exactly when the root sits below `w` and `alpha_1 < 1/5`.208pub fn chain_margin(base: u64) -> f64 {209    let q = base as f64;210    cubic(base, q.powf(BAR_A) * (1.0 - 1.0 / q))211}212213/// Returns the chain's wall, one past the last base up to [`CAP_BASE`] whose margin is not positive.214pub fn chain_wall() -> u64 {215    (3..=CAP_BASE)216        .rev()217        .find(|&q| chain_margin(q) <= 0.0)218        .map_or(3, |q| q + 1)219}220221/// Returns how the theorem reaches the set missing one digit: `proof` from the chain's wall, `certificate` at every base from [`DIGIT_WALL`] below it and at the [`CERTIFIED`] sets, `none` elsewhere.222///223/// ```224/// use mrlyrs::num::dissection::reach;225/// assert_eq!([reach(584, 3), reach(115, 57), reach(65, 0), reach(65, 32)], ["proof", "certificate", "certificate", "none"]);226/// ```227pub fn reach(base: u64, missing: u64) -> &'static str {228    if base >= chain_wall() {229        "proof"230    } else if base >= DIGIT_WALL || CERTIFIED.contains(&(base, missing)) {231        "certificate"232    } else {233        "none"234    }235}236237/// Returns `kappa_F = (base/phi(base)) #{f in F : gcd(f, base) = 1}/fill`, the main-term constant of the prime count, as a reduced fraction.238///239/// ```240/// assert_eq!(mrlyrs::num::dissection::kappa(10, &[1, 2, 3, 4, 6, 7, 8, 9]), (5, 4));241/// ```242pub fn kappa(base: u64, digits: &[u64]) -> (u64, u64) {243    let units = digits244        .iter()245        .filter(|&&f| gcd(u128::from(f), u128::from(base)) == 1)246        .count() as u64;247    let top = base * units;248    let bottom = totient(base as usize) as u64 * digits.len() as u64;249    let g = gcd(u128::from(top), u128::from(bottom)).max(1) as u64;250    (top / g, bottom / g)251}252253/// Returns whether the digit set keeps two consecutive digits, the hypothesis region C1 reads.254pub fn consecutive(digits: &[u64]) -> bool {255    digits.iter().any(|d| digits.contains(&(d + 1)))256}257258/// The set's own sums read on a log grid of `x`: the mass `A_F(x)`, the meter `M_F(x)` and the prime count `psi_F(x) = sum Lambda(n)` over the elements up to `x`.259pub struct Tally {260    /// The log of `x` at every sample, uniform from the first element to the last.261    pub log_x: Vec<f64>,262    /// The mass `A_F(x)`, the count of elements up to `x`.263    pub count: Vec<u64>,264    /// The meter `M_F(x)`, the sum of `mu(n)` over the elements up to `x`.265    pub meter: Vec<i64>,266    /// The prime count `psi_F(x)`, the sum of `Lambda(n)` over the elements up to `x`.267    pub primes: Vec<f64>,268}269270/// Tallies the set below `base^level` on a log grid of the given size, with the Mobius values and the primes sieved to the span.271pub fn tally(base: u64, digits: &[u64], level: usize, samples: usize) -> Tally {272    let values = design::elements(base, digits, level);273    let top = base.pow(level as u32) as usize;274    let mu = mobius_sieve(top);275    let prime = flags(top);276    let mut powers: Vec<(u64, f64)> = Vec::new();277    for p in (2..=top).take_while(|p| p * p <= top).filter(|&p| prime[p]) {278        let mut power = p * p;279        while power <= top {280            powers.push((power as u64, (p as f64).ln()));281            power = match power.checked_mul(p) {282                Some(next) => next,283                None => break,284            };285        }286    }287    powers.sort_by_key(|row| row.0);288    let mut meter = Vec::with_capacity(values.len());289    let mut psi = Vec::with_capacity(values.len());290    let (mut m, mut s) = (0i64, 0.0f64);291    for &n in &values {292        m += i64::from(mu[n as usize]);293        if prime[n as usize] {294            s += (n as f64).ln();295        } else if let Ok(slot) = powers.binary_search_by_key(&n, |row| row.0) {296            s += powers[slot].1;297        }298        meter.push(m);299        psi.push(s);300    }301    let log_x = design::log_grid(&values, samples);302    let slots: Vec<usize> = log_x303        .iter()304        .map(|&t| values.partition_point(|&v| (v as f64).ln() <= t))305        .collect();306    Tally {307        count: slots.iter().map(|&k| k as u64).collect(),308        meter: slots309            .iter()310            .map(|&k| if k == 0 { 0 } else { meter[k - 1] })311            .collect(),312        primes: slots313            .iter()314            .map(|&k| if k == 0 { 0.0 } else { psi[k - 1] })315            .collect(),316        log_x,317    }318}319320#[cfg(test)]321mod tests {322    use super::*;323324    fn missing(base: u64, digit: u64) -> Vec<u64> {325        (0..base).filter(|&d| d != digit).collect()326    }327328    #[test]329    fn the_paper_figure_cut_holds_732_202_26_40() {330        let cut = regions(10, 3, 8);331        let count = |r: Region| cut.iter().filter(|&&c| c == r).count();332        assert_eq!(333            [334                count(Region::A),335                count(Region::B),336                count(Region::C1),337                count(Region::C2)338            ],339            [732, 202, 26, 40]340        );341    }342343    #[test]344    fn the_mass_ratios_match_the_wall_readings() {345        let ratio = |base: u64, digit: u64, level: u32| {346            let c = masses(base, &missing(base, digit), level);347            format!("{:.4}", c[level as usize] / c[level as usize - 1])348        };349        assert_eq!(ratio(33, 16, 4), "74.1654");350        assert_eq!(ratio(33, 0, 4), "70.5661");351        assert_eq!(ratio(17, 8, 5), "35.5234");352    }353354    #[test]355    fn the_two_digit_column_reads_log_two_over_log_three() {356        let c = masses(3, &[0, 1], 1);357        assert_eq!((reading(3, 2, &c) * 1e6).floor(), 630929.0);358    }359360    #[test]361    fn the_chain_clears_the_bar_from_584() {362        assert_eq!(chain_wall(), 584);363        assert_eq!(format!("{:.4e}", chain_margin(584)), "6.0170e-3");364        assert_eq!(format!("{:.4e}", chain_margin(583)), "-8.3138e-3");365        assert_eq!((chain_exponent(584) * 1e6).ceil(), 199983.0);366    }367368    #[test]369    fn the_tally_meets_the_count_by_hand() {370        let t = tally(10, &missing(10, 7), 6, 64);371        assert_eq!(372            (t.count[63], t.meter[63], format!("{:.4}", t.primes[63])),373            (531440, -9, "441976.4885".to_string())374        );375    }376}