zeta.rs
22.1 kB · rust · 630 lines
1use crate::num::classics::primes;2use crate::num::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;13const NODES: usize = 4096;1415/// A complex number: a real and an imaginary part.16#[derive(Clone, Copy, Debug, Default, PartialEq)]17pub struct Complex {18 /// The real part.19 pub re: f64,20 /// The imaginary part.21 pub im: f64,22}2324impl Complex {25 /// Builds a complex number from its parts.26 pub const fn new(re: f64, im: f64) -> Complex {27 Complex { re, im }28 }29 /// Returns the modulus.30 pub fn abs(self) -> f64 {31 self.re.hypot(self.im)32 }33 /// Returns the principal argument.34 pub fn arg(self) -> f64 {35 self.im.atan2(self.re)36 }37 /// Returns the exponential.38 pub fn exp(self) -> Complex {39 let r = self.re.exp();40 Complex::new(r * self.im.cos(), r * self.im.sin())41 }42 /// Returns the principal logarithm.43 pub fn ln(self) -> Complex {44 Complex::new(self.abs().ln(), self.arg())45 }46 /// Returns a unit complex number at the given angle.47 pub fn turn(angle: f64) -> Complex {48 Complex::new(angle.cos(), angle.sin())49 }50}5152impl Add for Complex {53 type Output = Complex;54 fn add(self, other: Complex) -> Complex {55 Complex::new(self.re + other.re, self.im + other.im)56 }57}5859impl Add<f64> for Complex {60 type Output = Complex;61 fn add(self, other: f64) -> Complex {62 Complex::new(self.re + other, self.im)63 }64}6566impl Sub for Complex {67 type Output = Complex;68 fn sub(self, other: Complex) -> Complex {69 Complex::new(self.re - other.re, self.im - other.im)70 }71}7273impl Sub<f64> for Complex {74 type Output = Complex;75 fn sub(self, other: f64) -> Complex {76 Complex::new(self.re - other, self.im)77 }78}7980impl Mul for Complex {81 type Output = Complex;82 fn mul(self, other: Complex) -> Complex {83 Complex::new(84 self.re * other.re - self.im * other.im,85 self.re * other.im + self.im * other.re,86 )87 }88}8990impl Mul<f64> for Complex {91 type Output = Complex;92 fn mul(self, other: f64) -> Complex {93 Complex::new(self.re * other, self.im * other)94 }95}9697impl Div for Complex {98 type Output = Complex;99 fn div(self, other: Complex) -> Complex {100 let d = other.re * other.re + other.im * other.im;101 Complex::new(102 (self.re * other.re + self.im * other.im) / d,103 (self.im * other.re - self.re * other.im) / d,104 )105 }106}107108impl Neg for Complex {109 type Output = Complex;110 fn neg(self) -> Complex {111 Complex::new(-self.re, -self.im)112 }113}114115/// Returns a positive real base raised to a complex exponent.116pub fn raise(base: f64, exponent: Complex) -> Complex {117 (exponent * base.ln()).exp()118}119120/// Returns the Riemann-Siegel kernel, the cosine ratio that leads the remainder, in the form that stays finite at its removable points.121pub fn kernel(p: f64) -> f64 {122 let d = if p < 0.5 { 0.25 - p } else { p - 0.75 };123 let below = (2.0 * PI * d).sin();124 if below == 0.0 {125 return 0.5;126 }127 (PI * d * (1.0 + 2.0 * d)).sin() / below128}129130fn difference(order: usize, p: f64, step: f64) -> f64 {131 let mut sum = 0.0;132 let mut binomial = 1.0;133 for k in 0..=order {134 let offset = (order as f64 / 2.0 - k as f64) * step;135 let sign = if k.is_multiple_of(2) { 1.0 } else { -1.0 };136 sum += sign * binomial * kernel(p + offset);137 binomial = binomial * (order - k) as f64 / (k + 1) as f64;138 }139 sum / step.powi(order as i32)140}141142fn derivative(order: usize, p: f64) -> f64 {143 (4.0 * difference(order, p, STENCIL) - difference(order, p, 2.0 * STENCIL)) / 3.0144}145146/// Returns the first four Riemann-Siegel corrections at the fractional part p: the kernel and its derivatives by central differences with one Richardson step.147pub fn corrections(p: f64) -> [f64; 4] {148 let (pi2, pi4, pi6) = (PI * PI, PI.powi(4), PI.powi(6));149 [150 kernel(p),151 -derivative(3, p) / (96.0 * pi2),152 derivative(6, p) / (18_432.0 * pi4) + derivative(2, p) / (64.0 * pi2),153 -derivative(9, p) / (5_308_416.0 * pi6)154 - derivative(5, p) / (3_840.0 * pi4)155 - derivative(1, p) / (64.0 * pi2),156 ]157}158159/// The critical line: the Bernoulli numbers and the Euler-Maclaurin weights the two engines share, built once.160pub struct Line {161 bern: Vec<f64>,162 tail: Vec<f64>,163}164165impl Default for Line {166 fn default() -> Line {167 Line::new()168 }169}170171impl Line {172 /// Builds the line: the even Bernoulli numbers through the fourteenth and their Euler-Maclaurin weights.173 pub fn new() -> Line {174 let fractions = bernoulli(2 * TAIL + 1);175 let bern: Vec<f64> = (0..=TAIL)176 .map(|k| {177 let (num, den) = fractions[2 * k];178 num as f64 / den as f64179 })180 .collect();181 let mut factorial = 1.0;182 let tail = (0..=TAIL)183 .map(|k| {184 if k > 0 {185 factorial *= ((2 * k - 1) * 2 * k) as f64;186 }187 bern[k] / factorial188 })189 .collect();190 Line { bern, tail }191 }192 /// 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.193 pub fn theta(&self, t: f64) -> f64 {194 let z = Complex::new(0.25, 0.5 * t);195 let w = z + SHIFT as f64;196 let mut arg = ((w - 0.5) * w.ln() - w).im;197 for k in 0..SHIFT {198 arg -= (z + k as f64).arg();199 }200 let inverse = Complex::new(1.0, 0.0) / w;201 let square = inverse * inverse;202 let mut power = inverse;203 for k in 1..=TAIL {204 arg += self.bern[k] / ((2 * k * (2 * k - 1)) as f64) * power.im;205 power = power * square;206 }207 arg - 0.5 * t * PI.ln()208 }209 /// Returns zeta at one half plus i t by the complex Euler-Maclaurin sum: t plus ten terms and seven Bernoulli corrections.210 pub fn maclaurin(&self, t: f64) -> Complex {211 let s = Complex::new(0.5, t);212 let count = t.abs() as usize + SHIFT;213 let mut sum = Complex::default();214 for k in 1..=count {215 sum = sum + raise(k as f64, -s);216 }217 let base = count as f64;218 let mut out = sum + raise(base, -s + 1.0) / (s - 1.0) - raise(base, -s) * 0.5;219 let mut rising = s;220 let mut power = raise(base, -s - 1.0);221 let square = 1.0 / (base * base);222 for k in 1..=TAIL {223 out = out + rising * power * self.tail[k];224 rising = rising * (s + (2 * k - 1) as f64) * (s + (2 * k) as f64);225 power = power * square;226 }227 out228 }229 /// Returns zeta and its derivative together at any complex s but one, by the same Euler-Maclaurin sum: the modulus of t plus ten terms and seven Bernoulli corrections, each term differentiated in s.230 pub fn pair(&self, s: Complex) -> (Complex, Complex) {231 let one = Complex::new(1.0, 0.0);232 let count = s.im.abs() as usize + SHIFT;233 let mut value = Complex::default();234 let mut slope = Complex::default();235 for k in 1..=count {236 let kf = k as f64;237 let term = raise(kf, -s);238 value = value + term;239 slope = slope - term * kf.ln();240 }241 let base = count as f64;242 let log = base.ln();243 let head = raise(base, -s + 1.0) / (s - 1.0);244 let half = raise(base, -s) * 0.5;245 value = value + head - half;246 slope = slope - head * log - head / (s - 1.0) + half * log;247 let mut rising = s;248 let mut ratio = one / s;249 let mut power = raise(base, -s - 1.0);250 let square = 1.0 / (base * base);251 for (k, &weight) in self.tail.iter().enumerate().skip(1) {252 let term = rising * power * weight;253 value = value + term;254 slope = slope + term * (ratio - log);255 let (a, b) = (s + (2 * k - 1) as f64, s + (2 * k) as f64);256 rising = rising * a * b;257 ratio = ratio + one / a + one / b;258 power = power * square;259 }260 (value, slope)261 }262 /// Returns the wave coefficient of every zero at the given ordinates: F(rho) zeta(rho - 1) over zeta'(rho) at rho one half plus i gamma, F the Mellin transform of the bump.263 pub fn novelty_coefficients(&self, gammas: &[f64]) -> Vec<Complex> {264 gammas265 .iter()266 .map(|&g| {267 let rho = Complex::new(0.5, g);268 let (left, _) = self.pair(rho - 1.0);269 let (_, prime) = self.pair(rho);270 mellin(rho) * left / prime271 })272 .collect()273 }274 /// Returns Z(t) from the Euler-Maclaurin value turned onto the real axis.275 pub fn exact(&self, t: f64) -> f64 {276 (Complex::turn(self.theta(t)) * self.maclaurin(t)).re277 }278 /// Returns Z(t) by the Riemann-Siegel formula: the main sum and the first four corrections.279 pub fn siegel(&self, t: f64) -> f64 {280 let a = (t / (2.0 * PI)).sqrt();281 let whole = a.floor();282 let theta = self.theta(t);283 let mut sum = 0.0;284 for k in 1..=whole as usize {285 let kf = k as f64;286 sum += (theta - t * kf.ln()).cos() / kf.sqrt();287 }288 let sign = if (whole as u64).is_multiple_of(2) {289 -1.0290 } else {291 1.0292 };293 let mut weight = 1.0 / a.sqrt();294 let mut rest = 0.0;295 for c in corrections(a - whole) {296 rest += c * weight;297 weight /= a;298 }299 2.0 * sum + sign * rest300 }301 /// Returns Z(t): Euler-Maclaurin below the join, Riemann-Siegel above.302 pub fn z(&self, t: f64) -> f64 {303 if t < JOIN {304 self.exact(t)305 } else {306 self.siegel(t)307 }308 }309 /// Returns zeta on the line and Z(t) together, from the engine that serves the t.310 pub fn point(&self, t: f64) -> (Complex, f64) {311 if t < JOIN {312 let value = self.maclaurin(t);313 (value, (Complex::turn(self.theta(t)) * value).re)314 } else {315 let z = self.siegel(t);316 (Complex::turn(-self.theta(t)) * z, z)317 }318 }319 /// Returns the largest gap between the two engines over the t range on a grid.320 pub fn seam(&self, t0: f64, t1: f64, steps: usize) -> f64 {321 (0..=steps)322 .map(|k| {323 let t = t0 + (t1 - t0) * k as f64 / steps as f64;324 (self.siegel(t) - self.exact(t)).abs()325 })326 .fold(0.0, f64::max)327 }328 /// Returns the n-th Gram point, where theta is n pi, by Newton from the right.329 pub fn gram(&self, n: i64) -> f64 {330 let target = n as f64 * PI;331 let mut t = 2.0 * PI * (n as f64 + 2.0).max(3.0);332 for _ in 0..100 {333 let step = (self.theta(t) - target) / (0.5 * (t / (2.0 * PI)).ln());334 t -= step;335 if step.abs() < 1e-12 {336 break;337 }338 }339 t340 }341 fn brackets(&self, limit: f64, count: usize, exact: bool) -> Vec<(f64, f64)> {342 let z = |t: f64| if exact { self.exact(t) } else { self.z(t) };343 let mut out = Vec::new();344 let mut n = -1;345 let mut left = self.gram(n);346 let mut previous = left;347 let mut before = z(left);348 'walk: while out.len() < count && left < limit {349 let right = self.gram(n + 1);350 for k in 1..=STEPS {351 let t = (left + (right - left) * k as f64 / STEPS as f64).min(limit);352 let now = z(t);353 if before * now < 0.0 {354 out.push((previous, t));355 }356 previous = t;357 before = now;358 if t >= limit {359 break 'walk;360 }361 }362 left = right;363 n += 1;364 }365 out.truncate(count);366 out367 }368 fn bisect(&self, (mut a, mut b): (f64, f64)) -> f64 {369 let mut fa = self.exact(a);370 while b - a > TOLERANCE {371 let mid = 0.5 * (a + b);372 let fm = self.exact(mid);373 if fa * fm <= 0.0 {374 b = mid;375 } else {376 a = mid;377 fa = fm;378 }379 }380 0.5 * (a + b)381 }382 /// Returns the first zeros on the line: sign changes of Z between Gram points, refined by bisection on Euler-Maclaurin to a billionth.383 pub fn zeros(&self, count: usize) -> Vec<f64> {384 self.brackets(f64::INFINITY, count, true)385 .into_iter()386 .map(|pair| self.bisect(pair))387 .collect()388 }389 /// Counts the zeros on the line below t.390 pub fn count(&self, t: f64) -> usize {391 self.brackets(t, usize::MAX, false).len()392 }393}394395/// Returns the Chebyshev staircase at every whole number from one to x: the sum of ln p over the prime powers up to each.396pub fn psi_stair(x: usize) -> Vec<f64> {397 let mut jumps = vec![0.0; x + 1];398 for p in primes(x) {399 let mut q = p;400 while q <= x {401 jumps[q] += (p as f64).ln();402 q *= p;403 }404 }405 let mut sum = 0.0;406 jumps[1..]407 .iter()408 .map(|jump| {409 sum += jump;410 sum411 })412 .collect()413}414415/// 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.416pub fn psi_formula(x: f64, gammas: &[f64]) -> f64 {417 let log = x.ln();418 let waves: f64 = gammas419 .iter()420 .map(|&g| (0.5 * (g * log).cos() + g * (g * log).sin()) / (0.25 + g * g))421 .sum();422 x - 2.0 * x.sqrt() * waves - (2.0 * PI).ln() - 0.5 * (1.0 - 1.0 / (x * x)).ln()423}424425// NOVELTY426427/// The smooth window on [1, 2]: exp(4 - 1/((u - 1)(2 - u))) inside, zero outside, every derivative vanishing at the ends and a peak of one at u = 3/2.428pub fn bump(u: f64) -> f64 {429 if u <= 1.0 || u >= 2.0 {430 return 0.0;431 }432 (4.0 - 1.0 / ((u - 1.0) * (2.0 - u))).exp()433}434435/// Returns the Mellin transform of the bump at a complex s, the integral of bump(u) u^(s - 1) over [1, 2], by a 4096-node midpoint rule.436pub fn mellin(s: Complex) -> Complex {437 let step = 1.0 / NODES as f64;438 let mut sum = Complex::default();439 for k in 1..NODES {440 let u = 1.0 + k as f64 * step;441 sum = sum + raise(u, s - 1.0) * (bump(u) * step);442 }443 sum444}445446/// Returns the main term of the smoothed novelty: six over pi squared times the bump's transform at two.447pub fn novelty_main() -> f64 {448 6.0 / (PI * PI) * mellin(Complex::new(2.0, 0.0)).re449}450451/// Returns the smoothed novelty error at y: y squared times the totients weighed by the bump at n y, less the main term given; the totients must reach 2 over y.452pub fn smoothed_novelty(phi: &[u64], y: f64, main: f64) -> f64 {453 let lo = (1.0 / y).ceil() as usize;454 let hi = (2.0 / y).floor() as usize;455 let sum: f64 = phi[lo..=hi]456 .iter()457 .zip(lo..)458 .map(|(&p, n)| p as f64 * bump(n as f64 * y))459 .sum();460 y * y * sum - main461}462463/// Returns the sharp novelty error at y: y squared times the totient sum over the scales from 1 over y to 2 over y, both ends in, less nine over pi squared, from the prefix sums of the totients, which must reach 2 over y.464pub fn sharp_novelty(prefix: &[u64], y: f64) -> f64 {465 let lo = (1.0 / y).ceil() as usize;466 let hi = (2.0 / y).floor() as usize;467 y * y * (prefix[hi] - prefix[lo - 1]) as f64 - 9.0 / (PI * PI)468}469470/// Sums the waves of the zeros at log y: twice the real part of the coefficients times y to the minus i gamma, the smoothed error over y to the three halves that the zeros predict.471pub fn novelty_wave(gammas: &[f64], coef: &[Complex], log_y: f64) -> f64 {472 let sum: f64 = gammas473 .iter()474 .zip(coef)475 .map(|(g, c)| (*c * Complex::turn(-g * log_y)).re)476 .sum();477 2.0 * sum478}479480#[cfg(test)]481mod tests {482 use super::*;483 use crate::num::lattice::totients;484485 fn classic(t: f64) -> f64 {486 0.5 * t * (t / (2.0 * PI)).ln() - 0.5 * t - PI / 8.0487 + 1.0 / (48.0 * t)488 + 7.0 / (5760.0 * t.powi(3))489 + 31.0 / (80640.0 * t.powi(5))490 }491492 #[test]493 fn the_complex_arithmetic_round_trips() {494 let z = Complex::new(-1.5, 2.25);495 let back = z.ln().exp();496 assert!((back - z).abs() < 1e-14);497 assert!((z / z - Complex::new(1.0, 0.0)).abs() < 1e-15);498 assert!((raise(2.0, Complex::new(3.0, 0.0)).re - 8.0).abs() < 1e-13);499 assert!((raise(4.0, Complex::new(0.5, 0.0)) - Complex::new(2.0, 0.0)).abs() < 1e-14);500 }501502 #[test]503 fn theta_meets_the_asymptotic_series_and_the_first_gram_points() {504 let line = Line::new();505 for t in [20.0, 50.0, 100.0, 200.0] {506 assert!((line.theta(t) - classic(t)).abs() < 1e-9, "t {t}");507 }508 assert!(line.theta(0.0).abs() < 1e-12);509 assert!((line.gram(-1) - 9.666_908).abs() < 1e-6);510 assert!((line.gram(0) - 17.845_600).abs() < 1e-6);511 assert!((line.gram(1) - 23.170_283).abs() < 1e-6);512 assert!(line.theta(line.gram(2)).abs() - 2.0 * PI < 1e-10);513 }514515 #[test]516 fn maclaurin_meets_the_known_values_on_the_line() {517 let line = Line::new();518 assert!((line.maclaurin(0.0).re + 1.460_354_508_809_586_8).abs() < 1e-10);519 assert!(line.maclaurin(0.0).im.abs() < 1e-12);520 let one = line.maclaurin(1.0);521 assert!((one.re - 0.143_936_427_077_189).abs() < 1e-9);522 assert!((one.im + 0.722_099_743_531_673).abs() < 1e-9);523 assert!(line.maclaurin(14.134_725).abs() < 1e-5);524 for t in [3.0, 25.0, 140.0] {525 let value = line.maclaurin(t);526 assert!(527 (Complex::turn(line.theta(t)) * value).im.abs() < 1e-9,528 "t {t}"529 );530 }531 }532533 #[test]534 fn the_kernel_pins_its_centre_and_its_removable_points() {535 assert!((kernel(0.5) - (3.0 * PI / 8.0).cos()).abs() < 1e-15);536 assert!((kernel(0.0) - (PI / 8.0).cos()).abs() < 1e-15);537 assert_eq!(kernel(0.25), 0.5);538 assert_eq!(kernel(0.75), 0.5);539 assert!((kernel(0.25 + 1e-9) - 0.5).abs() < 1e-7);540 assert!((kernel(0.75 - 1e-9) - 0.5).abs() < 1e-7);541 let direct = |p: f64| (2.0 * PI * (p * p - p - 1.0 / 16.0)).cos() / (2.0 * PI * p).cos();542 for p in [0.05, 0.1, 0.4, 0.5, 0.6, 0.9, 0.95] {543 assert!((kernel(p) - direct(p)).abs() < 1e-13, "p {p}");544 }545 assert!(corrections(0.5)[1].abs() < 1e-9);546 }547548 #[test]549 fn siegel_meets_maclaurin_beyond_the_join() {550 let line = Line::new();551 assert!(line.seam(JOIN, 60.0, 800) < 5e-5);552 assert!(line.seam(60.0, 250.0, 1900) < 5e-6);553 assert!((line.z(JOIN) - line.exact(JOIN)).abs() < 5e-5);554 }555556 #[test]557 fn the_zeros_and_their_count_are_the_classic_ones() {558 let line = Line::new();559 let first = line.zeros(5);560 let known = [14.134_725, 21.022_040, 25.010_858, 30.424_876, 32.935_062];561 for (got, want) in first.iter().zip(known) {562 assert!((got - want).abs() < 1e-6, "{got} {want}");563 }564 assert_eq!(line.count(100.0), 29);565 assert_eq!(line.count(200.0), 79);566 assert_eq!(line.count(10.0), 0);567 let hundred = line.zeros(100);568 assert_eq!(hundred.len(), 100);569 assert!((hundred[99] - 236.524_230).abs() < 1e-5);570 assert!(hundred.windows(2).all(|w| w[1] > w[0]));571 }572573 #[test]574 fn the_pair_pins_zeta_off_the_line_and_its_slope() {575 let line = Line::new();576 let (below, _) = line.pair(Complex::new(-0.5, 0.0));577 assert!((below.re + 0.207_886_224_977_354_6).abs() < 1e-10);578 let (_, slope) = line.pair(Complex::new(0.5, 0.0));579 assert!((slope.re + 3.922_646_139_209_15).abs() < 1e-9);580 let (on, _) = line.pair(Complex::new(0.5, 30.0));581 assert!((on - line.maclaurin(30.0)).abs() < 1e-9);582 let (root, prime) = line.pair(Complex::new(0.5, 14.134_725_141_734_7));583 assert!(root.abs() < 1e-8);584 assert!((prime.abs() - 0.793_16).abs() < 1e-4);585 }586587 #[test]588 fn the_novelty_error_is_the_wave_of_the_first_zeros() {589 assert!((novelty_main() - 6.0 / (PI * PI) * 0.575_725_895_994).abs() < 1e-9);590 assert_eq!(bump(1.5), 1.0);591 assert_eq!(bump(1.0), 0.0);592 let line = Line::new();593 let gammas = line.zeros(10);594 let coef = line.novelty_coefficients(&gammas);595 assert!((coef[0].abs() - 0.1879).abs() < 5e-4);596 assert!((coef[9].abs() - 4.286e-3).abs() < 5e-6);597 let phi = totients(1 << 15);598 let mut prefix = vec![0u64; phi.len()];599 for n in 1..phi.len() {600 prefix[n] = prefix[n - 1] + phi[n];601 }602 let main = novelty_main();603 let (mut peak, mut miss) = (0.0f64, 0.0f64);604 for k in 0..=96 {605 let j = 8.0 + k as f64 / 16.0;606 let y = 2.0f64.powf(-j);607 let dot = smoothed_novelty(&phi, y, main) / y.powf(1.5);608 peak = peak.max(dot.abs());609 miss = miss.max((dot - novelty_wave(&gammas, &coef, y.ln())).abs());610 }611 assert!(miss / peak < 5e-2, "{miss} {peak}");612 let rough = sharp_novelty(&prefix, 2.0f64.powf(-10.0)) / 2.0f64.powf(-10.0);613 assert!(rough.abs() < 2.0 && rough != 0.0);614 }615616 #[test]617 fn psi_pins_the_staircase_and_the_smooth_guess() {618 let stair = psi_stair(100);619 assert!((stair[9] - 7.832_0).abs() < 1e-4);620 assert!((stair[99] - 94.045_311).abs() < 1e-5);621 assert_eq!(stair[0], 0.0);622 assert!((stair[7] - 3.0 * 2f64.ln() - 3f64.ln() - 5f64.ln() - 7f64.ln()).abs() < 1e-12);623 assert!((psi_formula(10.0, &[]) - 8.167_1).abs() < 1e-4);624 let line = Line::new();625 let zeros = line.zeros(100);626 let close = psi_formula(100.0, &zeros);627 assert!((close - stair[99]).abs() < 1.0, "{close}");628 assert!((close - stair[99]).abs() < (psi_formula(100.0, &[]) - stair[99]).abs());629 }630}