wiki-mahler-measure.rs

1use figures::{ink, save, Board};2use mrlyrs::core::error::Result;3use std::f64::consts::TAU;45const LEHMER: [f64; 11] = [1.0, 1.0, 0.0, -1.0, -1.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0];6const SALEM: f64 = 1.176_280_818_259_917_5;7const SAMPLES: usize = 4096;8const FINE: usize = 1 << 16;9const STRETCH: f64 = 0.34;10const FLOOR: f64 = 2.0;1112type Z = (f64, f64);1314fn mul(a: Z, b: Z) -> Z {15    (a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0)16}1718fn div(a: Z, b: Z) -> Z {19    let d = b.0 * b.0 + b.1 * b.1;20    ((a.0 * b.0 + a.1 * b.1) / d, (a.1 * b.0 - a.0 * b.1) / d)21}2223fn norm(z: Z) -> f64 {24    z.0.hypot(z.1)25}2627fn eval(z: Z) -> Z {28    LEHMER.iter().rev().fold((0.0, 0.0), |acc, c| {29        let m = mul(acc, z);30        (m.0 + c, m.1)31    })32}3334fn height(t: f64) -> f64 {35    norm(eval((t.cos(), t.sin()))).ln()36}3738fn roots() -> Vec<Z> {39    let degree = LEHMER.len() - 1;40    let seed = (0.4, 0.9);41    let mut zs: Vec<Z> = Vec::with_capacity(degree);42    let mut z = (1.0, 0.0);43    for _ in 0..degree {44        zs.push(z);45        z = mul(z, seed);46    }47    for _ in 0..500 {48        for k in 0..degree {49            let mut den = (1.0, 0.0);50            for j in 0..degree {51                if j != k {52                    den = mul(den, (zs[k].0 - zs[j].0, zs[k].1 - zs[j].1));53                }54            }55            let step = div(eval(zs[k]), den);56            zs[k] = (zs[k].0 - step.0, zs[k].1 - step.1);57        }58    }59    zs60}6162fn squeeze(v: f64) -> f64 {63    if v >= 0.0 {64        v65    } else {66        -FLOOR * (1.0 - (v / FLOOR).exp())67    }68}6970fn main() -> Result<()> {71    let zs = roots();72    let on: Vec<Z> = zs73        .iter()74        .copied()75        .filter(|z| (norm(*z) - 1.0).abs() < 1e-9)76        .collect();77    let off: Vec<Z> = zs78        .iter()79        .copied()80        .filter(|z| (norm(*z) - 1.0).abs() >= 1e-9)81        .collect();82    assert_eq!(on.len(), 8);83    assert_eq!(off.len(), 2);84    let outside: f64 = off.iter().map(|z| norm(*z).max(1.0)).product();85    assert!((outside - SALEM).abs() < 1e-12);86    assert!(off.iter().all(|z| z.1.abs() < 1e-12 && z.0 > 0.0));87    let mean = (0..FINE)88        .map(|i| height((i as f64 + 0.5) / FINE as f64 * TAU))89        .sum::<f64>()90        / FINE as f64;91    assert!((mean - SALEM.ln()).abs() < 1e-4);9293    let reach: Vec<f64> = (0..SAMPLES)94        .map(|i| 1.0 + STRETCH * squeeze(height(i as f64 / SAMPLES as f64 * TAU)))95        .collect();96    let mut angles: Vec<f64> = (0..SAMPLES)97        .map(|i| i as f64 / SAMPLES as f64 * TAU)98        .collect();99    angles.extend(on.iter().map(|z| z.1.atan2(z.0).rem_euclid(TAU)));100    angles.sort_by(|a, b| a.total_cmp(b));101    assert_eq!(angles.len(), SAMPLES + 8);102    let levels: Vec<f64> = angles.iter().map(|t| height(*t)).collect();103    let plane: Vec<(f64, f64)> = angles104        .iter()105        .zip(&levels)106        .map(|(t, v)| {107            let r = 1.0 + STRETCH * squeeze(*v);108            (r * t.cos(), r * t.sin())109        })110        .collect();111    let (mut x0, mut x1, mut y0, mut y1) = (f64::MAX, f64::MIN, f64::MAX, f64::MIN);112    for p in &plane {113        x0 = x0.min(p.0);114        x1 = x1.max(p.0);115        y0 = y0.min(p.1);116        y1 = y1.max(p.1);117    }118119    let mut board = Board::square();120    let frame = board.frame(0.08);121    let scale = frame.w / (x1 - x0).max(y1 - y0);122    let (fx, fy) = frame.center();123    let cx = fx - scale * (x0 + x1) / 2.0;124    let cy = fy + scale * (y0 + y1) / 2.0;125    let at = |x: f64, y: f64| (cx + scale * x, cy - scale * y);126127    let out = ink::fade(ink::blue(), 0.22);128    let inn = ink::fade(ink::orange(), 0.22);129    let (bx0, by0) = at(x0, y1);130    let (bx1, by1) = at(x1, y0);131    for py in (by0.floor() as usize)..(by1.ceil() as usize).min(board.height) {132        for px in (bx0.floor() as usize)..(bx1.ceil() as usize).min(board.width) {133            let x = (px as f64 + 0.5 - cx) / scale;134            let y = (cy - py as f64 - 0.5) / scale;135            let rho = x.hypot(y);136            let u = y.atan2(x).rem_euclid(TAU) / TAU * SAMPLES as f64;137            let i = (u.floor() as usize) % SAMPLES;138            let w = u - u.floor();139            let r = reach[i] * (1.0 - w) + reach[(i + 1) % SAMPLES] * w;140            if rho > 1.0 && rho < r {141                board.blend(px, py, out, 1.0);142            } else if rho < 1.0 && rho > r {143                board.blend(px, py, inn, 1.0);144            }145        }146    }147148    board.ring(cx, cy, scale, 3.0, ink::dim());149    let mut run: Vec<(f64, f64)> = Vec::new();150    let mut up = levels[0] >= 0.0;151    for i in 0..=angles.len() {152        let k = i % angles.len();153        let p = at(plane[k].0, plane[k].1);154        let side = levels[k] >= 0.0;155        if side != up {156            run.push(p);157            board.polyline(&run, 5.0, if up { ink::blue() } else { ink::orange() });158            run.clear();159            up = side;160        }161        run.push(p);162    }163    board.polyline(&run, 5.0, if up { ink::blue() } else { ink::orange() });164165    for z in &on {166        let (px, py) = at(z.0, z.1);167        board.disc(px, py, 8.0, ink::fg());168    }169    for z in &off {170        let (px, py) = at(z.0, z.1);171        board.disc(px, py, 13.0, ink::yellow());172    }173    save("wiki-mahler-measure", &board)?;174    Ok(())175}