qualia_core_db/solvers/special_functions/
zeta.rs

1//! Riemann zeta function `ζ(s)` for real `s > 1`, via Euler–Maclaurin acceleration:
2//! sum the first `N−1` terms directly, then add the integral tail and Bernoulli
3//! corrections. `None` for `s ≤ 1` (the series there needs analytic continuation —
4//! out of this function's honest domain).
5
6const N: u32 = 10;
7// Bernoulli numbers B₂, B₄, B₆, B₈ and the factorials (2k)!.
8const BERNOULLI: [f64; 4] = [1.0 / 6.0, -1.0 / 30.0, 1.0 / 42.0, -1.0 / 30.0];
9const FACT_2K: [f64; 4] = [2.0, 24.0, 720.0, 40320.0];
10
11/// `ζ(s)` for real `s > 1`. `None` otherwise (domain).
12pub fn zeta(s: f64) -> Option<f64> {
13    if s <= 1.0 {
14        return None;
15    }
16    let nf = N as f64;
17    // Σ_{n=1}^{N-1} n^{-s}
18    let mut sum: f64 = (1..N).map(|n| (n as f64).powf(-s)).sum();
19    // ∫ tail + ½ f(N)
20    sum += nf.powf(1.0 - s) / (s - 1.0);
21    sum += 0.5 * nf.powf(-s);
22    // Bernoulli corrections: Σ_k B_{2k}/(2k)! · (s)_{2k-1} · N^{-s-2k+1}
23    for k in 1..=4usize {
24        let mut poch = 1.0; // rising factorial (s)(s+1)…(s+2k-2), length 2k-1
25        for j in 0..(2 * k - 1) {
26            poch *= s + j as f64;
27        }
28        sum += BERNOULLI[k - 1] / FACT_2K[k - 1] * poch * nf.powf(-s - (2 * k - 1) as f64);
29    }
30    Some(sum)
31}
32
33#[cfg(test)]
34mod tests {
35    use super::*;
36    const TOL: f64 = 1e-9;
37
38    #[test]
39    fn exact_even_values() {
40        // ζ(2) = π²/6, ζ(4) = π⁴/90
41        let pi = core::f64::consts::PI;
42        assert!((zeta(2.0).unwrap() - pi * pi / 6.0).abs() < TOL);
43        assert!((zeta(4.0).unwrap() - pi.powi(4) / 90.0).abs() < TOL);
44    }
45
46    #[test]
47    fn apery_and_large_s() {
48        // ζ(3) = Apéry's constant
49        assert!((zeta(3.0).unwrap() - 1.202_056_903_159_594_3).abs() < 1e-8);
50        // ζ(s) → 1 as s → ∞
51        assert!((zeta(30.0).unwrap() - 1.0).abs() < 1e-9);
52    }
53
54    #[test]
55    fn domain_fails_closed() {
56        assert!(zeta(1.0).is_none());
57        assert!(zeta(0.5).is_none());
58    }
59}
qualia_core_db/solvers/special_functions/zeta.rs

qualia_core_db/solvers/special_functions/
zeta.rs
