research-zeta.rs
13.2 kB · rust · 437 lines
1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board, Frame};3use mrlynum::design::elements;4use mrlynum::zeta::{Complex, Line};56const PEEL: u32 = 8;7const DEPTH: usize = 10;8const CUT: f64 = 16.0;9const LMAX: usize = 16;10const NODES: usize = 34;11const HEIGHT: f64 = 40.0;12const WIDE: f64 = 1.20;13const REACH: f64 = 0.70;14const DISC: f64 = 0.45;15const ONE: Complex = Complex::new(1.0, 0.0);16const DATA: &str = "files/figures/data/research-zeta.json";1718struct Ladder {19 q: f64,20 lq: f64,21 k: f64,22 alpha: f64,23 gamma: Vec<f64>,24 qpow: Vec<f64>,25 head: Vec<f64>,26 blogs: Vec<f64>,27 logs: Vec<f64>,28 fall: Vec<f64>,29}3031fn power(log: f64, w: Complex) -> Complex {32 Complex::new(-w.re * log, -w.im * log).exp()33}3435impl Ladder {36 fn new(q: u64, digits: &[u64]) -> Ladder {37 let all = elements(q, digits, DEPTH);38 let split = q.pow(PEEL - 1);39 let stop = q.pow(PEEL);40 let head: Vec<f64> = all41 .iter()42 .filter(|&&n| n < split)43 .map(|&n| n as f64)44 .collect();45 let blogs: Vec<f64> = all46 .iter()47 .filter(|&&n| n >= split && n < stop)48 .map(|&n| (n as f64).ln())49 .collect();50 let logs: Vec<f64> = all51 .iter()52 .filter(|&&n| n >= split)53 .map(|&n| (n as f64).ln())54 .collect();55 let fall: Vec<f64> = all56 .iter()57 .filter(|&&n| n >= split)58 .map(|&n| 1.0 / n as f64)59 .collect();60 let k = digits.len() as f64;61 let mut gamma = vec![k];62 for l in 1..=LMAX {63 gamma.push(digits.iter().map(|&a| (a as f64).powi(l as i32)).sum());64 }65 let qpow = (0..=LMAX).map(|l| (q as f64).powi(-(l as i32))).collect();66 Ladder {67 q: q as f64,68 lq: (q as f64).ln(),69 k,70 alpha: k.ln() / (q as f64).ln(),71 gamma,72 qpow,73 head,74 blogs,75 logs,76 fall,77 }78 }7980 fn cofactor(&self, s: Complex) -> Complex {81 let top = ((CUT - s.re).ceil().max(1.0) as usize).min(NODES);82 let mut g = vec![Complex::default(); NODES + 1];83 let mut pw: Vec<Complex> = self84 .logs85 .iter()86 .map(|&l| power(l, s + top as f64))87 .collect();88 let mut acc = Complex::default();89 for v in &pw {90 acc = acc + *v;91 }92 g[top] = acc;93 for slot in g.iter_mut().take(NODES + 1).skip(top + 1) {94 let mut acc = Complex::default();95 for (p, &drop) in pw.iter_mut().zip(&self.fall) {96 *p = *p * drop;97 acc = acc + *p;98 }99 *slot = acc;100 }101 let mut numerator = Complex::default();102 for node in (0..top).rev() {103 let w = s + node as f64;104 let qw = power(self.lq, w);105 let mut num = Complex::default();106 for &log in &self.blogs {107 num = num + power(log, w);108 }109 let mut binom = ONE;110 for l in 1..=LMAX {111 binom = binom * (-(w + (l as f64 - 1.0))) * (1.0 / l as f64);112 if node + l > NODES {113 break;114 }115 num = num + binom * qw * (self.qpow[l] * self.gamma[l]) * g[node + l];116 }117 g[node] = num / (ONE - qw * self.k);118 if node == 0 {119 numerator = num;120 }121 }122 let mut d = Complex::default();123 for &n in &self.head {124 d = d + power(n.ln(), s);125 }126 (ONE - power(self.lq, s) * self.k) * d + numerator127 }128129 fn centres(&self) -> Vec<(Complex, bool)> {130 let mut out = Vec::new();131 for t in ordinates(self.lq, HEIGHT + 1.0) {132 let zero = Complex::new(self.alpha, t);133 let r0 = self.cofactor(zero) * (1.0 / self.lq);134 out.push((zero, r0.abs() > 1e-12));135 let one = Complex::new(self.alpha - 1.0, t);136 let r1 = one * (self.gamma[1] / (self.k * (self.q - 1.0))) * r0;137 out.push((one, r1.abs() > 1e-12));138 }139 out140 }141}142143fn ordinates(lq: f64, reach: f64) -> Vec<f64> {144 let step = std::f64::consts::TAU / lq;145 (0..)146 .map(|j| j as f64 * step)147 .take_while(|t| *t < reach)148 .collect()149}150151fn refine(l: &Ladder, seed: Complex) -> Option<Complex> {152 let h = 1e-6;153 let mut s = seed;154 for _ in 0..40 {155 let v = l.cofactor(s);156 if v.abs() < 1e-13 {157 return Some(s);158 }159 let a = l.cofactor(s + Complex::new(h, 0.0));160 let b = l.cofactor(s - Complex::new(h, 0.0));161 let d = (a - b) * (1.0 / (2.0 * h));162 if d.abs() < 1e-30 {163 return None;164 }165 let step = v / d;166 if step.abs() > 1.0 {167 return None;168 }169 s = s - step;170 }171 if l.cofactor(s).abs() < 1e-10 {172 Some(s)173 } else {174 None175 }176}177178fn hunt(l: &Ladder, lo: f64, hi: f64) -> Vec<Complex> {179 let cols = ((hi - lo) / 0.03).round() as usize;180 let rows = (HEIGHT / 0.03).round() as usize;181 let at = |i: usize, j: usize| {182 Complex::new(183 lo + (hi - lo) * i as f64 / cols as f64,184 0.02 + (HEIGHT - 0.04) * j as f64 / rows as f64,185 )186 };187 let mut grid = vec![0.0f64; (cols + 1) * (rows + 1)];188 for i in 0..=cols {189 for j in 0..=rows {190 grid[i * (rows + 1) + j] = l.cofactor(at(i, j)).abs();191 }192 }193 let mut found: Vec<Complex> = Vec::new();194 for i in 1..cols {195 for j in 1..rows {196 let here = grid[i * (rows + 1) + j];197 let mut least = true;198 for di in 0..3 {199 for dj in 0..3 {200 if (di, dj) != (1, 1) && grid[(i + di - 1) * (rows + 1) + j + dj - 1] <= here {201 least = false;202 }203 }204 }205 if !least {206 continue;207 }208 if let Some(root) = refine(l, at(i, j)) {209 if root.re < lo || root.re > hi || root.im < 0.02 || root.im > HEIGHT {210 continue;211 }212 if found.iter().all(|z| (*z - root).abs() > 1e-4) {213 found.push(root);214 }215 }216 }217 }218 found.sort_by(|a, b| a.im.partial_cmp(&b.im).unwrap());219 found220}221222fn split(l: &Ladder, zeros: &[Complex]) -> (Vec<Complex>, Vec<Complex>, Vec<Complex>) {223 let centres = l.centres();224 let mut teeth = Vec::new();225 let mut hollow = Vec::new();226 let mut family = Vec::new();227 for z in zeros {228 let mut best = f64::INFINITY;229 let mut live = false;230 for (c, alive) in ¢res {231 let d = (*z - *c).abs();232 if d < best {233 best = d;234 live = *alive;235 }236 }237 if best >= DISC {238 family.push(*z);239 } else if live {240 teeth.push(*z);241 } else {242 hollow.push(*z);243 }244 }245 (teeth, hollow, family)246}247248struct Panel {249 alpha: f64,250 lq: f64,251 teeth: Vec<Complex>,252 hollow: Vec<Complex>,253 family: Vec<Complex>,254}255256fn panel(board: &mut Board, frame: Frame, p: &Panel) {257 let lo = p.alpha - WIDE;258 let hi = p.alpha + REACH;259 let at = |s: Complex| {260 (261 frame.x + frame.w * (s.re - lo) / (hi - lo),262 frame.y + frame.h * (1.0 - s.im / HEIGHT),263 )264 };265 board.rect(frame.x, frame.y, frame.w, frame.h, ink::panel());266 plot::axis(board, frame, ink::line());267 board.segment(268 at(Complex::new(p.alpha, 0.0)),269 at(Complex::new(p.alpha, HEIGHT)),270 1.6,271 ink::fade(ink::dim(), 0.5),272 );273 board.segment(274 at(Complex::new(p.alpha - 1.0, 0.0)),275 at(Complex::new(p.alpha - 1.0, HEIGHT)),276 1.6,277 ink::fade(ink::dim(), 0.3),278 );279 for t in ordinates(p.lq, HEIGHT) {280 if t < 0.4 {281 continue;282 }283 for line in [p.alpha, p.alpha - 1.0] {284 let (x, y) = at(Complex::new(line, t));285 board.ring(x, y, 11.0, 1.6, ink::fade(ink::dim(), 0.9));286 }287 }288 for z in &p.family {289 let (x, y) = at(*z);290 board.disc(x, y, 7.0, ink::blue());291 }292 for z in &p.teeth {293 let (x, y) = at(*z);294 board.disc(x, y, 7.0, ink::yellow());295 }296 for z in &p.hollow {297 let (x, y) = at(*z);298 board.disc(x, y, 7.0, ink::fade(ink::fg(), 0.9));299 }300}301302fn flat(zs: &[Complex]) -> String {303 zs.iter()304 .flat_map(|z| [z.re, z.im])305 .map(|v| format!("{v:?}"))306 .collect::<Vec<String>>()307 .join(",")308}309310fn field(text: &str, key: &str) -> Vec<f64> {311 let head = format!("\"{key}\": [");312 let from = text.find(&head).expect("the data file lost a key") + head.len();313 let upto = from314 + text[from..]315 .find(']')316 .expect("the data file lost a bracket");317 text[from..upto]318 .split(',')319 .map(str::trim)320 .filter(|t| !t.is_empty())321 .map(|t| t.parse::<f64>().expect("the data file lost a number"))322 .collect()323}324325fn points(v: &[f64]) -> Vec<Complex> {326 v.chunks(2).map(|c| Complex::new(c[0], c[1])).collect()327}328329fn read(text: &str, name: &str) -> Panel {330 let axis = field(text, &format!("{name}_axis"));331 let teeth = points(&field(text, &format!("{name}_teeth")));332 let hollow = points(&field(text, &format!("{name}_hollow")));333 let family = points(&field(text, &format!("{name}_family")));334 assert!(axis.len() == 2 && axis[1] > 0.0 && !family.is_empty());335 Panel {336 alpha: axis[0],337 lq: axis[1],338 teeth,339 hollow,340 family,341 }342}343344fn compute() -> Result<()> {345 let design = Ladder::new(3, &[0, 1]);346 let full = Ladder::new(2, &[0, 1]);347 let known = Line::new().zeros(6);348 assert!(full.cofactor(Complex::new(0.5, known[0])).abs() < 1e-8);349350 let dz = hunt(&design, design.alpha - 0.92, design.alpha + 3.02);351 let fz = hunt(&full, full.alpha - 0.92, full.alpha + 3.02);352 let (fteeth, fhollow, ffamily) = split(&full, &fz);353 let (dteeth, dhollow, dfamily) = split(&design, &dz);354355 assert_eq!(fz.len(), 10);356 assert_eq!(fteeth.len(), 0);357 assert_eq!(fhollow.len(), 4);358 assert_eq!(ffamily.len(), 6);359 for z in &fhollow {360 assert!((z.re - full.alpha).abs() < 1e-6);361 }362 for (z, g) in ffamily.iter().zip(&known) {363 assert!((z.re - 0.5).abs() < 1e-6 && (z.im - g).abs() < 1e-6);364 }365366 assert_eq!(dz.len(), 13);367 assert_eq!(dhollow.len(), 0);368 assert_eq!(dteeth.len() + dfamily.len(), 13);369 assert!(dteeth.iter().all(|z| (z.re - design.alpha).abs() > 1e-3));370 assert!(dfamily.iter().any(|z| (z.re - 0.391038600).abs() < 1e-6));371 assert!(dfamily372 .iter()373 .any(|z| (z.re + 0.273079611).abs() < 1e-6 && (z.im - 39.262315320).abs() < 1e-6));374 let least = dfamily.iter().map(|z| z.re).fold(f64::INFINITY, f64::min);375 let most = dfamily376 .iter()377 .map(|z| z.re)378 .fold(f64::NEG_INFINITY, f64::max);379 assert!((least + 0.273079611).abs() < 1e-6 && (most - 0.391038600).abs() < 1e-6);380 assert!(dz381 .iter()382 .all(|z| z.re > design.alpha - WIDE && z.re < design.alpha + REACH));383 assert!(fz384 .iter()385 .all(|z| z.re > full.alpha - WIDE && z.re < full.alpha + REACH));386387 let daxis = format!("{:?},{:?}", design.alpha, design.lq);388 let faxis = format!("{:?},{:?}", full.alpha, full.lq);389 let text = format!(390 "{{\n \"design_axis\": [{}],\n \"design_teeth\": [{}],\n \"design_hollow\": [{}],\n \"design_family\": [{}],\n \"full_axis\": [{}],\n \"full_teeth\": [{}],\n \"full_hollow\": [{}],\n \"full_family\": [{}]\n}}\n",391 daxis,392 flat(&dteeth),393 flat(&dhollow),394 flat(&dfamily),395 faxis,396 flat(&fteeth),397 flat(&fhollow),398 flat(&ffamily)399 );400 let path = mrlyfig::out::root().join(DATA);401 let folder = path.parent().expect("the data path lost its folder");402 std::fs::create_dir_all(folder)403 .map_err(|e| mrlycore::MrlyError::Value(format!("cannot make {folder:?}: {e}")))?;404 std::fs::write(&path, text)405 .map_err(|e| mrlycore::MrlyError::Value(format!("cannot write {path:?}: {e}")))?;406 println!("data research-zeta {} zeros", dz.len() + fz.len());407 Ok(())408}409410fn render() -> Result<()> {411 let path = mrlyfig::out::root().join(DATA);412 let text = std::fs::read_to_string(&path).map_err(|e| {413 mrlycore::MrlyError::Value(format!(414 "cannot read {path:?}: {e}; run the example with compute"415 ))416 })?;417 let design = read(&text, "design");418 let full = read(&text, "full");419 let mut board = Board::square();420 let area = board.area(0.08);421 let wide = area.w * 0.45;422 let tall = area.h * 0.80;423 let top = Frame::new(area.x, area.y, wide, tall);424 let bottom = Frame::new(area.x + area.w - wide, area.y + area.h - tall, wide, tall);425 panel(&mut board, top, &design);426 panel(&mut board, bottom, &full);427 save("research-zeta", &board)?;428 Ok(())429}430431fn main() -> Result<()> {432 if std::env::args().nth(1).as_deref() == Some("compute") {433 compute()434 } else {435 render()436 }437}