main.rs

1use mrlyrs::math::three::sponge::{deep, dimension, distance, exact, profile, COVER, EDGE};2use std::time::Instant;34const PLUS: f64 = 7.0 / 27.0;5const SIDE: usize = 40;6const RUN: usize = 800;7const TOP: usize = 4;8const COARSE: usize = 72;9const CENTRE: usize = 120;10const STEPS: usize = 3000;11const DEPTH: i32 = 40;1213// CROSSING1415struct Hole {16    level: usize,17    start: f64,18    side: f64,19    rows: Vec<f64>,20}2122fn crossing(line: f64, top: usize) -> Vec<Hole> {23    let (mut local, mut start, mut rows) = (3.0 * line, 0.0, vec![0.0]);24    let mut out = Vec::new();25    for level in 1..=top {26        let side = 3f64.powi(-(level as i32) - 1);27        let digit = (3.0 * local).floor().clamp(0.0, 2.0);28        local = 3.0 * local - digit;29        if digit == 1.0 {30            let here = rows.iter().map(|q| q + side).collect();31            out.push(Hole {32                level,33                start: start + side,34                side,35                rows: here,36            });37        }38        let choices: &[f64] = if digit == 1.0 {39            &[0.0, 2.0]40        } else {41            &[0.0, 1.0, 2.0]42        };43        rows = rows44            .iter()45            .flat_map(|q| choices.iter().map(move |b| q + b * side))46            .collect();47        start += digit * side;48    }49    out50}5152fn integrand(a: f64, radius: f64, hole: &Hole, line: f64) -> f64 {53    let (half, reach) = (hole.side / 2.0, line - hole.start);54    let fall = (radius * radius - (hole.start + a).powi(2)).max(0.0).sqrt();55    let width = reach.min(hole.side - fall) - a.max(fall);56    if fall >= half || width <= 0.0 {57        return 0.0;58    }59    4.0 * (half - fall) * width60}6162fn simpson(63    f: &dyn Fn(f64) -> f64,64    a: f64,65    b: f64,66    fa: f64,67    fm: f64,68    fb: f64,69    whole: f64,70    depth: usize,71) -> f64 {72    let m = (a + b) / 2.0;73    let (lm, rm) = ((a + m) / 2.0, (m + b) / 2.0);74    let (flm, frm) = (f(lm), f(rm));75    let left = (m - a) * (fa + 4.0 * flm + fm) / 6.0;76    let right = (b - m) * (fm + 4.0 * frm + fb) / 6.0;77    if depth == 0 || (left + right - whole).abs() <= 1e-22 {78        return left + right + (left + right - whole) / 15.0;79    }80    simpson(f, a, m, fa, flm, fm, left, depth - 1) + simpson(f, m, b, fm, frm, fb, right, depth - 1)81}8283fn quadrature(radius: f64, hole: &Hole, line: f64) -> f64 {84    let f = |a: f64| integrand(a, radius, hole, line);85    let reach = line - hole.start;86    let low = ((radius * radius - hole.side * hole.side / 4.0)87        .max(0.0)88        .sqrt()89        - hole.start)90        .clamp(0.0, reach);91    let cuts = 64;92    let mut total = 0.0;93    for k in 0..cuts {94        let a = low + (reach - low) * k as f64 / cuts as f64;95        let b = low + (reach - low) * (k + 1) as f64 / cuts as f64;96        let (fa, fb, fm) = (f(a), f(b), f((a + b) / 2.0));97        total += simpson(98            &f,99            a,100            b,101            fa,102            fm,103            fb,104            (b - a) * (fa + 4.0 * fm + fb) / 6.0,105            40,106        );107    }108    total109}110111// RASTER112113fn box_low(radius: f64, hole: &Hole, line: f64) -> f64 {114    (radius * radius - hole.side * hole.side / 4.0)115        .max(0.0)116        .sqrt()117        .max(hole.start)118        .min(line)119}120121fn rows_raster(radius: f64, hole: &Hole, line: f64) -> [f64; 3] {122    let low = box_low(radius, hole, line);123    let (step, pace) = ((line - low) / SIDE as f64, hole.side / RUN as f64);124    let half = (2.0 * step * step + pace * pace).sqrt() / 2.0;125    let (mut mid, mut inner, mut outer) = (0.0, 0.0, 0.0);126    for row in &hole.rows {127        for a in 0..SIDE {128            for b in 0..SIDE {129                for c in 0..RUN {130                    let u = 1.0 / 3.0 + low + (a as f64 + 0.5) * step;131                    let v = 1.0 / 3.0 + low + (b as f64 + 0.5) * step;132                    let d = distance([row + (c as f64 + 0.5) * pace, u, v]);133                    mid += f64::from(u8::from(d > radius));134                    inner += f64::from(u8::from(d > radius + half));135                    outer += f64::from(u8::from(d > radius - half));136                }137            }138        }139    }140    let cell = step * step * pace;141    [mid * cell, inner * cell, outer * cell]142}143144fn outside(radius: f64, holes: &[Hole], line: f64) -> f64 {145    let cell = 1.0 / (3.0 * COARSE as f64);146    let mut worst: f64 = 0.0;147    for a in 0..COARSE {148        for b in 0..COARSE {149            for c in 0..COARSE {150                let along = (a as f64 + 0.5) * cell;151                let (y, z) = ((b as f64 + 0.5) * cell, (c as f64 + 0.5) * cell);152                let (u, v) = (y.min(1.0 / 3.0 - y), z.min(1.0 / 3.0 - z));153                if u > line || v > line {154                    continue;155                }156                let inside = holes.iter().any(|h| {157                    let low = box_low(radius, h, line);158                    u >= low && v >= low && h.rows.iter().any(|q| along > *q && along < q + h.side)159                });160                if !inside {161                    worst = worst.max(distance([along, 1.0 / 3.0 + y, 1.0 / 3.0 + z]));162                }163            }164        }165    }166    worst167}168169fn centre_raster(radius: f64) -> [f64; 3] {170    let lap = (radius * radius - EDGE * EDGE).max(0.0).sqrt();171    let step = (EDGE - lap) / CENTRE as f64;172    let half = step * 3f64.sqrt() / 2.0;173    let (mut mid, mut inner, mut outer) = (0.0, 0.0, 0.0);174    for a in 0..CENTRE {175        for b in 0..CENTRE {176            for c in 0..CENTRE {177                let at = [a, b, c].map(|k| 1.0 / 3.0 + lap + (k as f64 + 0.5) * step);178                let d = distance(at);179                mid += f64::from(u8::from(d > radius));180                inner += f64::from(u8::from(d > radius + half));181                outer += f64::from(u8::from(d > radius - half));182            }183        }184    }185    let cell = 8.0 * step.powi(3);186    [mid * cell, inner * cell, outer * cell]187}188189// PERIOD190191fn periodic(radius: f64) -> f64 {192    let (mut total, mut weight, mut reach) = (20.0 / 27.0, 1.0, radius);193    for _ in 0..=DEPTH {194        total += weight * exact(reach);195        weight *= 27.0 / 20.0;196        reach /= 3.0;197    }198    radius.powf(dimension() - 3.0) * total199}200201fn golden(f: &dyn Fn(f64) -> f64, mut a: f64, mut b: f64, sign: f64) -> f64 {202    let ratio = (5f64.sqrt() - 1.0) / 2.0;203    for _ in 0..200 {204        let (x, y) = (b - ratio * (b - a), a + ratio * (b - a));205        if sign * f(x) > sign * f(y) {206            b = y;207        } else {208            a = x;209        }210    }211    (a + b) / 2.0212}213214fn main() {215    let clock = Instant::now();216    println!("THRESHOLDS");217    for m in 1..=5 {218        let half = 3f64.powi(-m - 1) / 2.0;219        println!(220            "delta_{m} = sqrt(1/36 + s_{m}^2/4) = {:.12}",221            (1.0 / 36.0 + half * half).sqrt()222        );223    }224    println!(225        "arms swallowed at sqrt(10)/18 = {:.12}",226        10f64.sqrt() / 18.0227    );228    println!("plus swallowed at sqrt(2)/6 = {:.12}", COVER);229    println!();230    println!(231        "LEVELS: closed deep, lab quadrature by level, row raster by level (midpoint, bracket)"232    );233    let radii = [234        0.22,235        0.2,236        0.18,237        0.174,238        0.17,239        0.168,240        0.1675,241        0.167,242        0.16675,243        0.1667,244        0.16667,245        EDGE,246        0.148,247        1.0 / 8.0,248        1.0 / 12.0,249    ];250    let mut worst_ratio: f64 = 0.0;251    let mut worst_centre: f64 = 0.0;252    let mut worst_gap: f64 = 0.0;253    for &radius in &radii {254        let line = radius.min(EDGE);255        let holes = crossing(line, 12);256        let parts: Vec<f64> = holes257            .iter()258            .map(|h| h.rows.len() as f64 * quadrature(radius, h, line))259            .collect();260        let sum: f64 = parts.iter().sum();261        worst_gap = worst_gap.max((deep(radius) - sum).abs());262        println!(263            "radius {radius:.6}  T = {:.12}  deep = {:.6e}  lab sum = {:.6e}  gap = {:.1e}",264            exact(radius),265            deep(radius),266            sum,267            (deep(radius) - sum).abs()268        );269        for (h, part) in holes.iter().zip(&parts).filter(|(h, _)| h.level <= TOP) {270            if *part == 0.0 {271                println!(272                    "  level {} x {}: closed 0, raster skipped",273                    h.level,274                    h.rows.len()275                );276                continue;277            }278            let [mid, inner, outer] = rows_raster(radius, h, line);279            if radius >= EDGE {280                worst_ratio = worst_ratio.max((mid / part - 1.0).abs());281            }282            assert!(283                inner <= *part && *part <= outer,284                "bracket misses level {} at {radius}",285                h.level286            );287            println!(288                "  level {} x {}: lab {:.6e}  raster {:.6e}  in [{:.6e}, {:.6e}]",289                h.level,290                h.rows.len(),291                part,292                mid,293                inner,294                outer295            );296        }297        if radius >= EDGE {298            let far = outside(radius, &holes, line);299            assert!(300                far <= radius,301                "an uncovered point outside the boxes at {radius}"302            );303            let [mid, inner, outer] = centre_raster(radius);304            let closed = PLUS - exact(radius) - 24.0 * deep(radius);305            assert!(306                inner <= closed && closed <= outer,307                "centre bracket misses at {radius}"308            );309            worst_centre = worst_centre.max((mid / closed - 1.0).abs());310            println!("  arm cells outside the boxes: largest distance {far:.9} <= radius");311            println!("  centre deficit: closed {closed:.9e}  raster {mid:.9e}  in [{inner:.6e}, {outer:.6e}]");312        }313    }314    println!("largest relative gap on [1/6, sqrt(2)/6], row raster midpoint against the level integral: {worst_ratio:.2e}");315    println!(316        "largest gap at all {} radii, crate deep against the lab quadrature: {worst_gap:.1e}",317        radii.len()318    );319    println!("largest relative gap on [1/6, sqrt(2)/6], centre raster midpoint against the closed deficit: {worst_centre:.2e}");320    println!();321    println!("PERIOD: p from exact T on u = ln(1/eps) over one period from ln(1/(sqrt(2)/6))");322    let start = (1.0 / COVER).ln();323    let span = 3f64.ln();324    let at = |u: f64| periodic((-u).exp());325    let grid = |n: usize| -> Vec<(f64, f64)> {326        (0..n)327            .map(|i| start + span * (i as f64 + 0.5) / n as f64)328            .map(|u| (u, at(u)))329            .collect()330    };331    let walk = grid(STEPS);332    let step = span / STEPS as f64;333    let (top_u, _) =334        walk.iter().copied().fold(335            (0.0, f64::NEG_INFINITY),336            |b, x| if x.1 > b.1 { x } else { b },337        );338    let (bottom_u, _) =339        walk.iter()340            .copied()341            .fold((0.0, f64::INFINITY), |b, x| if x.1 < b.1 { x } else { b });342    let crest = golden(&at, top_u - step, top_u + step, 1.0);343    let trough = golden(&at, bottom_u - step, bottom_u + step, -1.0);344    let (high, low) = (at(crest), at(trough));345    println!("maximum p = {high:.9} at eps = {:.9}", (-crest).exp());346    println!("minimum p = {low:.9} at eps = {:.9}", (-trough).exp());347    println!(348        "swing (max - min)/min = {:.6} %",349        100.0 * (high - low) / low350    );351    let mean = |n: usize| grid(n).iter().map(|x| x.1).sum::<f64>() / n as f64;352    println!(353        "logarithmic mean of p over the period: {:.9} ({} steps), {:.9} ({} steps)",354        mean(STEPS),355        STEPS,356        mean(2 * STEPS),357        2 * STEPS358    );359    let window: Vec<f64> = walk360        .iter()361        .filter(|x| (-x.0).exp() > EDGE)362        .map(|x| x.1)363        .collect();364    let window_high = window.iter().copied().fold(f64::NEG_INFINITY, f64::max);365    let window_low = window.iter().copied().fold(f64::INFINITY, f64::min);366    println!("on the phases (1/6, sqrt(2)/6]: p runs over [{window_low:.9}, {window_high:.9}], {} of {STEPS} steps", window.len());367    let falling = window.windows(2).all(|pair| pair[1] > pair[0]);368    println!("p falls as eps grows, at every step across (1/6, sqrt(2)/6]: {falling}");369    for eps in [1.0 / 12.0, 1.0 / 8.0, EDGE, 0.2, COVER] {370        let edge = profile(eps).map_or("none".to_string(), |v| format!("{v:.9}"));371        println!(372            "p({eps:.9}) = {:.9}   upper edge (Deep = 0) {edge}",373            periodic(eps)374        );375    }376    println!();377    println!("runtime {:.1} s", clock.elapsed().as_secs_f64());378}