paper-novelty-meter.rs
3.9 kB · rust · 119 lines
1use mrlycore::errors::Result;2use mrlycore::Color;3use mrlyfig::{ink, plot, save, Board, Frame};4use mrlynum::lattice::totients;5use mrlynum::zeta::{novelty_main, novelty_wave, sharp_novelty, smoothed_novelty, Line};6use std::f64::consts::PI;78const REACH: usize = 3_000_000;9const LOW: f64 = 8.0;10const HIGH: f64 = 20.5;11const PER_OCTAVE: usize = 16;12const ZEROS: usize = 29;13const CURVE: usize = 2400;1415fn beads(board: &mut Board, area: Frame, values: &[f64], span: f64, r: f64, color: Color) {16 let last = (values.len() - 1) as f64;17 let pts: Vec<(f64, f64)> = values18 .iter()19 .enumerate()20 .map(|(i, v)| {21 (22 area.x + area.w * i as f64 / last,23 area.y + area.h * (0.5 - 0.5 * v / span),24 )25 })26 .collect();27 plot::dots(board, &pts, r, color);28}2930fn reach(values: &[f64]) -> f64 {31 values.iter().fold(0.0f64, |a, v| a.max(v.abs())) * 1.1232}3334fn main() -> Result<()> {35 let phi = totients(REACH);36 assert_eq!(&phi[1..=10], &[1, 1, 2, 2, 4, 2, 6, 4, 6, 4]);37 let mut prefix = vec![0u64; REACH + 1];38 for n in 1..=REACH {39 prefix[n] = prefix[n - 1] + phi[n];40 }41 let total = prefix[REACH] as f64;42 let expect = 3.0 * (REACH as f64).powi(2) / (PI * PI);43 assert!((total - expect).abs() < REACH as f64 * (REACH as f64).ln());4445 let main_bump = novelty_main();46 assert!((main_bump - 6.0 / (PI * PI) * 0.575_725_895_994).abs() < 1e-9);4748 let samples = ((HIGH - LOW) * PER_OCTAVE as f64).round() as usize + 1;49 assert_eq!(samples, 201);50 let js: Vec<f64> = (0..samples)51 .map(|k| LOW + k as f64 / PER_OCTAVE as f64)52 .collect();53 let smooth: Vec<f64> = js54 .iter()55 .map(|&j| {56 let y = 2.0f64.powf(-j);57 smoothed_novelty(&phi, y, main_bump) / y.powf(1.5)58 })59 .collect();60 let rough: Vec<f64> = js61 .iter()62 .map(|&j| {63 let y = 2.0f64.powf(-j);64 sharp_novelty(&prefix, y) / y65 })66 .collect();6768 let line = Line::new();69 let gammas = line.zeros(ZEROS);70 assert!((gammas[0] - 14.134_725).abs() < 1e-6);71 assert!((gammas[ZEROS - 1] - 98.831_194).abs() < 1e-5);72 let coef = line.novelty_coefficients(&gammas);73 assert!((coef[0].abs() - 0.1879).abs() < 5e-4);74 assert!((coef[9].abs() - 4.286e-3).abs() < 5e-6);75 let wave = |j: f64| novelty_wave(&gammas, &coef, -j * 2.0f64.ln());76 let peak = smooth.iter().fold(0.0f64, |a, v| a.max(v.abs()));77 let miss = js78 .iter()79 .zip(&smooth)80 .map(|(&j, &v)| (v - wave(j)).abs())81 .fold(0.0f64, f64::max);82 assert!(miss / peak < 1e-2);8384 let mut board = Board::square();85 let frame = board.frame(0.08);86 let gap = 28.0;87 let top_h = frame.h * 0.62 - gap / 2.0;88 let top = Frame::new(frame.x, frame.y, frame.w, top_h);89 let foot = Frame::new(90 frame.x,91 frame.y + top_h + gap,92 frame.w,93 frame.h - top_h - gap,94 );95 for area in [top, foot] {96 plot::axis(&mut board, area, ink::line());97 let mid = area.y + area.h / 2.0;98 board.segment((area.x, mid), (area.x + area.w, mid), 1.0, ink::line());99 }100 let inner_top = top.inset(16.0);101 let inner_foot = foot.inset(16.0);102 let span_top = reach(&smooth);103 let span_foot = reach(&rough);104 let curve: Vec<(f64, f64)> = (0..=CURVE)105 .map(|k| {106 let t = k as f64 / CURVE as f64;107 let v = wave(LOW + (HIGH - LOW) * t);108 (109 inner_top.x + inner_top.w * t,110 inner_top.y + inner_top.h * (0.5 - 0.5 * v / span_top),111 )112 })113 .collect();114 board.polyline(&curve, 2.0, ink::orange());115 beads(&mut board, inner_top, &smooth, span_top, 3.4, ink::blue());116 beads(&mut board, inner_foot, &rough, span_foot, 3.4, ink::blue());117 save("paper-novelty-meter", &board)?;118 Ok(())119}