tube.rs

1use crate::{checked, Fault};2use mrlymath::two;3use wasm_bindgen::prelude::*;45const REACH: usize = 729;6const FLOOR: f64 = 3.0;78// TRANSFORM910fn lower(f: &[f64], d: &mut [f64], hull: &mut [usize], edge: &mut [f64]) {11    let n = f.len();12    let mut k = 0;13    hull[0] = 0;14    edge[0] = f64::NEG_INFINITY;15    edge[1] = f64::INFINITY;16    for q in 1..n {17        let mut cut;18        loop {19            let p = hull[k];20            cut = ((f[q] + (q * q) as f64) - (f[p] + (p * p) as f64)) / (2 * (q - p)) as f64;21            if k > 0 && cut <= edge[k] {22                k -= 1;23            } else {24                break;25            }26        }27        k += 1;28        hull[k] = q;29        edge[k] = cut;30        edge[k + 1] = f64::INFINITY;31    }32    k = 0;33    for (q, slot) in d.iter_mut().enumerate() {34        while edge[k + 1] < q as f64 {35            k += 1;36        }37        let gap = q as f64 - hull[k] as f64;38        *slot = gap * gap + f[hull[k]];39    }40}4142fn transform(types: &[u8], side: usize) -> Vec<f64> {43    let far = (4 * side * side) as f64;44    let mut square: Vec<f64> = types45        .iter()46        .map(|&b| if b != 0 { 0.0 } else { far })47        .collect();48    let mut lane = vec![0.0; side];49    let mut out = vec![0.0; side];50    let mut hull = vec![0usize; side];51    let mut edge = vec![0.0; side + 1];52    for col in 0..side {53        for row in 0..side {54            lane[row] = square[row * side + col];55        }56        lower(&lane, &mut out, &mut hull, &mut edge);57        for row in 0..side {58            square[row * side + col] = out[row];59        }60    }61    for row in 0..side {62        lane.copy_from_slice(&square[row * side..(row + 1) * side]);63        lower(&lane, &mut out, &mut hull, &mut edge);64        square[row * side..(row + 1) * side].copy_from_slice(&out);65    }66    square.iter().map(|v| v.sqrt()).collect()67}6869// FIELD7071struct Field {72    side: usize,73    number: usize,74    dimension: f64,75    dist: Vec<f64>,76}7778impl Field {79    fn read(code: &str, number: usize, level: usize, base: usize) -> Result<Field, Fault> {80        if number < 2 {81            return Err(Fault::new("the side must be at least two."));82        }83        if level < 1 {84            return Err(Fault::new("the level must be at least one."));85        }86        let side = number87            .checked_pow(level as u32)88            .filter(|&side| side <= REACH)89            .ok_or_else(|| {90                Fault::new(format!(91                    "side {number} at level {level} passes the {REACH} cells a side the tube allows."92                ))93            })?;94        let code = checked(code, 2, base)?;95        let tile = two::create(code, number, 1, 0, base)?;96        let digits = tile.types().bytes().iter().filter(|&&b| b != 0).count();97        if digits == 0 {98            return Err(Fault::new("the empty design has no tube."));99        }100        let cell = two::create(code, number, level, 0, base)?;101        if cell.width() != side || cell.height() != side {102            return Err(Fault::new("the design is not a square grid."));103        }104        let dist = transform(cell.types().bytes(), side);105        Ok(Field {106            side,107            number,108            dimension: (digits as f64).ln() / (number as f64).ln(),109            dist,110        })111    }112113    fn volume(&self, eps: f64) -> f64 {114        self.dist115            .iter()116            .map(|&d| (eps - d + 1.0).clamp(0.0, 1.0))117            .sum()118    }119}120121// EXPORTS122123/// Reads the distance field of the design's level-L grid: the Euclidean distance in cell widths from every cell centre to the nearest filled cell centre, row-major, zero on the filled cells.124///125/// The transform is exact, the two-pass lower envelope of parabolas, so a cell inside a square126/// hole of side `s` reads the whole number of cells to the hole's wall. The grid is refused past127/// 729 cells a side.128#[wasm_bindgen]129pub fn tube_distance(130    code: &str,131    number: usize,132    level: usize,133    base: usize,134) -> Result<Vec<f32>, Fault> {135    let field = Field::read(code, number, level, base)?;136    Ok(field.dist.iter().map(|&v| v as f32).collect())137}138139/// Measures the inner tube at radius `eps_cells`: the area within that distance of the design, in cell widths squared, the filled cells counted whole.140///141/// A cell carries the share of itself the radius reaches, `eps - dist + 1` clamped to the unit142/// interval, which is exact at whole radii on a design whose holes are squares: the carpet at143/// level 6 and radius 21 cells reads `5912/6561` of the unit square.144#[wasm_bindgen]145pub fn tube_volume(146    code: &str,147    number: usize,148    level: usize,149    base: usize,150    eps_cells: f64,151) -> Result<f64, Fault> {152    Ok(Field::read(code, number, level, base)?.volume(eps_cells))153}154155/// Walks the Minkowski profile across the radii the grid resolves: the pairs `u = ln(1/eps)` and `M = eps^(d-2) V(eps)` in unit-square units, flat, two numbers a sample.156///157/// The radius runs from `1/q` down to three cells, evenly in `u`, and the pairs come back empty158/// when the grid is too coarse to hold that range. `d` is `log_q k`, the design's fractal159/// dimension, so `M` is the Minkowski content reading whose limit the closed profile is.160#[wasm_bindgen]161pub fn tube_profile(162    code: &str,163    number: usize,164    level: usize,165    base: usize,166    samples: usize,167) -> Result<Vec<f32>, Fault> {168    let field = Field::read(code, number, level, base)?;169    let side = field.side as f64;170    let top = side / field.number as f64;171    if top <= FLOOR || samples < 2 {172        return Ok(Vec::new());173    }174    let (low, high) = ((side / top).ln(), (side / FLOOR).ln());175    let mut pairs = Vec::with_capacity(2 * samples);176    for step in 0..samples {177        let u = low + (high - low) * step as f64 / (samples - 1) as f64;178        let cells = side * (-u).exp();179        let eps = cells / side;180        let volume = field.volume(cells) / (side * side);181        pairs.push(u as f32);182        pairs.push((eps.powf(field.dimension - 2.0) * volume) as f32);183    }184    Ok(pairs)185}186187/// Reads the closed limit profile `G(t)` of the interior-hole class at base `q` with `k` filled cells, the hole sum in closed form.188///189/// Every hole of such a design is an isolated open square whose boundary lies in the design, so190/// the inner tube is the exact sum over the levels, `V(eps) = sum k^(m-1) h(q^-m, eps)`, and191/// `M(eps) = eps^(d-2) V(eps)` runs onto `G`, log-periodic with period `ln q`. The carpet is192/// `q = 3`, `k = 8`, where `G(1/3) = G(1) = 379/280` and the swing above zero is the whole proof193/// that the design is not Minkowski measurable.194#[wasm_bindgen]195pub fn tube_closed(q: usize, k: usize, t: f64) -> Result<f64, Fault> {196    if q < 2 || k <= q || k >= q * q {197        return Err(Fault::new(format!(198            "the class needs {q} < k < {}, not {k}.",199            q * q200        )));201    }202    if !(t > 0.0 && t <= 1.0) {203        return Err(Fault::new("the phase must lie in zero to one."));204    }205    let (base, mass) = (q as f64, k as f64);206    let mut step = 0i32;207    while t * base.powi(step) < 0.5 && step < 512 {208        step += 1;209    }210    while step > -512 && t * base.powi(step - 1) >= 0.5 {211        step -= 1;212    }213    let tail = (mass / (base * base)).powi(step - 1);214    let phase = t * base.powi(step) / base;215    let dimension = mass.ln() / base.ln();216    let sum = 1.0 / (base * base - mass) + 4.0 * phase / (mass - base)217        - 4.0 * phase * phase / (mass - 1.0);218    Ok(t.powf(dimension - 2.0) * tail * sum)219}220221/// Decides whether the design's holes are the isolated interior squares the closed profile needs: every empty cell of the level-one tile off the border, and no two of them touching, even at a corner.222#[wasm_bindgen]223pub fn tube_class(code: &str, number: usize, base: usize) -> Result<bool, Fault> {224    if number < 3 {225        return Ok(false);226    }227    let tile = two::create(checked(code, 2, base)?, number, 1, 0, base)?;228    let types = tile.types().bytes().to_vec();229    if types.len() != number * number {230        return Ok(false);231    }232    let holes: Vec<(usize, usize)> = (0..types.len())233        .filter(|&at| types[at] == 0)234        .map(|at| (at / number, at % number))235        .collect();236    if holes.is_empty() {237        return Ok(false);238    }239    if holes240        .iter()241        .any(|&(row, col)| row == 0 || col == 0 || row + 1 == number || col + 1 == number)242    {243        return Ok(false);244    }245    for (at, &(row, col)) in holes.iter().enumerate() {246        for &(other, next) in &holes[at + 1..] {247            if row.abs_diff(other) <= 1 && col.abs_diff(next) <= 1 {248                return Ok(false);249            }250        }251    }252    Ok(true)253}