paper-erdos-125-upper-density.rs

1use figures::{ink, save, Board};2use mrlyrs::core::error::Result;3use mrlyrs::num::sumset::Pair;45const LEVEL: u32 = 6;6const SAMPLES: usize = 3072;78// THE MEASURE910fn cantor(mut y: f64) -> f64 {11    let mut out = 0.0;12    let mut w = 1.0;13    for _ in 0..40 {14        if y < 0.0 {15            return out;16        }17        if y >= 0.5 {18            return out + w;19        }20        y *= 3.0;21        w /= 2.0;22        if y >= 0.5 {23            out += w;24            y -= 1.0;25        }26    }27    out + w / 2.028}2930fn energy(tau: f64, level: u32) -> f64 {31    let t = tau * 3f64.powi(level as i32);32    let digits = (t.log(4.0).ceil() as i32).max(0) + 7;33    let mut atoms = vec![0.0f64];34    for l in 1..=digits {35        let step = t * 4f64.powi(-l);36        let more: Vec<f64> = atoms.iter().map(|c| c + step).collect();37        atoms.extend(more);38    }39    let tail = t * 4f64.powi(-digits) / 6.0;40    let weight = 1.0 / atoms.len() as f64;41    let mut p = vec![0.0f64; (t / 3.0).ceil() as usize + 3];42    for c in atoms {43        let c = c + tail;44        let n = c.floor();45        let g = cantor(n + 1.0 - c);46        p[n as usize] += weight * g;47        p[n as usize + 1] += weight * (1.0 - g);48    }49    for j in 0..level {50        let s = 3usize.pow(j);51        let mut q = vec![0.0f64; p.len() + s];52        for (i, v) in p.iter().enumerate() {53            q[i] += v / 2.0;54            q[i + s] += v / 2.0;55        }56        p = q;57    }58    3f64.powi(level as i32) * p.iter().map(|v| v * v).sum::<f64>()59}6061fn chain(k: u32) -> Pair {62    let p = 3u64.pow(k);63    let mut m = 0;64    while 4u64.pow(m) < p {65        m += 1;66    }67    Pair { three: k, four: m }68}6970// THE CHECKS7172fn lattice(pair: Pair) -> f64 {73    3f64.powi(pair.three as i32) * pair.energy() as f64 / 4f64.powi((pair.three + pair.four) as i32)74}7576fn main() -> Result<()> {77    let taus: Vec<f64> = (0..=SAMPLES)78        .map(|i| 4f64.powf(i as f64 / SAMPLES as f64))79        .collect();80    let curves: Vec<Vec<f64>> = (0..=LEVEL)81        .map(|k| taus.iter().map(|&t| energy(t, k)).collect())82        .collect();83    assert_eq!(curves[0][0], 1.0);84    assert_eq!(curves[0][SAMPLES], 0.5);85    for k in 0..LEVEL as usize {86        assert!((0..=SAMPLES).all(|i| curves[k + 1][i] >= curves[k][i] - 1e-12));87    }88    for hs in &curves {89        let area: f64 = (0..SAMPLES)90            .map(|i| (taus[i + 1] - taus[i]) * (hs[i] + hs[i + 1]) / 2.0)91            .sum();92        assert!(area > 2.0 && area < 4.0);93    }9495    let orbit: Vec<(f64, f64)> = (0..=LEVEL)96        .map(|k| {97            if k == 0 {98                return (1.0, energy(1.0, 0));99            }100            let pair = chain(k);101            let tau = pair.scale();102            assert!((1.0..4.0).contains(&tau));103            let h = energy(tau, k);104            assert!((h - lattice(pair)).abs() < 1e-12 * h);105            (tau, h)106        })107        .collect();108    for k in 0..LEVEL as usize {109        let step = orbit[k + 1].1 / orbit[k].1;110        assert!((9.0 / 16.0..=4.5).contains(&step));111    }112    let least = chain(6);113    assert_eq!(least.four, 5);114    assert!((least.ratio(least.energy()) - 858849.0 / 524288.0).abs() < 1e-15);115116    let mut board = Board::square();117    let frame = board.frame(0.08);118    let lo = 0.5 * 0.92;119    let hi = curves[LEVEL as usize]120        .iter()121        .copied()122        .fold(0.0f64, f64::max)123        * 1.03;124    let at = |tau: f64, h: f64| {125        (126            frame.x + frame.w * tau.log(4.0),127            frame.y + frame.h * (hi - h) / (hi - lo),128        )129    };130    let foot = frame.y + frame.h;131    for &(tau, h) in &orbit {132        let (x, y) = at(tau, h);133        board.segment((x, foot), (x, y), 1.4, ink::line());134    }135    board.segment((frame.x, foot), (frame.x + frame.w, foot), 1.4, ink::line());136    for (k, hs) in curves.iter().enumerate() {137        let pts: Vec<(f64, f64)> = taus.iter().zip(hs).map(|(&t, &h)| at(t, h)).collect();138        let t = k as f64 / LEVEL as f64;139        let color = ink::mix(ink::dim(), ink::blue(), t * t);140        board.polyline(&pts, 1.3 + 1.3 * t, color);141    }142    for &(tau, h) in &orbit {143        let (x, y) = at(tau, h);144        board.disc(x, y, 9.0, ink::ground());145        board.disc(x, y, 6.5, ink::orange());146    }147    save("paper-erdos-125-upper-density", &board)?;148    Ok(())149}