design.rs

1use crate::factor::mobius_sieve;2use crate::fft::fft;3use std::f64::consts::PI;45/// The ordinates of the first fourteen nontrivial zeros of the Riemann zeta function, the imaginary parts of the zeros on the critical line in ascending order.6pub const ZETA_ORDINATES: [f64; 14] = [7    14.134725141734693,8    21.022039638771556,9    25.01085758014569,10    30.424876125859512,11    32.93506158773919,12    37.58617815882567,13    40.9187190121475,14    43.327073280915,15    48.00515088116716,16    49.7738324776723,17    52.97032147771446,18    56.44624769706339,19    59.34704400260235,20    60.83177852460981,21];2223/// Returns the digits a bitmask names inside the base, ascending.24///25/// ```26/// assert_eq!(mrlynum::design::digits_of(0b1011, 10), vec![0, 1, 3]);27/// ```28pub fn digits_of(mask: u32, base: u64) -> Vec<u64> {29    (0..base).filter(|&d| mask & (1u32 << d) != 0).collect()30}3132/// Returns the elements of the digit design below the base raised to the depth, ascending: the whole numbers of at most that many base digits, every digit drawn from the set and the leading digit nonzero.33///34/// ```35/// assert_eq!(mrlynum::design::elements(3, &[0, 1], 3), vec![1, 3, 4, 9, 10, 12, 13]);36/// ```37pub fn elements(base: u64, digits: &[u64], depth: usize) -> Vec<u64> {38    let lead: Vec<u64> = digits.iter().copied().filter(|&d| d > 0).collect();39    if depth == 0 || lead.is_empty() {40        return Vec::new();41    }42    let mut out = lead.clone();43    let mut level = lead;44    for _ in 1..depth {45        let mut next = Vec::with_capacity(level.len() * digits.len());46        for value in &level {47            for digit in digits {48                next.push(value * base + digit);49            }50        }51        out.extend_from_slice(&next);52        level = next;53    }54    out55}5657/// Returns the count of elements the design holds at the depth, the length [`elements`] returns without building them.58///59/// ```60/// assert_eq!(mrlynum::design::size(&[0, 1], 20), 1048575);61/// ```62pub fn size(digits: &[u64], depth: usize) -> u128 {63    let lead = digits.iter().filter(|&&d| d > 0).count() as u128;64    let k = digits.len() as u128;65    if depth == 0 || lead == 0 {66        return 0;67    }68    let mut run = 1u128;69    let mut total = 0u128;70    for _ in 0..depth {71        total += lead * run;72        run *= k;73    }74    total75}7677/// Returns the Mobius value of every number by trial division over the primes below the square root of the largest.78///79/// ```80/// assert_eq!(mrlynum::design::mobius_of(&[1, 2, 3, 4, 5, 6]), vec![1, -1, -1, 0, -1, 1]);81/// ```82pub fn mobius_of(values: &[u64]) -> Vec<i8> {83    let top = values.iter().copied().max().unwrap_or(0);84    let root = (top as f64).sqrt() as usize + 2;85    let primes: Vec<u64> = crate::classics::primes(root)86        .into_iter()87        .map(|p| p as u64)88        .collect();89    values90        .iter()91        .map(|&value| {92            let mut rest = value;93            let mut sign = 1i8;94            for &p in &primes {95                if p * p > rest {96                    break;97                }98                if rest % p == 0 {99                    rest /= p;100                    if rest % p == 0 {101                        return 0;102                    }103                    sign = -sign;104                }105            }106            if rest > 1 {107                sign = -sign;108            }109            sign110        })111        .collect()112}113114/// Returns the running design Mobius meter, the partial sums of the Mobius values along the elements.115///116/// ```117/// assert_eq!(mrlynum::design::meter(&[1, -1, -1, 0, -1]), vec![1, 0, -1, -1, -2]);118/// ```119pub fn meter(mu: &[i8]) -> Vec<i64> {120    let mut total = 0i64;121    mu.iter()122        .map(|&value| {123            total += i64::from(value);124            total125        })126        .collect()127}128129/// Returns the log grid uniform over the span of the elements, from the log of the first to the log of the last.130pub fn log_grid(values: &[u64], samples: usize) -> Vec<f64> {131    if values.is_empty() || samples == 0 {132        return Vec::new();133    }134    let lo = (values[0] as f64).ln();135    let hi = (values[values.len() - 1] as f64).ln();136    let step = if samples > 1 {137        (hi - lo) / (samples - 1) as f64138    } else {139        0.0140    };141    (0..samples).map(|i| lo + step * i as f64).collect()142}143144/// Reads the running meter at every point of the log grid and divides by x to the exponent.145pub fn resample(values: &[u64], running: &[i64], exponent: f64, log_x: &[f64]) -> Vec<f64> {146    log_x147        .iter()148        .map(|&t| {149            let slot = values.partition_point(|&v| (v as f64).ln() <= t);150            let held = if slot == 0 {151                0.0152            } else {153                running[slot - 1] as f64154            };155            held / (exponent * t).exp()156        })157        .collect()158}159160/// Returns the density echo, the sum of mu(n) A_F(n)/n over the whole numbers up to each grid point divided by x to the exponent, sieving the Mobius values to the largest element.161pub fn echo_series(values: &[u64], log_x: &[f64], exponent: f64) -> Vec<f64> {162    let top = match values.last() {163        Some(&value) => value,164        None => return vec![0.0; log_x.len()],165    };166    let mu = mobius_sieve(top as usize);167    let mut out = Vec::with_capacity(log_x.len());168    let mut total = 0.0f64;169    let mut seen = 0u64;170    let mut at = 0usize;171    let mut n = 1u64;172    for &t in log_x {173        let mut bound = t.exp().floor() as u64;174        while bound < top && ((bound + 1) as f64).ln() <= t {175            bound += 1;176        }177        while bound > 0 && (bound as f64).ln() > t {178            bound -= 1;179        }180        while n <= bound {181            if at < values.len() && values[at] == n {182                seen += 1;183                at += 1;184            }185            total += f64::from(mu[n as usize]) * seen as f64 / n as f64;186            n += 1;187        }188        out.push(total / (exponent * t).exp());189    }190    out191}192193/// Returns the frequency axis and the power spectrum of the series: the mean removed, a Hann window laid on, a real transform taken, and bin j read as the ordinate 2 pi j over the log range.194pub fn spectrum(log_x: &[f64], series: &[f64]) -> (Vec<f64>, Vec<f64>) {195    let n = series.len();196    if n < 2 || !n.is_power_of_two() {197        return (Vec::new(), Vec::new());198    }199    let step = (log_x[n - 1] - log_x[0]) / (n - 1) as f64;200    let mean = series.iter().sum::<f64>() / n as f64;201    let mut re: Vec<f64> = series202        .iter()203        .enumerate()204        .map(|(i, &value)| {205            let window = 0.5 - 0.5 * (2.0 * PI * i as f64 / (n - 1) as f64).cos();206            (value - mean) * window207        })208        .collect();209    let mut im = vec![0.0f64; n];210    fft(&mut re, &mut im, false);211    let range = n as f64 * step;212    let bins = n / 2 + 1;213    let gamma = (0..bins).map(|j| 2.0 * PI * j as f64 / range).collect();214    let power = (0..bins).map(|j| re[j] * re[j] + im[j] * im[j]).collect();215    (gamma, power)216}217218/// Returns the running median of the power over a window of the given width, the window clamped at the ends.219pub fn median_floor(power: &[f64], width: usize) -> Vec<f64> {220    let n = power.len();221    if n == 0 || width == 0 {222        return vec![0.0; n];223    }224    let half = width / 2;225    let mut window = vec![0.0f64; width];226    (0..n)227        .map(|i| {228            for (slot, cell) in window.iter_mut().enumerate() {229                let at = (i + slot).saturating_sub(half).min(n - 1);230                *cell = power[at];231            }232            let mid = width / 2;233            let (_, value, _) =234                window.select_nth_unstable_by(mid, |a, b| a.partial_cmp(b).unwrap());235            *value236        })237        .collect()238}239240/// Returns the power over its local median floor, the score a peak is read against.241pub fn score(power: &[f64], width: usize) -> Vec<f64> {242    let floor = median_floor(power, width);243    power244        .iter()245        .zip(floor.iter())246        .map(|(&value, &base)| value / base.max(f64::MIN_POSITIVE))247        .collect()248}249250/// Returns the bins inside the band that rise above both neighbours and clear the score threshold, strongest first.251pub fn peaks(gamma: &[f64], score: &[f64], band: (f64, f64), threshold: f64) -> Vec<usize> {252    let mut found: Vec<usize> = (1..score.len().saturating_sub(1))253        .filter(|&i| {254            gamma[i] > band.0255                && gamma[i] < band.1256                && score[i] > score[i - 1]257                && score[i] > score[i + 1]258                && score[i] > threshold259        })260        .collect();261    found.sort_by(|&a, &b| score[b].partial_cmp(&score[a]).unwrap());262    found263}264265/// Returns the design's pole lattice below the top, the ordinates 2 pi j over log q of the poles its Dirichlet series carries.266///267/// ```268/// let lines = mrlynum::design::pole_lattice(3, 13.0);269/// assert_eq!(format!("{:.4} {:.4} {}", lines[0], lines[1], lines.len()), "5.7192 11.4384 2");270/// ```271pub fn pole_lattice(base: u64, top: f64) -> Vec<f64> {272    let step = 2.0 * PI / (base as f64).ln();273    let mut out = Vec::new();274    let mut j = 1;275    while step * (j as f64) < top {276        out.push(step * j as f64);277        j += 1;278    }279    out280}281282/// Returns the distance from the ordinate to the nearest entry of the list, infinite when the list is empty.283pub fn nearest(value: f64, list: &[f64]) -> f64 {284    list.iter()285        .map(|&entry| (entry - value).abs())286        .fold(f64::INFINITY, f64::min)287}288289/// Returns the root mean square of the upper half of the series, the size the echo and the meter are compared at.290pub fn upper_rms(series: &[f64]) -> f64 {291    let half = series.len() / 2;292    let tail = &series[half..];293    if tail.is_empty() {294        return 0.0;295    }296    (tail.iter().map(|value| value * value).sum::<f64>() / tail.len() as f64).sqrt()297}298299#[cfg(test)]300mod tests {301    use super::*;302    use crate::factor::mobius;303304    #[test]305    fn the_design_is_the_digit_strings_in_order() {306        assert_eq!(elements(10, &[0, 1], 2), vec![1, 10, 11]);307        assert_eq!(elements(3, &[1, 2], 2), vec![1, 2, 4, 5, 7, 8]);308        for depth in 1..=8 {309            let built = elements(3, &[0, 1], depth);310            assert_eq!(built.len() as u128, size(&[0, 1], depth));311            assert!(built.windows(2).all(|pair| pair[0] < pair[1]));312        }313        let full = elements(10, &[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], 4);314        assert_eq!(full, (1..10000).collect::<Vec<u64>>());315    }316317    #[test]318    fn the_two_mobius_paths_agree() {319        let values = elements(3, &[0, 1], 12);320        let quick = mobius_of(&values);321        let slow: Vec<i8> = values.iter().map(|&v| mobius(v as usize)).collect();322        assert_eq!(quick, slow);323    }324325    #[test]326    fn the_frequency_axis_is_two_pi_over_the_log_range() {327        let values = elements(3, &[0, 1], 10);328        let running = meter(&mobius_of(&values));329        let log_x = log_grid(&values, 1024);330        let series = resample(&values, &running, 0.5 * 2f64.ln() / 3f64.ln(), &log_x);331        let (gamma, power) = spectrum(&log_x, &series);332        assert_eq!(gamma.len(), 513);333        assert_eq!(power.len(), 513);334        let step = (log_x[1023] - log_x[0]) / 1023.0;335        let range = 1024.0 * step;336        assert_eq!(gamma[0], 0.0);337        for j in [1usize, 7, 100, 512] {338            assert_eq!(gamma[j], 2.0 * PI * j as f64 / range);339        }340    }341}