ladder.rs

1use crate::design::elements;2use crate::zeta::{raise, Complex};3use mrlycore::errors::{value_error, Result};45/// The relative rounding allowance the double-precision ladder charges against the scale it carries.6///7/// The truncation bound is the proved one, carried by the same recursion that carries the value; this constant is a measured allowance for the arithmetic itself and is not a proof.8pub const ROUNDING: f64 = 1e-13;910const PEEL_TARGET: f64 = 100.0;11const SHIFT_START: f64 = 6.0;12const SHIFT_STEP: f64 = 6.0;13const SHIFT_CAP: f64 = 400.0;14const CUT_START: usize = 8;15const CUT_STEP: usize = 6;16const CUT_CAP: usize = 120;17const RATIO_CAP: f64 = 0.9;1819/// A digit design: the base `q`, the digit set `F` its elements are written with, and the peel depth `P` its ladder starts at.20///21/// The elements are the whole numbers whose base-`q` digits all lie in `F` with a nonzero leading digit, and the object is their Dirichlet series `zeta_F(s) = sum n^(-s)`, of abscissa `alpha = log_q k` for `k = card F`.22/// With `E_j` the sum over the elements of exactly `j` digits and `G_P = sum_(j >= P) E_j`, splitting an element on its last digit gives `(1 - k q^(-w)) G_P(w) = E_P(w) + sum_(l >= 1) binom(-w, l) q^(-w-l) gamma_l G_P(w+l)` with `gamma_l = sum_(a in F) a^l`, and `zeta_F = D_(P-1) + G_P` for the finite Dirichlet polynomial `D_(P-1)` over the elements below `q^(P-1)`.23/// The `l`-series has ratio `max F / q^P`, so the peel depth buys the convergence and the small quantity `G_P` is carried directly, never as a difference of two large ones.24#[derive(Clone, Debug, PartialEq)]25pub struct Design {26    base: u64,27    digits: Vec<u64>,28    peel: usize,29    lead: usize,30    largest: u64,31    low: Vec<f64>,32    mid: Vec<f64>,33    gamma: Vec<f64>,34}3536impl Design {37    /// Builds a design on the base and the digit set, choosing the peel depth, or an error when the base or the digits are out of range.38    ///39    /// ```40    /// let design = mrlynum::ladder::Design::new(3, &[0, 1]).unwrap();41    /// assert_eq!(design.peel(), 7);42    /// assert!((design.abscissa() - 0.630_929_753_571_457).abs() < 1e-14);43    /// ```44    pub fn new(base: u64, digits: &[u64]) -> Result<Design> {45        let mut set: Vec<u64> = digits.to_vec();46        set.sort_unstable();47        set.dedup();48        if set.len() < 2 {49            return value_error("a design needs at least two digits.");50        }51        let lead = set.iter().filter(|&&d| d > 0).count();52        if lead == 0 {53            return value_error("a design needs a nonzero digit.");54        }55        let size = set.len() as f64;56        let mut peel = 2usize;57        while lead as f64 * size.powi(peel as i32) <= PEEL_TARGET {58            peel += 1;59        }60        Design::with_peel(base, &set, peel)61    }6263    /// Builds a design at an explicit peel depth, at least two, or an error when the base, the digits or the depth are out of range.64    ///65    /// ```66    /// let shallow = mrlynum::ladder::Design::with_peel(2, &[0, 1], 4).unwrap();67    /// assert_eq!(shallow.peel(), 4);68    /// ```69    pub fn with_peel(base: u64, digits: &[u64], peel: usize) -> Result<Design> {70        let mut set: Vec<u64> = digits.to_vec();71        set.sort_unstable();72        set.dedup();73        if base < 2 {74            return value_error(format!("base {base} is under two."));75        }76        if set.len() < 2 {77            return value_error("a design needs at least two digits.");78        }79        if set.iter().any(|&d| d >= base) {80            return value_error(format!("a digit sits outside the base {base}."));81        }82        let lead = set.iter().filter(|&&d| d > 0).count();83        if lead == 0 {84            return value_error("a design needs a nonzero digit.");85        }86        if peel < 2 {87            return value_error(format!("peel depth {peel} is under two."));88        }89        if (base as f64).powi(peel as i32) > 9.007_199_254_740_992e15 {90            return value_error(format!("peel depth {peel} overruns the exact integers."));91        }92        let below = elements(base, &set, peel - 1);93        let through = elements(base, &set, peel);94        let low: Vec<f64> = below.iter().map(|&n| (n as f64).ln()).collect();95        let mid: Vec<f64> = through[below.len()..]96            .iter()97            .map(|&n| (n as f64).ln())98            .collect();99        let largest = *set.last().unwrap();100        let gamma: Vec<f64> = (0..=CUT_CAP)101            .map(|l| set.iter().map(|&a| (a as f64).powi(l as i32)).sum())102            .collect();103        Ok(Design {104            base,105            digits: set,106            peel,107            lead,108            largest,109            low,110            mid,111            gamma,112        })113    }114115    /// Returns the base.116    pub fn base(&self) -> u64 {117        self.base118    }119120    /// Returns the digit set, ascending.121    pub fn digits(&self) -> &[u64] {122        &self.digits123    }124125    /// Returns the peel depth.126    pub fn peel(&self) -> usize {127        self.peel128    }129130    /// Returns the abscissa `alpha = log_q k`.131    pub fn abscissa(&self) -> f64 {132        (self.digits.len() as f64).ln() / (self.base as f64).ln()133    }134135    /// Returns the pole spacing `2 pi / log q`.136    pub fn period(&self) -> f64 {137        2.0 * std::f64::consts::PI / (self.base as f64).ln()138    }139140    /// Returns the pole `s_(m,j) = alpha - m + 2 pi i j / log q`.141    pub fn pole(&self, m: usize, j: i64) -> Complex {142        Complex::new(self.abscissa() - m as f64, self.period() * j as f64)143    }144}145146// LADDER147148fn poly(logs: &[f64], w: Complex) -> Complex {149    let mut acc = Complex::new(0.0, 0.0);150    for &l in logs {151        acc = acc + ((-w) * l).exp();152    }153    acc154}155156fn poly_scale(logs: &[f64], sigma: f64) -> f64 {157    logs.iter().map(|&l| (-sigma * l).exp()).sum()158}159160fn tail(design: &Design, sigma: f64) -> f64 {161    let ratio = design.digits.len() as f64 * (design.base as f64).powf(-sigma);162    if ratio >= 1.0 {163        return f64::INFINITY;164    }165    design.lead as f64 * ratio.powi(design.peel as i32 - 1) / (1.0 - ratio)166}167168fn log_tail(design: &Design, sigma: f64) -> f64 {169    let ratio = design.digits.len() as f64 * (design.base as f64).powf(-sigma);170    if ratio >= 1.0 {171        return f64::INFINITY;172    }173    (design.lead as f64).ln() + (design.peel as f64 - 1.0) * ratio.ln() - (-ratio).ln_1p()174}175176fn log_binomial(top: f64, pick: usize) -> f64 {177    (1..=pick)178        .map(|i| ((top - pick as f64 + i as f64) / i as f64).ln())179        .sum()180}181182fn cut(design: &Design, w: Complex, depth: usize) -> f64 {183    let (span, sigma) = (w.abs(), w.re);184    let base = design.base as f64;185    let largest = design.largest as f64;186    let ratio = ((span + depth as f64 + 1.0) / (depth as f64 + 2.0)).max(1.0) * largest187        / base.powi(design.peel as i32);188    if ratio >= RATIO_CAP {189        return f64::INFINITY;190    }191    let rest = log_tail(design, sigma + depth as f64 + 1.0);192    if rest.is_infinite() {193        return f64::INFINITY;194    }195    let log = log_binomial(span + depth as f64, depth + 1)196        + (-sigma - depth as f64 - 1.0) * base.ln()197        + (design.digits.len() as f64).ln()198        + (depth as f64 + 1.0) * largest.ln()199        + rest;200    log.exp() / (1.0 - ratio)201}202203struct Rung {204    value: Complex,205    bound: f64,206    scale: f64,207}208209fn rung(design: &Design, s: Complex, shift: f64, depth: usize, whole: bool) -> Rung {210    let base = design.base as f64;211    let size = design.digits.len() as f64;212    let levels = (shift - s.re).ceil().max(0.0) as usize;213    if !whole && levels == 0 {214        let front = Complex::new(1.0, 0.0) - raise(base, -s) * size;215        let rest = front.abs() * tail(design, s.re);216        return Rung {217            value: Complex::new(0.0, 0.0),218            bound: rest,219            scale: rest,220        };221    }222    let mut value = vec![Complex::new(0.0, 0.0); levels + depth + 2];223    let mut bound = vec![0.0f64; levels + depth + 2];224    let mut scale = vec![0.0f64; levels + depth + 2];225    for j in levels..levels + depth + 2 {226        bound[j] = tail(design, s.re + j as f64);227        scale[j] = bound[j];228    }229    for j in (0..levels).rev() {230        let w = s + j as f64;231        let mut acc = poly(&design.mid, w);232        let mut err = cut(design, w, depth);233        let mut wide = poly_scale(&design.mid, w.re) + err;234        let mut coefficient = Complex::new(1.0, 0.0);235        for l in 1..=depth {236            coefficient = coefficient * ((-w) - (l as f64 - 1.0)) * (1.0 / l as f64);237            let weight = raise(base, (-w) - l as f64) * coefficient * design.gamma[l];238            acc = acc + weight * value[j + l];239            err += weight.abs() * bound[j + l];240            wide += weight.abs() * scale[j + l];241        }242        if !whole && j == 0 {243            return Rung {244                value: acc,245                bound: err,246                scale: wide,247            };248        }249        let den = Complex::new(1.0, 0.0) - raise(base, -w) * size;250        value[j] = acc / den;251        bound[j] = err / den.abs();252        scale[j] = wide / den.abs();253    }254    Rung {255        value: value[0],256        bound: bound[0],257        scale: scale[0],258    }259}260261fn crossing(design: &Design, s: Complex, whole: bool) -> Option<(i64, i64)> {262    let turn = s.im / design.period();263    let j = turn.round();264    if (turn - j).abs() > 1e-9 {265        return None;266    }267    let drop = design.abscissa() - s.re;268    let m = drop.round();269    let first = if whole { 0.0 } else { 1.0 };270    if (drop - m).abs() > 1e-9 || m < first {271        return None;272    }273    Some((m as i64, j as i64))274}275276fn tune(277    design: &Design,278    s: Complex,279    tolerance: f64,280    whole: bool,281    head: f64,282) -> Result<(Complex, f64)> {283    if let Some((m, j)) = crossing(design, s, whole) {284        return value_error(format!(285            "the ladder walks through the pole s_({m},{j}) = alpha - {m} + 2 pi i {j} / log q and cannot read that point."286        ));287    }288    let mut shift = SHIFT_START;289    let mut depth = CUT_START;290    let mut best = f64::INFINITY;291    while shift <= SHIFT_CAP && depth <= CUT_CAP {292        let step = rung(design, s, shift, depth, whole);293        let carried = step.bound + ROUNDING * (step.scale + head);294        if carried < tolerance {295            return Ok((step.value, carried));296        }297        best = best.min(carried);298        shift += SHIFT_STEP;299        depth += CUT_STEP;300    }301    value_error(format!(302        "tolerance {tolerance:e} is out of reach at s = {} + {}i, best bound {best:e}.",303        s.re, s.im304    ))305}306307// READINGS308309/// Returns `zeta_F(s)` and the bound it is known to, or an error when the tolerance is out of reach or the point sits on the blind lattice.310///311/// The bound is the propagated truncation bound, which is proved, plus [`ROUNDING`] times the scale the ladder carries, which is a measured allowance; the acceptance test charges every term the returned bound carries, the Dirichlet polynomial's own scale included, so the returned bound is never above the tolerance asked.312/// The recursion walks `w = s + j` upward and divides by `1 - k q^(-w)`, so it cannot read `s = alpha - m + 2 pi i j / log q` for a whole `m >= 0`: at those points it raises and names the pole. On a full digit set that lattice is `s = 1, 0, -1, -2, ...`, where `zeta_F` is `zeta` and only `s = 1` is singular.313///314/// ```315/// let design = mrlynum::ladder::Design::new(2, &[0, 1]).unwrap();316/// let s = mrlynum::zeta::Complex::new(2.0, 0.0);317/// let (value, bound) = mrlynum::ladder::zeta(&design, s, 1e-10).unwrap();318/// assert!((value.re - std::f64::consts::PI * std::f64::consts::PI / 6.0).abs() < bound);319/// ```320pub fn zeta(design: &Design, s: Complex, tolerance: f64) -> Result<(Complex, f64)> {321    let (value, bound) = tune(design, s, tolerance, true, poly_scale(&design.low, s.re))?;322    Ok((poly(&design.low, s) + value, bound))323}324325/// Returns the Lyndon cofactor `Z(s) = zeta_F(s) (1 - k q^(-s))` and the bound it is known to, or an error when the tolerance is out of reach or the point sits on the blind lattice.326///327/// The cofactor is the ladder numerator over the finite polynomial, so it is analytic on `Re s > alpha - 1` and a zero census on it needs no pole-free strip.328/// The blind lattice here is `s = alpha - m + 2 pi i j / log q` for a whole `m >= 1`: the numerator itself is read at `m = 0`, which is what makes the residue available.329pub fn cofactor(design: &Design, s: Complex, tolerance: f64) -> Result<(Complex, f64)> {330    let base = design.base as f64;331    let size = design.digits.len() as f64;332    let front = Complex::new(1.0, 0.0) - raise(base, -s) * size;333    let head = front.abs() * poly_scale(&design.low, s.re);334    let (value, bound) = tune(design, s, tolerance, false, head)?;335    Ok((front * poly(&design.low, s) + value, bound))336}337338/// Returns the residue of `zeta_F` at `s_(m,j) = alpha - m + 2 pi i j / log q` and the bound it is known to, or an error when the tolerance is out of reach.339///340/// At `m = 0` the factor `1 - k q^(-s)` has derivative `log q`, so the residue is the ladder numerator over `log q`; the column below comes off the same recursion, `(1 - q^m) R_m = sum_(l = 1)^m binom(-s_(m,j), l) q^(-s_(m,j)-l) gamma_l R_(m-l)`.341/// The head term `(1 - k q^(-s)) D_(P-1)(s)` is dropped, since the factor vanishes at the pole; in double precision the factor is not exactly zero, and the residue of the dropped term, about `1e-15` at base 10, is not charged separately and sits inside the carried scale.342pub fn residue(design: &Design, m: usize, j: i64, tolerance: f64) -> Result<(Complex, f64)> {343    let base = design.base as f64;344    let (head, err) = tune(design, design.pole(0, j), tolerance, false, 0.0)?;345    let mut values = vec![head * (1.0 / base.ln())];346    let mut bounds = vec![err / base.ln()];347    for level in 1..=m {348        let w = design.pole(level, j);349        let mut acc = Complex::new(0.0, 0.0);350        let mut err = 0.0;351        let mut coefficient = Complex::new(1.0, 0.0);352        for l in 1..=level {353            coefficient = coefficient * ((-w) - (l as f64 - 1.0)) * (1.0 / l as f64);354            let weight = raise(base, (-w) - l as f64) * coefficient * design.gamma[l];355            acc = acc + weight * values[level - l];356            err += weight.abs() * (bounds[level - l] + ROUNDING * values[level - l].abs());357        }358        let den = 1.0 - base.powi(level as i32);359        values.push(acc * (1.0 / den));360        bounds.push(err / den.abs());361    }362    Ok((values[m], bounds[m]))363}364365#[cfg(test)]366mod tests {367    use super::*;368    use crate::zeta::Line;369370    const PINNED: f64 = 1e-10;371372    fn near(value: Complex, re: f64, im: f64, bound: f64) {373        assert!(374            (value.re - re).abs() < bound && (value.im - im).abs() < bound,375            "{value:?} is not {re} + {im}i inside {bound:e}"376        );377    }378379    #[test]380    fn the_full_base_two_design_is_the_zeta_function() {381        let design = Design::new(2, &[0, 1]).unwrap();382        assert_eq!(design.peel(), 7);383        assert_eq!(design.abscissa(), 1.0);384        let pinned = [385            (Complex::new(2.0, 0.0), 1.644_934_066_848_226_4, 0.0),386            (387                Complex::new(0.3, 40.0),388                0.748_775_209_504_225_8,389                -1.440_885_440_634_440_5,390            ),391            (392                Complex::new(-1.0, 2.0),393                0.168_915_669_770_834_4,394                -0.070_515_988_908_254_42,395            ),396        ];397        for (s, re, im) in pinned {398            let (value, bound) = zeta(&design, s, 1e-8).unwrap();399            near(value, re, im, bound);400        }401    }402403    #[test]404    fn the_ladder_meets_the_critical_line_engine() {405        let design = Design::new(2, &[0, 1]).unwrap();406        let line = Line::new();407        for t in [12.0, 18.5] {408            let (value, bound) = zeta(&design, Complex::new(0.5, t), PINNED).unwrap();409            let want = line.maclaurin(t);410            assert!(411                (value - want).abs() < bound + 1e-9,412                "t {t} {value:?} {want:?}"413            );414        }415    }416417    #[test]418    fn the_residue_is_one_at_the_abscissa_and_zero_above_it() {419        let design = Design::new(2, &[0, 1]).unwrap();420        let (one, bound) = residue(&design, 0, 0, PINNED).unwrap();421        near(one, 1.0, 0.0, bound);422        let (off, bound) = residue(&design, 0, 1, PINNED).unwrap();423        near(off, 0.0, 0.0, bound);424    }425426    #[test]427    fn the_base_three_residues_meet_the_harvest_midpoints() {428        let pinned = [429            (430                vec![0u64, 1],431                [432                    (0.799_023_642_655_835, 0.0),433                    (0.231_891_517_689_918_2, -0.501_067_414_481_068_3),434                    (0.003_963_750_587_374_229, -0.033_400_689_771_018_48),435                ],436            ),437            (438                vec![0, 2],439                [440                    (0.515_977_601_099_115_1, 0.0),441                    (0.135_292_875_345_483_3, 0.329_874_091_244_466_17),442                    (-0.021_699_535_708_644_655, -0.000_946_804_031_498_436_2),443                ],444            ),445        ];446        for (digits, wanted) in pinned {447            let design = Design::new(3, &digits).unwrap();448            for (j, (re, im)) in wanted.iter().enumerate() {449                let (value, bound) = residue(&design, 0, j as i64, PINNED).unwrap();450                near(value, *re, *im, bound);451            }452        }453    }454455    #[test]456    fn the_second_digit_set_scales_the_first() {457        let one = Design::new(3, &[0, 1]).unwrap();458        let two = Design::new(3, &[0, 2]).unwrap();459        for s in [Complex::new(2.0, 3.0), Complex::new(0.8, 23.0)] {460            let (a, ba) = zeta(&two, s, PINNED).unwrap();461            let (b, bb) = zeta(&one, s, PINNED).unwrap();462            assert!((a - raise(2.0, -s) * b).abs() < ba + bb);463        }464    }465466    #[test]467    fn the_base_ten_design_misses_its_last_digit() {468        let design = Design::new(10, &[0, 1, 2, 3, 4, 5, 6, 7, 8]).unwrap();469        assert_eq!(design.peel(), 2);470        let (value, bound) = zeta(&design, Complex::new(2.0, 0.0), PINNED).unwrap();471        near(value, 1.623_924_927_084_640_5, 0.0, bound);472        let (value, bound) = zeta(&design, Complex::new(1.1, 3.0), PINNED).unwrap();473        near(474            value,475            0.618_774_761_630_085,476            -0.025_931_036_103_864_726,477            bound,478        );479        let (value, bound) = residue(&design, 0, 0, PINNED).unwrap();480        near(value, 1.022_344_896_589_484_1, 0.0, bound);481    }482483    #[test]484    fn the_cofactor_is_the_series_times_the_lyndon_factor() {485        let design = Design::new(2, &[0, 1]).unwrap();486        for s in [Complex::new(2.0, 0.0), Complex::new(0.3, 40.0)] {487            let (value, bound) = cofactor(&design, s, 1e-8).unwrap();488            let (series, other) = zeta(&design, s, 1e-8).unwrap();489            let front = Complex::new(1.0, 0.0) - raise(2.0, -s) * 2.0;490            assert!((value - front * series).abs() < bound + front.abs() * other);491        }492        let (value, bound) = cofactor(&design, Complex::new(2.0, 0.0), PINNED).unwrap();493        near(value, 0.822_467_033_424_113_2, 0.0, bound);494    }495496    #[test]497    fn the_bound_survives_a_change_of_peel_depth() {498        for s in [Complex::new(2.0, 0.0), Complex::new(0.3, 40.0)] {499            let shallow = Design::with_peel(2, &[0, 1], 4).unwrap();500            let deep = Design::with_peel(2, &[0, 1], 8).unwrap();501            let (a, ba) = zeta(&shallow, s, 1e-8).unwrap();502            let (b, bb) = zeta(&deep, s, 1e-8).unwrap();503            assert!((a - b).abs() < ba + bb, "{s:?} {a:?} {b:?} {ba:e} {bb:e}");504        }505    }506507    #[test]508    fn the_column_below_a_pole_is_the_contour_average_around_it() {509        let design = Design::new(3, &[0, 1]).unwrap();510        let (want, bound) = residue(&design, 1, 1, PINNED).unwrap();511        let centre = design.pole(1, 1);512        let (radius, nodes) = (0.05, 8);513        let mut got = Complex::new(0.0, 0.0);514        let mut carried = 0.0;515        for k in 0..nodes {516            let turn = Complex::turn(std::f64::consts::TAU * k as f64 / nodes as f64);517            let (value, other) = zeta(&design, centre + turn * radius, 1e-9).unwrap();518            got = got + turn * radius * value * (1.0 / nodes as f64);519            carried += radius * other / nodes as f64;520        }521        assert!(522            (got - want).abs() < bound + carried + 1e-9,523            "{got:?} {want:?}"524        );525    }526527    #[test]528    fn the_ladder_raises_when_the_tolerance_is_out_of_reach() {529        let design = Design::new(2, &[0, 1]).unwrap();530        assert!(zeta(&design, Complex::new(-1.0, 2.0), 1e-11).is_err());531        assert!(zeta(&design, Complex::new(-1.0, 2.0), 1e-9).is_ok());532    }533534    #[test]535    fn the_constructor_rejects_a_design_it_cannot_peel() {536        assert!(Design::new(5, &[3]).is_err());537        assert!(Design::new(5, &[0]).is_err());538        assert!(Design::new(1, &[0, 1]).is_err());539        assert!(Design::with_peel(5, &[0, 7], 3).is_err());540        assert!(Design::with_peel(2, &[0, 1], 1).is_err());541    }542543    #[test]544    fn the_returned_bound_never_exceeds_the_tolerance_asked() {545        let design = Design::new(2, &[0, 1]).unwrap();546        for s in [547            Complex::new(-2.0, 2.0),548            Complex::new(-6.0, 2.0),549            Complex::new(0.3, 40.0),550            Complex::new(2.0, 0.0),551        ] {552            for tolerance in [1e-1, 1e-5, 1e-8] {553                if let Ok((_, bound)) = zeta(&design, s, tolerance) {554                    assert!(bound < tolerance, "zeta {s:?} {tolerance:e} {bound:e}");555                }556                if let Ok((_, bound)) = cofactor(&design, s, tolerance) {557                    assert!(bound < tolerance, "cofactor {s:?} {tolerance:e} {bound:e}");558                }559            }560        }561    }562563    #[test]564    fn the_cofactor_bound_carries_the_lyndon_factor_on_an_empty_rung() {565        let design = Design::new(2, &[0, 1]).unwrap();566        let s = Complex::new(6.0, std::f64::consts::PI / 2f64.ln());567        let front = Complex::new(1.0, 0.0) - raise(2.0, -s) * 2.0;568        assert!((front.abs() - 1.031_25).abs() < 1e-12);569        let (_, bound) = cofactor(&design, s, 1e-7).unwrap();570        let (_, other) = zeta(&design, s, 1e-7).unwrap();571        assert!(bound > front.abs() * 9.6e-10, "{bound:e}");572        assert!(bound > other * 1.03, "{bound:e} {other:e}");573    }574575    #[test]576    fn the_ladder_names_the_pole_it_walks_through() {577        let design = Design::new(2, &[0, 1]).unwrap();578        for m in 0..3 {579            let s = design.pole(m, 0);580            let message = zeta(&design, s, 1e-7).unwrap_err().to_string();581            assert!(message.contains("walks through the pole"), "{message}");582        }583        assert!(zeta(&design, Complex::new(-1.0, 1e-6), 1e-7).is_ok());584    }585}