zeta.rs

1use crate::classics::primes;2use crate::series::bernoulli;3use std::f64::consts::PI;4use std::ops::{Add, Div, Mul, Neg, Sub};56/// The t where the walk hands over from Euler-Maclaurin to Riemann-Siegel.7pub const JOIN: f64 = 20.0;8const SHIFT: usize = 10;9const TAIL: usize = 7;10const STEPS: usize = 10;11const STENCIL: f64 = 0.02;12const TOLERANCE: f64 = 1e-9;1314/// A complex number: a real and an imaginary part.15#[derive(Clone, Copy, Debug, Default, PartialEq)]16pub struct Complex {17    /// The real part.18    pub re: f64,19    /// The imaginary part.20    pub im: f64,21}2223impl Complex {24    /// Builds a complex number from its parts.25    pub const fn new(re: f64, im: f64) -> Complex {26        Complex { re, im }27    }28    /// Returns the modulus.29    pub fn abs(self) -> f64 {30        self.re.hypot(self.im)31    }32    /// Returns the principal argument.33    pub fn arg(self) -> f64 {34        self.im.atan2(self.re)35    }36    /// Returns the exponential.37    pub fn exp(self) -> Complex {38        let r = self.re.exp();39        Complex::new(r * self.im.cos(), r * self.im.sin())40    }41    /// Returns the principal logarithm.42    pub fn ln(self) -> Complex {43        Complex::new(self.abs().ln(), self.arg())44    }45    /// Returns a unit complex number at the given angle.46    pub fn turn(angle: f64) -> Complex {47        Complex::new(angle.cos(), angle.sin())48    }49}5051impl Add for Complex {52    type Output = Complex;53    fn add(self, other: Complex) -> Complex {54        Complex::new(self.re + other.re, self.im + other.im)55    }56}5758impl Add<f64> for Complex {59    type Output = Complex;60    fn add(self, other: f64) -> Complex {61        Complex::new(self.re + other, self.im)62    }63}6465impl Sub for Complex {66    type Output = Complex;67    fn sub(self, other: Complex) -> Complex {68        Complex::new(self.re - other.re, self.im - other.im)69    }70}7172impl Sub<f64> for Complex {73    type Output = Complex;74    fn sub(self, other: f64) -> Complex {75        Complex::new(self.re - other, self.im)76    }77}7879impl Mul for Complex {80    type Output = Complex;81    fn mul(self, other: Complex) -> Complex {82        Complex::new(83            self.re * other.re - self.im * other.im,84            self.re * other.im + self.im * other.re,85        )86    }87}8889impl Mul<f64> for Complex {90    type Output = Complex;91    fn mul(self, other: f64) -> Complex {92        Complex::new(self.re * other, self.im * other)93    }94}9596impl Div for Complex {97    type Output = Complex;98    fn div(self, other: Complex) -> Complex {99        let d = other.re * other.re + other.im * other.im;100        Complex::new(101            (self.re * other.re + self.im * other.im) / d,102            (self.im * other.re - self.re * other.im) / d,103        )104    }105}106107impl Neg for Complex {108    type Output = Complex;109    fn neg(self) -> Complex {110        Complex::new(-self.re, -self.im)111    }112}113114/// Returns a positive real base raised to a complex exponent.115pub fn raise(base: f64, exponent: Complex) -> Complex {116    (exponent * base.ln()).exp()117}118119/// Returns the Riemann-Siegel kernel, the cosine ratio that leads the remainder, in the form that stays finite at its removable points.120pub fn kernel(p: f64) -> f64 {121    let d = if p < 0.5 { 0.25 - p } else { p - 0.75 };122    let below = (2.0 * PI * d).sin();123    if below == 0.0 {124        return 0.5;125    }126    (PI * d * (1.0 + 2.0 * d)).sin() / below127}128129fn difference(order: usize, p: f64, step: f64) -> f64 {130    let mut sum = 0.0;131    let mut binomial = 1.0;132    for k in 0..=order {133        let offset = (order as f64 / 2.0 - k as f64) * step;134        let sign = if k.is_multiple_of(2) { 1.0 } else { -1.0 };135        sum += sign * binomial * kernel(p + offset);136        binomial = binomial * (order - k) as f64 / (k + 1) as f64;137    }138    sum / step.powi(order as i32)139}140141fn derivative(order: usize, p: f64) -> f64 {142    (4.0 * difference(order, p, STENCIL) - difference(order, p, 2.0 * STENCIL)) / 3.0143}144145/// Returns the first four Riemann-Siegel corrections at the fractional part p: the kernel and its derivatives by central differences with one Richardson step.146pub fn corrections(p: f64) -> [f64; 4] {147    let (pi2, pi4, pi6) = (PI * PI, PI.powi(4), PI.powi(6));148    [149        kernel(p),150        -derivative(3, p) / (96.0 * pi2),151        derivative(6, p) / (18_432.0 * pi4) + derivative(2, p) / (64.0 * pi2),152        -derivative(9, p) / (5_308_416.0 * pi6)153            - derivative(5, p) / (3_840.0 * pi4)154            - derivative(1, p) / (64.0 * pi2),155    ]156}157158/// The critical line: the Bernoulli numbers and the Euler-Maclaurin weights the two engines share, built once.159pub struct Line {160    bern: Vec<f64>,161    tail: Vec<f64>,162}163164impl Default for Line {165    fn default() -> Line {166        Line::new()167    }168}169170impl Line {171    /// Builds the line: the even Bernoulli numbers through the fourteenth and their Euler-Maclaurin weights.172    pub fn new() -> Line {173        let fractions = bernoulli(2 * TAIL + 1);174        let bern: Vec<f64> = (0..=TAIL)175            .map(|k| {176                let (num, den) = fractions[2 * k];177                num as f64 / den as f64178            })179            .collect();180        let mut factorial = 1.0;181        let tail = (0..=TAIL)182            .map(|k| {183                if k > 0 {184                    factorial *= ((2 * k - 1) * 2 * k) as f64;185                }186                bern[k] / factorial187            })188            .collect();189        Line { bern, tail }190    }191    /// Returns the Riemann-Siegel theta: the argument of gamma at one quarter plus i t over two, less t ln pi over two, by Stirling's series after a shift of ten.192    pub fn theta(&self, t: f64) -> f64 {193        let z = Complex::new(0.25, 0.5 * t);194        let w = z + SHIFT as f64;195        let mut arg = ((w - 0.5) * w.ln() - w).im;196        for k in 0..SHIFT {197            arg -= (z + k as f64).arg();198        }199        let inverse = Complex::new(1.0, 0.0) / w;200        let square = inverse * inverse;201        let mut power = inverse;202        for k in 1..=TAIL {203            arg += self.bern[k] / ((2 * k * (2 * k - 1)) as f64) * power.im;204            power = power * square;205        }206        arg - 0.5 * t * PI.ln()207    }208    /// Returns zeta at one half plus i t by the complex Euler-Maclaurin sum: t plus ten terms and seven Bernoulli corrections.209    pub fn maclaurin(&self, t: f64) -> Complex {210        let s = Complex::new(0.5, t);211        let count = t.abs() as usize + SHIFT;212        let mut sum = Complex::default();213        for k in 1..=count {214            sum = sum + raise(k as f64, -s);215        }216        let base = count as f64;217        let mut out = sum + raise(base, -s + 1.0) / (s - 1.0) - raise(base, -s) * 0.5;218        let mut rising = s;219        let mut power = raise(base, -s - 1.0);220        let square = 1.0 / (base * base);221        for k in 1..=TAIL {222            out = out + rising * power * self.tail[k];223            rising = rising * (s + (2 * k - 1) as f64) * (s + (2 * k) as f64);224            power = power * square;225        }226        out227    }228    /// Returns Z(t) from the Euler-Maclaurin value turned onto the real axis.229    pub fn exact(&self, t: f64) -> f64 {230        (Complex::turn(self.theta(t)) * self.maclaurin(t)).re231    }232    /// Returns Z(t) by the Riemann-Siegel formula: the main sum and the first four corrections.233    pub fn siegel(&self, t: f64) -> f64 {234        let a = (t / (2.0 * PI)).sqrt();235        let whole = a.floor();236        let theta = self.theta(t);237        let mut sum = 0.0;238        for k in 1..=whole as usize {239            let kf = k as f64;240            sum += (theta - t * kf.ln()).cos() / kf.sqrt();241        }242        let sign = if (whole as u64).is_multiple_of(2) {243            -1.0244        } else {245            1.0246        };247        let mut weight = 1.0 / a.sqrt();248        let mut rest = 0.0;249        for c in corrections(a - whole) {250            rest += c * weight;251            weight /= a;252        }253        2.0 * sum + sign * rest254    }255    /// Returns Z(t): Euler-Maclaurin below the join, Riemann-Siegel above.256    pub fn z(&self, t: f64) -> f64 {257        if t < JOIN {258            self.exact(t)259        } else {260            self.siegel(t)261        }262    }263    /// Returns zeta on the line and Z(t) together, from the engine that serves the t.264    pub fn point(&self, t: f64) -> (Complex, f64) {265        if t < JOIN {266            let value = self.maclaurin(t);267            (value, (Complex::turn(self.theta(t)) * value).re)268        } else {269            let z = self.siegel(t);270            (Complex::turn(-self.theta(t)) * z, z)271        }272    }273    /// Returns the largest gap between the two engines over the t range on a grid.274    pub fn seam(&self, t0: f64, t1: f64, steps: usize) -> f64 {275        (0..=steps)276            .map(|k| {277                let t = t0 + (t1 - t0) * k as f64 / steps as f64;278                (self.siegel(t) - self.exact(t)).abs()279            })280            .fold(0.0, f64::max)281    }282    /// Returns the n-th Gram point, where theta is n pi, by Newton from the right.283    pub fn gram(&self, n: i64) -> f64 {284        let target = n as f64 * PI;285        let mut t = 2.0 * PI * (n as f64 + 2.0).max(3.0);286        for _ in 0..100 {287            let step = (self.theta(t) - target) / (0.5 * (t / (2.0 * PI)).ln());288            t -= step;289            if step.abs() < 1e-12 {290                break;291            }292        }293        t294    }295    fn brackets(&self, limit: f64, count: usize, exact: bool) -> Vec<(f64, f64)> {296        let z = |t: f64| if exact { self.exact(t) } else { self.z(t) };297        let mut out = Vec::new();298        let mut n = -1;299        let mut left = self.gram(n);300        let mut previous = left;301        let mut before = z(left);302        'walk: while out.len() < count && left < limit {303            let right = self.gram(n + 1);304            for k in 1..=STEPS {305                let t = (left + (right - left) * k as f64 / STEPS as f64).min(limit);306                let now = z(t);307                if before * now < 0.0 {308                    out.push((previous, t));309                }310                previous = t;311                before = now;312                if t >= limit {313                    break 'walk;314                }315            }316            left = right;317            n += 1;318        }319        out.truncate(count);320        out321    }322    fn bisect(&self, (mut a, mut b): (f64, f64)) -> f64 {323        let mut fa = self.exact(a);324        while b - a > TOLERANCE {325            let mid = 0.5 * (a + b);326            let fm = self.exact(mid);327            if fa * fm <= 0.0 {328                b = mid;329            } else {330                a = mid;331                fa = fm;332            }333        }334        0.5 * (a + b)335    }336    /// Returns the first zeros on the line: sign changes of Z between Gram points, refined by bisection on Euler-Maclaurin to a billionth.337    pub fn zeros(&self, count: usize) -> Vec<f64> {338        self.brackets(f64::INFINITY, count, true)339            .into_iter()340            .map(|pair| self.bisect(pair))341            .collect()342    }343    /// Counts the zeros on the line below t.344    pub fn count(&self, t: f64) -> usize {345        self.brackets(t, usize::MAX, false).len()346    }347}348349/// Returns the Chebyshev staircase at every whole number from one to x: the sum of ln p over the prime powers up to each.350pub fn psi_stair(x: usize) -> Vec<f64> {351    let mut jumps = vec![0.0; x + 1];352    for p in primes(x) {353        let mut q = p;354        while q <= x {355            jumps[q] += (p as f64).ln();356            q *= p;357        }358    }359    let mut sum = 0.0;360    jumps[1..]361        .iter()362        .map(|jump| {363            sum += jump;364            sum365        })366        .collect()367}368369/// Returns the von Mangoldt explicit formula at x over the zeros at the given ordinates and their mirrors: x less the sum of x to the rho over rho, less ln two pi, less half the ln of one minus x to the minus two.370pub fn psi_formula(x: f64, gammas: &[f64]) -> f64 {371    let log = x.ln();372    let waves: f64 = gammas373        .iter()374        .map(|&g| (0.5 * (g * log).cos() + g * (g * log).sin()) / (0.25 + g * g))375        .sum();376    x - 2.0 * x.sqrt() * waves - (2.0 * PI).ln() - 0.5 * (1.0 - 1.0 / (x * x)).ln()377}378379#[cfg(test)]380mod tests {381    use super::*;382383    fn classic(t: f64) -> f64 {384        0.5 * t * (t / (2.0 * PI)).ln() - 0.5 * t - PI / 8.0385            + 1.0 / (48.0 * t)386            + 7.0 / (5760.0 * t.powi(3))387            + 31.0 / (80640.0 * t.powi(5))388    }389390    #[test]391    fn the_complex_arithmetic_round_trips() {392        let z = Complex::new(-1.5, 2.25);393        let back = z.ln().exp();394        assert!((back - z).abs() < 1e-14);395        assert!((z / z - Complex::new(1.0, 0.0)).abs() < 1e-15);396        assert!((raise(2.0, Complex::new(3.0, 0.0)).re - 8.0).abs() < 1e-13);397        assert!((raise(4.0, Complex::new(0.5, 0.0)) - Complex::new(2.0, 0.0)).abs() < 1e-14);398    }399400    #[test]401    fn theta_meets_the_asymptotic_series_and_the_first_gram_points() {402        let line = Line::new();403        for t in [20.0, 50.0, 100.0, 200.0] {404            assert!((line.theta(t) - classic(t)).abs() < 1e-9, "t {t}");405        }406        assert!(line.theta(0.0).abs() < 1e-12);407        assert!((line.gram(-1) - 9.666_908).abs() < 1e-6);408        assert!((line.gram(0) - 17.845_600).abs() < 1e-6);409        assert!((line.gram(1) - 23.170_283).abs() < 1e-6);410        assert!(line.theta(line.gram(2)).abs() - 2.0 * PI < 1e-10);411    }412413    #[test]414    fn maclaurin_meets_the_known_values_on_the_line() {415        let line = Line::new();416        assert!((line.maclaurin(0.0).re + 1.460_354_508_809_586_8).abs() < 1e-10);417        assert!(line.maclaurin(0.0).im.abs() < 1e-12);418        let one = line.maclaurin(1.0);419        assert!((one.re - 0.143_936_427_077_189).abs() < 1e-9);420        assert!((one.im + 0.722_099_743_531_673).abs() < 1e-9);421        assert!(line.maclaurin(14.134_725).abs() < 1e-5);422        for t in [3.0, 25.0, 140.0] {423            let value = line.maclaurin(t);424            assert!(425                (Complex::turn(line.theta(t)) * value).im.abs() < 1e-9,426                "t {t}"427            );428        }429    }430431    #[test]432    fn the_kernel_pins_its_centre_and_its_removable_points() {433        assert!((kernel(0.5) - (3.0 * PI / 8.0).cos()).abs() < 1e-15);434        assert!((kernel(0.0) - (PI / 8.0).cos()).abs() < 1e-15);435        assert_eq!(kernel(0.25), 0.5);436        assert_eq!(kernel(0.75), 0.5);437        assert!((kernel(0.25 + 1e-9) - 0.5).abs() < 1e-7);438        assert!((kernel(0.75 - 1e-9) - 0.5).abs() < 1e-7);439        let direct = |p: f64| (2.0 * PI * (p * p - p - 1.0 / 16.0)).cos() / (2.0 * PI * p).cos();440        for p in [0.05, 0.1, 0.4, 0.5, 0.6, 0.9, 0.95] {441            assert!((kernel(p) - direct(p)).abs() < 1e-13, "p {p}");442        }443        assert!(corrections(0.5)[1].abs() < 1e-9);444    }445446    #[test]447    fn siegel_meets_maclaurin_beyond_the_join() {448        let line = Line::new();449        assert!(line.seam(JOIN, 60.0, 800) < 5e-5);450        assert!(line.seam(60.0, 250.0, 1900) < 5e-6);451        assert!((line.z(JOIN) - line.exact(JOIN)).abs() < 5e-5);452    }453454    #[test]455    fn the_zeros_and_their_count_are_the_classic_ones() {456        let line = Line::new();457        let first = line.zeros(5);458        let known = [14.134_725, 21.022_040, 25.010_858, 30.424_876, 32.935_062];459        for (got, want) in first.iter().zip(known) {460            assert!((got - want).abs() < 1e-6, "{got} {want}");461        }462        assert_eq!(line.count(100.0), 29);463        assert_eq!(line.count(200.0), 79);464        assert_eq!(line.count(10.0), 0);465        let hundred = line.zeros(100);466        assert_eq!(hundred.len(), 100);467        assert!((hundred[99] - 236.524_230).abs() < 1e-5);468        assert!(hundred.windows(2).all(|w| w[1] > w[0]));469    }470471    #[test]472    fn psi_pins_the_staircase_and_the_smooth_guess() {473        let stair = psi_stair(100);474        assert!((stair[9] - 7.832_0).abs() < 1e-4);475        assert!((stair[99] - 94.045_311).abs() < 1e-5);476        assert_eq!(stair[0], 0.0);477        assert!((stair[7] - 3.0 * 2f64.ln() - 3f64.ln() - 5f64.ln() - 7f64.ln()).abs() < 1e-12);478        assert!((psi_formula(10.0, &[]) - 8.167_1).abs() < 1e-4);479        let line = Line::new();480        let zeros = line.zeros(100);481        let close = psi_formula(100.0, &zeros);482        assert!((close - stair[99]).abs() < 1.0, "{close}");483        assert!((close - stair[99]).abs() < (psi_formula(100.0, &[]) - stair[99]).abs());484    }485}