formulas.rs

1use crate::classics::primes;2use crate::factor::mobius_sieve;3use crate::series::harmonic;45pub use crate::series::li;67// SIEVES89fn primal(limit: usize) -> Vec<bool> {10    let mut flag = vec![limit >= 2; limit + 1];11    for slot in flag.iter_mut().take(2.min(limit + 1)) {12        *slot = false;13    }14    let mut p = 2;15    while p * p <= limit {16        if flag[p] {17            let mut m = p * p;18            while m <= limit {19                flag[m] = false;20                m += p;21            }22        }23        p += 1;24    }25    flag26}2728fn pairs(flag: &[bool], number: usize) -> usize {29    (2..=number / 2)30        .filter(|&p| flag[p] && flag[number - p])31        .count()32}3334// CONSTANTS BY SUMMATION3536/// Returns the Wallis product taken to n paired factors, four k squared over four k squared less one, walking to pi over two.37///38/// ```39/// assert!((mrlynum::formulas::wallis(1) - 4.0 / 3.0).abs() < 1e-15);40/// ```41pub fn wallis(n: usize) -> f64 {42    let mut out = 1.0;43    for k in 1..=n {44        let square = ((2 * k) * (2 * k)) as f64;45        out *= square / (square - 1.0);46    }47    out48}4950/// Returns the Leibniz alternating sum of the odd reciprocals over n terms, walking to pi over four.51///52/// ```53/// assert!((mrlynum::formulas::leibniz(2) - 2.0 / 3.0).abs() < 1e-15);54/// ```55pub fn leibniz(n: usize) -> f64 {56    let mut out = 0.0;57    for k in 0..n {58        let term = 1.0 / (2 * k + 1) as f64;59        out += if k.is_multiple_of(2) { term } else { -term };60    }61    out62}6364/// Returns the Basel sum of the reciprocal squares over n terms, walking to pi squared over six.65///66/// ```67/// assert!((mrlynum::formulas::basel(3) - 49.0 / 36.0).abs() < 1e-15);68/// ```69pub fn basel(n: usize) -> f64 {70    (1..=n).map(|k| 1.0 / (k * k) as f64).sum()71}7273/// Returns the harmonic sum of n terms less the logarithm of n, walking to the Euler-Mascheroni constant.74///75/// ```76/// assert!((mrlynum::formulas::euler_gamma_partial(1) - 1.0).abs() < 1e-15);77/// ```78pub fn euler_gamma_partial(n: usize) -> f64 {79    if n == 0 {80        return 0.0;81    }82    harmonic(n) - (n as f64).ln()83}8485/// Returns one plus one over n raised to the n, walking to the natural base.86///87/// ```88/// assert!((mrlynum::formulas::e_partial(1) - 2.0).abs() < 1e-15);89/// ```90pub fn e_partial(n: usize) -> f64 {91    if n == 0 {92        return 1.0;93    }94    (1.0 + 1.0 / n as f64).powf(n as f64)95}9697// THE PRIMES COUNTED9899/// Returns the count of primes at or below n.100///101/// ```102/// assert_eq!(mrlynum::formulas::prime_count(100), 25);103/// ```104pub fn prime_count(n: usize) -> usize {105    primes(n).len()106}107108/// Returns the count of unordered pairs of primes summing to the number, zero below four.109///110/// ```111/// assert_eq!(mrlynum::formulas::goldbach(100), 6);112/// ```113pub fn goldbach(number: usize) -> usize {114    if number < 4 {115        return 0;116    }117    pairs(&primal(number), number)118}119120/// Returns the count of prime pairs at every even number from four up to the top, one entry per even number.121///122/// ```123/// assert_eq!(mrlynum::formulas::goldbach_record(10), vec![1, 1, 1, 2]);124/// ```125pub fn goldbach_record(top: usize) -> Vec<usize> {126    if top < 4 {127        return Vec::new();128    }129    let flag = primal(top);130    (2..=top / 2).map(|k| pairs(&flag, 2 * k)).collect()131}132133/// Returns the Mertens function at n, the Mobius values of one through n summed.134///135/// ```136/// assert_eq!(mrlynum::formulas::mertens(100), 1);137/// ```138pub fn mertens(n: usize) -> i64 {139    mobius_sieve(n).iter().skip(1).map(|&v| i64::from(v)).sum()140}141142#[cfg(test)]143mod tests {144    use super::*;145    use crate::prime::splits;146    use std::f64::consts::{E, PI};147148    #[test]149    fn the_small_partials_are_their_exact_fractions() {150        assert!((wallis(1) - 4.0 / 3.0).abs() < 1e-15);151        assert!((wallis(2) - 64.0 / 45.0).abs() < 1e-15);152        assert!((leibniz(1) - 1.0).abs() < 1e-15);153        assert!((leibniz(2) - 2.0 / 3.0).abs() < 1e-15);154        assert!((basel(3) - 49.0 / 36.0).abs() < 1e-15);155        assert!((euler_gamma_partial(1) - 1.0).abs() < 1e-15);156        assert!((e_partial(2) - 2.25).abs() < 1e-15);157    }158159    #[test]160    fn the_five_partials_walk_to_their_constants() {161        assert!((wallis(200_000) - PI / 2.0).abs() < 1e-5);162        assert!((leibniz(200_000) - PI / 4.0).abs() < 1e-5);163        assert!((basel(200_000) - PI * PI / 6.0).abs() < 1e-4);164        assert!((euler_gamma_partial(200_000) - crate::series::EULER).abs() < 1e-5);165        assert!((e_partial(200_000) - E).abs() < 1e-4);166    }167168    #[test]169    fn li_matches_its_own_series_at_two() {170        assert!((li(2.0) - 1.045_163_780_117_493).abs() < 1e-12);171    }172173    #[test]174    fn the_prime_count_is_the_sieve() {175        assert_eq!(prime_count(1), 0);176        assert_eq!(prime_count(100), 25);177        assert_eq!(prime_count(1000), 168);178        assert_eq!(prime_count(10_000), 1229);179    }180181    #[test]182    fn goldbach_counts_the_prime_pairs_the_long_way() {183        assert_eq!(goldbach(4), 1);184        assert_eq!(goldbach(100), 6);185        assert_eq!(goldbach(1000), 28);186        for even in (4..=600).step_by(2) {187            assert_eq!(goldbach(even), splits(even).len());188        }189    }190191    #[test]192    fn the_goldbach_record_never_reaches_zero_below_ten_thousand() {193        let record = goldbach_record(10_000);194        assert_eq!(record.len(), 4999);195        assert_eq!(record[0], 1);196        assert_eq!(record[498], 28);197        assert_eq!(record.iter().copied().min(), Some(1));198    }199200    #[test]201    fn mertens_sums_the_mobius_values() {202        assert_eq!(mertens(0), 0);203        assert_eq!(mertens(1), 1);204        assert_eq!(mertens(100), 1);205        assert_eq!(mertens(1000), 2);206        assert_eq!(mertens(10_000), -23);207    }208}