series.rs

1use crate::classics::primes;2use crate::factor::mobius_sieve;3use std::f64::consts::PI;45/// The Basel constant, pi squared over six, the value zeta takes at two.6pub const BASEL: f64 = PI * PI / 6.0;78/// The visible density, six over pi squared, the share of lattice pairs that are coprime.9pub const VISIBLE: f64 = 6.0 / (PI * PI);1011/// The Catalan constant, the value the Dirichlet beta function takes at two.12pub const CATALAN: f64 = 0.915_965_594_177_219;1314/// The Apery constant, the value zeta takes at three.15pub const APERY: f64 = 1.202_056_903_159_594;1617/// The Euler constant, the limit of the harmonic sum less the logarithm.18pub const EULER: f64 = 0.577_215_664_901_532_9;1920/// Returns the logarithmic integral of a positive x by the Ramanujan series, the smooth count of the primes below x.21///22/// ```23/// assert!((mrlynum::series::li(1_000_000.0) - 78_627.549).abs() < 1e-3);24/// ```25pub fn li(x: f64) -> f64 {26    let log = x.ln();27    let mut sum = 0.0f64;28    let mut term = 1.0f64;29    let mut odds = 0.0f64;30    for n in 1..200u32 {31        term *= log / n as f64;32        if n > 1 {33            term /= 2.0;34        }35        if !n.is_multiple_of(2) {36            odds += 1.0 / n as f64;37        }38        let piece = term * odds;39        sum += if n.is_multiple_of(2) { -piece } else { piece };40        if piece.abs() < 1e-17 * sum.abs() {41            break;42        }43    }44    EULER + log.abs().ln() + x.sqrt() * sum45}4647/// Returns the partial harmonic sum, the reciprocals of one through the term count.48pub fn harmonic(terms: usize) -> f64 {49    (1..=terms).map(|k| 1.0 / k as f64).sum()50}5152/// Returns the zeta value above one: the partial sum closed by its Euler-Maclaurin tail.53///54/// Panics at an s of one or below, where the sum does not converge.55pub fn zeta(s: f64, terms: usize) -> f64 {56    assert!(s > 1.0, "zeta needs an s above one");57    let mut sum = 0.0;58    for k in 1..=terms {59        sum += (k as f64).powf(-s);60    }61    let n = terms as f64;62    sum + n.powf(1.0 - s) / (s - 1.0) - n.powf(-s) / 2.0 + s * n.powf(-s - 1.0) / 12.063}6465/// Returns the Euler product of zeta, one over one minus p to the minus s over the primes up to the limit.66pub fn euler_product(s: f64, limit: usize) -> f64 {67    primes(limit)68        .iter()69        .map(|&p| 1.0 / (1.0 - (p as f64).powf(-s)))70        .product()71}7273/// Returns the Dirichlet beta value, the alternating odd-denominator sum averaged over its last two partial sums.74pub fn beta(s: f64, terms: usize) -> f64 {75    let mut sum = 0.0;76    let mut previous = 0.0;77    for k in 0..terms {78        previous = sum;79        let term = ((2 * k + 1) as f64).powf(-s);80        sum += if k.is_multiple_of(2) { term } else { -term };81    }82    0.5 * (previous + sum)83}8485/// Returns the Dirichlet lambda value, one minus two to the minus s times zeta.86pub fn lambda(s: f64, terms: usize) -> f64 {87    (1.0 - 2f64.powf(-s)) * zeta(s, terms)88}8990/// Returns the mod-eight rhythm of the number, the discriminant minus-eight character: one on one and three, minus one on five and seven, zero on the evens.91///92/// ```93/// assert_eq!(mrlynum::series::chi8(3), 1);94/// assert_eq!(mrlynum::series::chi8(5), -1);95/// ```96pub fn chi8(number: usize) -> i8 {97    [0, 1, 0, 1, 0, -1, 0, -1][number % 8]98}99100/// Returns the mod-four rhythm of the number: zero, one, zero, minus one.101///102/// ```103/// assert_eq!(mrlynum::series::chi4(7), -1);104/// ```105pub fn chi4(number: usize) -> i8 {106    [0, 1, 0, -1][number % 4]107}108109/// Returns the mod-three rhythm of the number: zero, one, minus one.110pub fn chi3(number: usize) -> i8 {111    [0, 1, -1][number % 3]112}113114/// Returns the L-series partial sum with a periodic rhythm painted on the terms.115pub fn dirichlet(s: f64, rhythm: &[i8], terms: usize) -> f64 {116    if rhythm.is_empty() {117        return 0.0;118    }119    let mut sum = 0.0;120    for n in 1..=terms {121        let paint = rhythm[n % rhythm.len()];122        if paint != 0 {123            sum += f64::from(paint) * (n as f64).powf(-s);124        }125    }126    sum127}128129/// Counts the lattice points of the dimension-cube of the limit whose coordinates share no divisor, by Mobius inversion.130///131/// Panics at a zero dimension, and wraps once the limit to the dimension passes a signed hundred and twenty-eight bits.132pub fn visible(limit: usize, dimension: u32) -> u128 {133    assert!(dimension > 0, "visible needs a dimension above zero");134    let mu = mobius_sieve(limit);135    let mut total: i128 = 0;136    for (k, &value) in mu.iter().enumerate().skip(1) {137        if value == 0 {138            continue;139        }140        let block = (limit / k) as i128;141        total += i128::from(value) * block.pow(dimension);142    }143    total as u128144}145146/// Returns the Wallis product of one minus one over the odd squares, walking to pi over four.147pub fn wallis(factors: usize) -> f64 {148    let mut out = 1.0;149    for n in 1..=factors {150        let odd = (2 * n + 1) as f64;151        out *= 1.0 - 1.0 / (odd * odd);152    }153    out154}155156fn reduce(num: i128, den: i128) -> (i128, i128) {157    let (mut a, mut b) = (num.abs(), den.abs());158    while b != 0 {159        (a, b) = (b, a % b);160    }161    let sign = if den < 0 { -1 } else { 1 };162    if a == 0 {163        return (0, 1);164    }165    (sign * num / a, sign * den / a)166}167168fn binomial(n: usize, k: usize) -> i128 {169    let mut out: i128 = 1;170    for i in 0..k {171        out = out * (n - i) as i128 / (i + 1) as i128;172    }173    out174}175176/// Builds the first Bernoulli numbers as exact reduced fractions on the minus one half convention.177///178/// Panics past a count of thirty-two, where the exact fractions overflow a signed hundred and twenty-eight bits.179pub fn bernoulli(count: usize) -> Vec<(i128, i128)> {180    assert!(count <= 32, "the exact fractions overflow past thirty-two");181    let mut out: Vec<(i128, i128)> = Vec::with_capacity(count);182    for m in 0..count {183        if m == 0 {184            out.push((1, 1));185            continue;186        }187        let mut sum = (0i128, 1i128);188        for (j, &(num, den)) in out.iter().enumerate() {189            let weight = binomial(m + 1, j) * num;190            sum = reduce(sum.0 * den + weight * sum.1, sum.1 * den);191        }192        out.push(reduce(-sum.0, sum.1 * (m + 1) as i128));193    }194    out195}196197#[cfg(test)]198mod tests {199    use super::*;200    use crate::lattice::coprime_pairs;201202    #[test]203    fn zeta_meets_the_basel_the_apery_and_the_quartic_sum() {204        assert!((zeta(2.0, 10_000) - BASEL).abs() < 1e-9);205        assert!((zeta(3.0, 10_000) - APERY).abs() < 1e-12);206        assert!((zeta(4.0, 10_000) - PI.powi(4) / 90.0).abs() < 1e-9);207    }208209    #[test]210    fn the_euler_product_meets_the_basel_and_the_apery_sum() {211        assert!((euler_product(2.0, 100_000) - BASEL).abs() < 1e-5);212        assert!((euler_product(3.0, 100_000) - APERY).abs() < 1e-6);213    }214215    #[test]216    fn beta_walks_to_catalan_and_to_a_quarter_turn() {217        assert!((beta(2.0, 1_000_000) - CATALAN).abs() < 1e-9);218        assert!((beta(1.0, 1_000_000) - PI / 4.0).abs() < 1e-6);219    }220221    #[test]222    fn half_lambda_is_the_grid_fluctuation_constant() {223        assert!((0.5 * lambda(4.0, 10_000) - PI.powi(4) / 192.0).abs() < 1e-9);224    }225226    #[test]227    fn the_mod_four_rhythm_paints_the_beta_series() {228        let painted = dirichlet(2.0, &[0, 1, 0, -1], 1_000_000);229        assert!((painted - beta(2.0, 1_000_000)).abs() < 1e-9);230    }231232    #[test]233    fn the_mod_three_rhythm_paints_the_l_series() {234        let painted = dirichlet(1.0, &[0, 1, -1], 1_000_000);235        assert!((painted - PI / (3.0 * 3f64.sqrt())).abs() < 1e-5);236    }237238    #[test]239    fn visible_counts_the_coprime_pairs_of_a_window() {240        for n in 1..=2_000 {241            assert_eq!(visible(n, 2), u128::from(coprime_pairs(n)), "window {n}");242        }243    }244245    #[test]246    fn the_visible_density_is_one_over_zeta() {247        let flat = visible(10_000, 2) as f64 / 1e8;248        assert!((flat - VISIBLE).abs() < 1e-3);249        let cube = visible(1_000, 3) as f64 / 1e9;250        assert!((cube - 1.0 / zeta(3.0, 100_000)).abs() < 1e-2);251    }252253    #[test]254    fn the_harmonic_walk_leaves_the_euler_mascheroni_gap() {255        let gap = harmonic(1_000_000) - 1_000_000f64.ln();256        assert!((gap - 0.577_215_664_9).abs() < 1e-5);257    }258259    #[test]260    fn bernoulli_pins_the_known_fractions() {261        let list = bernoulli(32);262        assert_eq!(list[0], (1, 1));263        assert_eq!(list[1], (-1, 2));264        assert_eq!(list[2], (1, 6));265        assert_eq!(list[3], (0, 1));266        assert_eq!(list[4], (-1, 30));267        assert_eq!(list[12], (-691, 2730));268    }269270    #[test]271    fn every_odd_bernoulli_past_the_first_is_zero() {272        for (index, &(num, den)) in bernoulli(32).iter().enumerate().skip(3) {273            if !index.is_multiple_of(2) {274                assert_eq!((num, den), (0, 1), "index {index}");275            }276        }277    }278279    #[test]280    #[should_panic(expected = "visible needs a dimension above zero")]281    fn visible_refuses_a_zero_dimension() {282        let _ = visible(10, 0);283    }284285    #[test]286    #[should_panic(expected = "the exact fractions overflow past thirty-two")]287    fn bernoulli_refuses_a_count_past_thirty_two() {288        let _ = bernoulli(33);289    }290291    #[test]292    #[should_panic(expected = "zeta needs an s above one")]293    fn zeta_refuses_an_s_of_one() {294        let _ = zeta(1.0, 10);295    }296297    #[test]298    fn wallis_walks_to_a_quarter_turn() {299        assert!((wallis(1_000_000) - PI / 4.0).abs() < 1e-6);300    }301302    #[test]303    fn the_rhythms_repeat_over_their_full_period() {304        let four: Vec<i8> = (0..8).map(chi4).collect();305        assert_eq!(four, vec![0, 1, 0, -1, 0, 1, 0, -1]);306        let three: Vec<i8> = (0..6).map(chi3).collect();307        assert_eq!(three, vec![0, 1, -1, 0, 1, -1]);308        let eight: Vec<i8> = (0..16).map(chi8).collect();309        assert_eq!(310            eight,311            vec![0, 1, 0, 1, 0, -1, 0, -1, 0, 1, 0, 1, 0, -1, 0, -1]312        );313    }314315    #[test]316    fn the_mod_eight_rhythm_multiplies_across_the_odd_numbers() {317        for a in (1..=99usize).step_by(2) {318            for b in (1..=99usize).step_by(2) {319                assert_eq!(chi8(a) * chi8(b), chi8(a * b), "{a} {b}");320            }321        }322    }323324    #[test]325    fn li_pins_the_smooth_prime_counts() {326        assert!((li(2.0) - 1.045_163_780_1).abs() < 1e-9);327        assert!((li(1_000.0) - 177.609_657_990_2).abs() < 1e-8);328        assert!((li(10_000.0) - 1_246.137_215_9).abs() < 1e-6);329        assert!((li(100_000.0) - 9_629.809_001_1).abs() < 1e-6);330        assert!((li(1_000_000.0) - 78_627.549_159_5).abs() < 1e-6);331    }332}