novelty.rs

1use crate::Fault;2use mrlynum::lattice::totients;3use mrlynum::zeta::{novelty_main, novelty_wave, sharp_novelty, smoothed_novelty, Complex, Line};4use wasm_bindgen::prelude::*;56const LOW: f64 = 8.0;7const HIGHEST: f64 = 21.0;8const PER_OCTAVE: usize = 32;9const ZEROS: usize = 138;1011/// The novelty meter: the totients sieved once, the smoothed and the sharp novelty error read on a log grid of y, and the first zeros with the coefficients of their waves.12#[wasm_bindgen]13pub struct Novelty {14    heights: Vec<f64>,15    smooth: Vec<f64>,16    rough: Vec<f64>,17    gammas: Vec<f64>,18    coef: Vec<Complex>,19    sieve: usize,20}2122#[wasm_bindgen]23impl Novelty {24    /// Builds the meter on y = 2^-j for j from eight to the height at the given samples per octave, the totients sieved to two to the height plus one, with the given count of zeros, at most 138.25    #[wasm_bindgen(constructor)]26    pub fn new(high: f64, per_octave: usize, zeros: usize) -> Result<Novelty, Fault> {27        if !(LOW + 1.0..=HIGHEST).contains(&high) {28            return Err(Fault::new(format!(29                "the height runs from {} to {HIGHEST}.",30                LOW + 1.031            )));32        }33        if !(1..=PER_OCTAVE).contains(&per_octave) {34            return Err(Fault::new(format!(35                "an octave takes 1 to {PER_OCTAVE} samples."36            )));37        }38        if zeros > ZEROS {39            return Err(Fault::new(format!("the meter holds {ZEROS} zeros.")));40        }41        let sieve = 2.0f64.powf(high + 1.0).ceil() as usize;42        let phi = totients(sieve);43        let mut prefix = vec![0u64; sieve + 1];44        for n in 1..=sieve {45            prefix[n] = prefix[n - 1] + phi[n];46        }47        let main = novelty_main();48        let samples = ((high - LOW) * per_octave as f64).round() as usize + 1;49        let heights: Vec<f64> = (0..samples)50            .map(|k| LOW + k as f64 / per_octave as f64)51            .collect();52        let smooth = heights53            .iter()54            .map(|&j| {55                let y = 2.0f64.powf(-j);56                smoothed_novelty(&phi, y, main) / y.powf(1.5)57            })58            .collect();59        let rough = heights60            .iter()61            .map(|&j| {62                let y = 2.0f64.powf(-j);63                sharp_novelty(&prefix, y) / y64            })65            .collect();66        let line = Line::new();67        let gammas = line.zeros(zeros);68        let coef = line.novelty_coefficients(&gammas);69        Ok(Novelty {70            heights,71            smooth,72            rough,73            gammas,74            coef,75            sieve,76        })77    }78    fn first(&self, count: usize) -> Result<usize, Fault> {79        if count > self.gammas.len() {80            return Err(Fault::new(format!(81                "the meter holds {} zeros.",82                self.gammas.len()83            )));84        }85        Ok(count)86    }87    /// Returns j at every sample, the log base two of one over y.88    pub fn heights(&self) -> Vec<f64> {89        self.heights.clone()90    }91    /// Returns the error at every sample: the smoothed error over y to the three halves, or the sharp window's error over y.92    pub fn dots(&self, sharp: bool) -> Vec<f64> {93        if sharp {94            self.rough.clone()95        } else {96            self.smooth.clone()97        }98    }99    /// Returns the ordinate of every zero the meter holds.100    pub fn gammas(&self) -> Vec<f64> {101        self.gammas.clone()102    }103    /// Returns the modulus of every zero's wave coefficient.104    pub fn amplitudes(&self) -> Vec<f64> {105        self.coef.iter().map(|c| c.abs()).collect()106    }107    /// Sums the waves of the first zeros at the given heights j, on the scale of the dots: over y to the three halves, or over y when sharp, where the waves shrink by the root of y.108    pub fn wave(&self, count: usize, sharp: bool, heights: &[f64]) -> Result<Vec<f64>, Fault> {109        let count = self.first(count)?;110        let log2 = 2.0f64.ln();111        Ok(heights112            .iter()113            .map(|&j| {114                let sum = novelty_wave(&self.gammas[..count], &self.coef[..count], -j * log2);115                if sharp {116                    sum * 2.0f64.powf(-0.5 * j)117                } else {118                    sum119                }120            })121            .collect())122    }123    /// Returns the largest gap between the smoothed dots and the wave of the first zeros, over the largest dot.124    pub fn miss(&self, count: usize) -> Result<f64, Fault> {125        let wave = self.wave(count, false, &self.heights)?;126        let peak = self.smooth.iter().fold(0.0f64, |a, v| a.max(v.abs()));127        let gap = self128            .smooth129            .iter()130            .zip(&wave)131            .map(|(d, w)| (d - w).abs())132            .fold(0.0f64, f64::max);133        Ok(gap / peak)134    }135    /// Returns the reach of the totient sieve.136    pub fn sieve(&self) -> usize {137        self.sieve138    }139}