paper-first-base-below-a-quarter.rs

1use mrlycore::errors::Result;2use mrlycore::MrlyError;3use mrlyfig::out::root;4use mrlyfig::{ink, save, Board};5use std::path::PathBuf;67const NAME: &str = "paper-first-base-below-a-quarter";8const LOW: usize = 10;9const HIGH: usize = 40;10const SUB: usize = 4;11const CAP: f64 = 6.0e6;12const QUARTER: f64 = 0.25;13const RUNGS: usize = 395;1415// LADDER1617fn weights(q: usize, a0: usize, nd: u32, m: usize) -> Vec<f64> {18    let w = q.pow(nd);19    let inv = 1.0 / (q as f64 - 1.0);20    let c = 2.0 * a0 as f64 - q as f64 + 1.0;21    let mass: usize = (0..q).filter(|a| *a != a0).sum();22    let lip = 2.0 * std::f64::consts::PI * mass as f64 * inv;23    let slack = lip / (2.0 * w as f64 * m as f64) + 1e-12;24    let step = 1.0 / (w as f64 * m as f64);25    let mut g = vec![0.0; w];26    for cell in 0..w.div_ceil(2) {27        let mut top = 0.0f64;28        for r in 0..m {29            let t = ((cell * m + r) as f64 + 0.5) * step;30            let d = (std::f64::consts::PI * q as f64 * t).sin() / (std::f64::consts::PI * t).sin();31            let e = d * d - 2.0 * d * (std::f64::consts::PI * c * t).cos() + 1.0;32            let v = if e > 0.0 { e.sqrt() * inv } else { 0.0 };33            top = top.max(v);34        }35        let u = (top + slack).min(1.0);36        g[cell] = u;37        g[w - 1 - cell] = u;38    }39    g40}4142fn sweep(g: &[f64], q: usize, nd: u32, y: &[f64]) -> Vec<f64> {43    let s = q.pow(nd - 1);44    let p = q.pow(nd - 2);45    (0..s)46        .map(|v| {47            let b = (v % p) * q;48            (0..q).map(|c| g[v * q + c] * y[b + c]).sum()49        })50        .collect()51}5253fn power(g: &[f64], q: usize, nd: u32, seed: Option<Vec<f64>>, steps: usize) -> Vec<f64> {54    let s = q.pow(nd - 1);55    let mut y = seed.unwrap_or_else(|| vec![1.0; s]);56    for _ in 0..steps {57        let z = sweep(g, q, nd, &y);58        let top = z.iter().copied().fold(0.0f64, f64::max);59        if top <= 0.0 {60            return y;61        }62        y = z.into_iter().map(|x| x / top).collect();63    }64    y65}6667fn bound(g: &[f64], q: usize, nd: u32, y: &[f64]) -> f64 {68    let z = sweep(g, q, nd, y);69    let mu = z70        .iter()71        .zip(y.iter())72        .map(|(a, b)| a / b)73        .fold(0.0f64, f64::max);74    mu.ln() / (q as f64).ln()75}7677fn lift(y: &[f64], q: usize) -> Vec<f64> {78    (0..y.len() * q).map(|v| y[v / q]).collect()79}8081fn exponent(q: usize, a0: usize) -> f64 {82    let mut nd = 3u32;83    let start = weights(q, a0, nd, SUB);84    let mut y = power(&start, q, nd, None, 300);85    let mut e = bound(&start, q, nd, &y);86    while e >= QUARTER && nd < 5 && (q as f64).powi(nd as i32 + 1) <= CAP {87        let g = weights(q, a0, nd + 1, SUB);88        let z = power(&g, q, nd + 1, Some(lift(&y, q)), 40);89        let step = bound(&g, q, nd + 1, &z);90        let gain = e - step;91        nd += 1;92        y = z;93        e = step;94        if e >= QUARTER && gain <= e - QUARTER {95            break;96        }97    }98    e99}100101fn ladder() -> Vec<(usize, f64)> {102    let mut rungs: Vec<(usize, f64)> = Vec::new();103    for q in LOW..=HIGH {104        for a0 in 0..q.div_ceil(2) {105            rungs.push((q, exponent(q, a0)));106        }107    }108    rungs109}110111// DATA112113fn path() -> PathBuf {114    root()115        .join("files")116        .join("figures")117        .join("data")118        .join(format!("{NAME}.json"))119}120121fn write_data(rungs: &[(usize, f64)]) -> Result<PathBuf> {122    let file = path();123    let folder = file.parent().unwrap().to_path_buf();124    std::fs::create_dir_all(&folder)125        .map_err(|e| MrlyError::Value(format!("cannot make {folder:?}: {e}")))?;126    let bases: Vec<String> = rungs.iter().map(|rung| rung.0.to_string()).collect();127    let exponents: Vec<String> = rungs.iter().map(|rung| format!("{}", rung.1)).collect();128    let text = format!(129        "{{\"base\":[{}],\"exponent\":[{}]}}",130        bases.join(","),131        exponents.join(",")132    );133    std::fs::write(&file, text)134        .map_err(|e| MrlyError::Value(format!("cannot write {file:?}: {e}")))?;135    Ok(file)136}137138fn field(text: &str, name: &str) -> Vec<f64> {139    let head = format!("\"{name}\":[");140    let Some(start) = text.find(&head) else {141        return Vec::new();142    };143    let body = &text[start + head.len()..];144    let end = body.find(']').unwrap_or(0);145    body[..end]146        .split(',')147        .filter_map(|token| token.parse().ok())148        .collect()149}150151fn read_data() -> Result<Vec<(usize, f64)>> {152    let file = path();153    let raw = std::fs::read_to_string(&file)154        .map_err(|e| MrlyError::Value(format!("cannot read {file:?}: {e}; run -- compute")))?;155    let text: String = raw.chars().filter(|c| !c.is_whitespace()).collect();156    let bases = field(&text, "base");157    let exponents = field(&text, "exponent");158    Ok(bases159        .iter()160        .zip(exponents.iter())161        .map(|(base, exponent)| (*base as usize, *exponent))162        .collect())163}164165// PRESS166167fn compute() -> Result<()> {168    let rungs = ladder();169    let below = |q: usize| rungs.iter().filter(|r| r.0 == q && r.1 < QUARTER).count();170    let sets = |q: usize| q.div_ceil(2);171    assert_eq!(rungs.len(), RUNGS, "the ladder has 395 rungs");172    assert!((LOW..21).all(|q| below(q) == 0), "no base under 21 clears");173    assert_eq!(below(21), 1, "base 21 clears at exactly one digit");174    assert!(175        (34..=HIGH).all(|q| below(q) == sets(q)),176        "every base from 34 clears at every digit"177    );178    let gold = rungs.iter().filter(|r| r.1 < QUARTER).count();179    assert_eq!(gold, 163, "163 rungs sit under the quarter");180    let file = write_data(&rungs)?;181    println!(182        "{NAME} {} rungs {gold} under the quarter -> {file:?}",183        RUNGS184    );185    Ok(())186}187188fn draw() -> Result<()> {189    let rungs = read_data()?;190    assert_eq!(rungs.len(), RUNGS);191    let lo = rungs.iter().map(|r| r.1).fold(f64::MAX, f64::min);192    let hi = rungs.iter().map(|r| r.1).fold(f64::MIN, f64::max);193    let pad = (hi - lo) * 0.06;194    let (foot, head) = (lo - pad, hi + pad);195    assert!(196        foot < QUARTER && QUARTER < head,197        "the quarter line is inside"198    );199200    let mut board = Board::square();201    let frame = board.frame(0.08);202    let span = (HIGH - LOW + 1) as f64;203    let slot = frame.w / span;204    let dash = slot * 0.72;205    let place = |e: f64| frame.y + frame.h * (head - e) / (head - foot);206207    let line = place(QUARTER);208    board.rect(209        frame.x,210        line - 1.0,211        frame.w,212        2.0,213        ink::fade(ink::blue(), 0.9),214    );215216    for (q, e) in &rungs {217        let x = frame.x + ((q - LOW) as f64 + 0.5) * slot - dash / 2.0;218        let y = place(*e);219        let (color, thick) = if *e < QUARTER {220            (ink::yellow(), 3.8)221        } else {222            (ink::fade(ink::dim(), 0.7), 2.6)223        };224        board.rect(x, y - thick / 2.0, dash, thick, color);225    }226    save(NAME, &board)?;227    Ok(())228}229230fn main() -> Result<()> {231    if std::env::args().nth(1).as_deref() == Some("compute") {232        return compute();233    }234    draw()235}