star.rs

1use crate::bang::code_to_corners;2use mrlycore::errors::{value_error, Result};3use mrlynum::series::{CATALAN, EULER};45/// The exact reading of a cut layer: how many cells were inked out of how many were read.6#[derive(Clone, Copy, Debug, PartialEq, Eq)]7pub struct Share {8    /// The count of inked cells.9    pub inked: i64,10    /// The count of cells read.11    pub cells: i64,12}1314impl Share {15    /// The share as a real number.16    pub fn value(self) -> f64 {17        self.inked as f64 / self.cells as f6418    }1920    /// The share in lowest terms, numerator then denominator.21    pub fn reduced(self) -> (i64, i64) {22        let (mut a, mut b) = (self.inked.abs(), self.cells.abs());23        while b != 0 {24            let next = a % b;25            a = b;26            b = next;27        }28        let divisor = a.max(1);29        (self.inked / divisor, self.cells / divisor)30    }31}3233/// The real character mod 8 of `Q(sqrt 2)`: `+1` at `n = 1, 7`, `-1` at `n = 3, 5`, zero at even `n`.34///35/// ```36/// use mrlymath::six::star::chi8;37/// assert_eq!([chi8(1), chi8(3), chi8(5), chi8(7), chi8(4)], [1, -1, -1, 1, 0]);38/// ```39pub fn chi8(number: usize) -> i64 {40    [0, 1, 0, -1, 0, -1, 0, 1][number % 8]41}4243/// The closed form of the star arm's ink at odd `n`, `1/2 + chi_8(n)/(2n)`, as `n + chi_8(n)` cells of `2n`.44///45/// The arm `x = y` of the cut is a diameter of `2n` cells with no width to choose, and on it the two46/// block parities agree, so the space rule collapses to `floor(x/4)` even and the inked count is47/// `f(3n) - f(n) = n + chi_8(n)` with `f(m) = 4 floor(m/8) + min(m mod 8, 4)`.48///49/// ```50/// use mrlymath::six::star::arm_law;51/// assert_eq!(arm_law(7).unwrap().reduced(), (4, 7));52/// ```53pub fn arm_law(number: usize) -> Result<Share> {54    if number == 0 || number.is_multiple_of(2) {55        return value_error("the layer number must be odd.");56    }57    Ok(Share {58        inked: number as i64 + chi8(number),59        cells: 2 * number as i64,60    })61}6263/// The constant beside the decay, `ln(1 + sqrt 2)/(2 sqrt 2) - G/8 - gamma/4 - (ln 2)/2`.64///65/// Every limb is a closed form the subtraction of the two ink laws hands over: the character part is66/// `L(1, chi_8)/2` by the class number formula for `Q(sqrt 2)`, the square tail is `-G/8` with `G`67/// Catalan's constant, and the harmonic tail carries `-gamma/4 - (ln 2)/2`.68pub fn constant() -> f64 {69    (1.0 + 2f64.sqrt()).ln() / (2.0 * 2f64.sqrt()) - CATALAN / 8.0 - EULER / 4.0 - 2f64.ln() / 2.070}7172/// The cell-frame decay coefficient of a band of half-width `W` cells, `-(K + b)/(4(2K + 1))` for `K = floor(W/2)`.73///74/// `b` is the block parity of the band edge, one when `floor(K/2)` is even and zero when it is odd.75/// The band's per-layer excess is `kappa chi + (m + q chi_8(n))/n + chi/(8 n^2)`, and summing the76/// `1/n` term over the first `L` odd layers hands this coefficient as `m/2`, so it is a theorem and77/// not a fit. At even `W` it reads `-(W + 2b)/(8(W + 1))` and at odd `W` it reads78/// `-(W - 1 + 2b)/(8W)`, since `x - y` is even on the cut and an odd width is never a new band:79/// `W = 1` returns the arm's `-1/4` and `W = 3` returns `-1/6`. As `W` grows it tends to `-1/8`.80///81/// ```82/// use mrlymath::six::star::width_law;83/// assert_eq!([width_law(0), width_law(2), width_law(4)], [-0.25, -1.0 / 6.0, -0.1]);84/// assert_eq!([width_law(1), width_law(3)], [width_law(0), width_law(2)]);85/// ```86pub fn width_law(half: usize) -> f64 {87    let reach = (half / 2) as f64;88    let parity = f64::from(u8::from((half / 4).is_multiple_of(2)));89    -(reach + parity) / (4.0 * (2.0 * reach + 1.0))90}9192/// The three classes of layer count the `1/L^2` term of the decay reads.93#[derive(Clone, Copy, Debug, PartialEq, Eq)]94pub enum Branch {95    /// `L` divisible by four, where the residual is `-23/192`.96    Zero,97    /// `L` twice an odd number, where the residual is `+25/192`.98    Two,99    /// Odd `L`, where the constant shifts by `1/8` and the error is `O(1/L)`.100    Odd,101}102103impl Branch {104    /// The branch of a layer count.105    pub fn of(layers: usize) -> Branch {106        match layers % 4 {107            0 => Branch::Zero,108            2 => Branch::Two,109            _ => Branch::Odd,110        }111    }112113    /// The exact `1/L^2` coefficient at even `L`, absent at odd `L`.114    pub fn residual(self) -> Option<f64> {115        match self {116            Branch::Zero => Some(-23.0 / 192.0),117            Branch::Two => Some(25.0 / 192.0),118            Branch::Odd => None,119        }120    }121122    /// The constant the ladder converges on, `C` at even `L` and `C + 1/8` at odd `L`.123    pub fn constant(self) -> f64 {124        constant() + if self == Branch::Odd { 0.125 } else { 0.0 }125    }126127    /// The name of the branch.128    pub fn name(self) -> &'static str {129        match self {130            Branch::Zero => "0 mod 4",131            Branch::Two => "2 mod 4",132            Branch::Odd => "odd",133        }134    }135}136137/// The `L`-layer reading of the ghost star's decay in the cell frame.138#[derive(Clone, Copy, Debug, PartialEq)]139pub struct Decay {140    /// The layer count `L`.141    pub layers: usize,142    /// The mean excess of the star over the background across the first `L` odd layers.143    pub excess: f64,144    /// That excess times `L`.145    pub scaled: f64,146    /// The scaled excess plus `(ln L)/4`, which settles on the branch constant.147    pub logged: f64,148    /// The miss of the settled value against the branch constant.149    pub miss: f64,150    /// The miss times `L`, the reading odd `L` leaves behind.151    pub linear: f64,152    /// The miss times `L` squared, the reading even `L` leaves behind.153    pub residual: f64,154    /// The slope of the scaled excess against `ln L`, read from `L/2` to `L`, absent unless `L` is divisible by four.155    pub slope: Option<f64>,156}157158/// The ghost star of a coded cube's hexagonal cut stack, read in the cell frame.159///160/// The cut of the level-one cube of side `n` lives on the cube's own `4n` cell lattice, one layer at161/// a time and nothing resampled, so the half-cell displacement `1/(2n)` between consecutive layers162/// is carried rather than averaged away. A cell `(x, y, z)` of the plane `x + y + z = 6n - 2` is163/// inked when the design's corner mask holds the three block parities `floor(c/4) mod 2`.164pub struct Star {165    corners: [bool; 8],166}167168impl Star {169    /// Reads the star of a base-2 space code, the carpet being `23`.170    pub fn new(code: u128) -> Result<Star> {171        let filled = code_to_corners(code, 3, 2)?;172        let mut corners = [false; 8];173        for corner in filled {174            corners[usize::from(corner[0]) << 2175                | usize::from(corner[1]) << 1176                | usize::from(corner[2])] = true;177        }178        Ok(Star { corners })179    }180181    fn inked(&self, x: i64, y: i64, z: i64) -> bool {182        let bit = |c: i64| (c.div_euclid(4) & 1) as usize;183        self.corners[bit(x) << 2 | bit(y) << 1 | bit(z)]184    }185186    /// The ink of the cut cell at column `x` and even height `z` of the layer at odd `n`.187    ///188    /// The answer is absent where the plane `x + y + z = 6n - 2` leaves the cube, which is what draws189    /// the hexagon's own edge.190    pub fn cell(&self, number: usize, x: i64, z: i64) -> Option<bool> {191        let n = number as i64;192        let y = 6 * n - 2 - x - z;193        let inside = (0..4 * n).contains(&x) && (0..4 * n).contains(&z) && (0..4 * n).contains(&y);194        inside.then(|| self.inked(x, y, z))195    }196197    /// The exact ink share of the whole hexagonal cut at odd `n`, the background the star is read against.198    ///199    /// The rows are the `2n` even heights `z`, and each runs over the `x` the plane leaves inside the200    /// cube. The fill depends on `x` only through `x mod 8`, so a row is counted in eight steps and201    /// the whole layer in `O(n)`.202    ///203    /// ```204    /// use mrlymath::six::star::Star;205    /// assert_eq!(Star::new(23).unwrap().hexagon(3).unwrap().reduced(), (7, 9));206    /// ```207    pub fn hexagon(&self, number: usize) -> Result<Share> {208        if number == 0 || number.is_multiple_of(2) {209            return value_error("the layer number must be odd.");210        }211        let n = number as i64;212        let size = 4 * n;213        let (mut cells, mut inked) = (0i64, 0i64);214        for index in 0..2 * n {215            let z = 2 * index;216            let target = 6 * n - 2 - z;217            let low = 0.max(target - size + 1);218            let high = (size - 1).min(target);219            if high < low {220                continue;221            }222            for rest in 0..8 {223                let count = (high - rest).div_euclid(8) - (low - 1 - rest).div_euclid(8);224                if count <= 0 {225                    continue;226                }227                cells += count;228                if self.inked(rest, (target - rest).rem_euclid(8), z) {229                    inked += count;230                }231            }232        }233        Ok(Share { inked, cells })234    }235236    /// The exact ink share of the band of half-width `W` cells about the arm `x = y` at odd `n`.237    ///238    /// At `W = 0` the band is the arm itself, the `2n` cells of the diameter, whose share is the239    /// closed form of [`arm_law`].240    ///241    /// ```242    /// use mrlymath::six::star::{arm_law, Star};243    /// assert_eq!(Star::new(23).unwrap().arm(9, 0).unwrap(), arm_law(9).unwrap());244    /// ```245    pub fn arm(&self, number: usize, half: usize) -> Result<Share> {246        if number == 0 || number.is_multiple_of(2) {247            return value_error("the layer number must be odd.");248        }249        let (n, half) = (number as i64, half as i64);250        let size = 4 * n;251        let (mut cells, mut inked) = (0i64, 0i64);252        for index in 0..2 * n {253            let z = 2 * index;254            let target = 6 * n - 2 - z;255            let base = target / 2;256            for x in base - half - 1..=base + half + 1 {257                let y = target - x;258                if !(0..size).contains(&x) || !(0..size).contains(&y) || (x - y).abs() > half {259                    continue;260                }261                cells += 1;262                inked += i64::from(self.inked(x, y, z));263            }264        }265        Ok(Share { inked, cells })266    }267268    /// The per-layer excess of the star band over the hexagon across the first `L` odd layers.269    pub fn excesses(&self, layers: usize, half: usize) -> Result<Vec<f64>> {270        if layers == 0 {271            return value_error("the layer count must be at least one.");272        }273        (0..layers)274            .map(|step| {275                let number = 2 * step + 1;276                Ok(self.arm(number, half)?.value() - self.hexagon(number)?.value())277            })278            .collect()279    }280}281282/// The decay read off the per-layer excesses at a layer count, the slope taken from `L/2` to `L`.283///284/// At `W = 0` and even `L` the reading is `L * excess_L = -(ln L)/4 + C + O(1/L^2)`, so the slope285/// settles on `-1/4` and the residual on the branch's exact `1/L^2` coefficient. The sliding window286/// carries a residue trap: the per-layer excess opens on a term in the ink law's character whose sum287/// vanishes at even `L` alone, so a window from `L/2` to `L` cancels it only when `L` is divisible by288/// four. At `L = 2 mod 4` the half is odd and the window converges somewhere else entirely, so the289/// slope is refused there rather than reported wrong.290pub fn decay(excesses: &[f64], layers: usize) -> Result<Decay> {291    if layers < 2 || layers > excesses.len() {292        return value_error("the layer count must be at least two and within the excesses read.");293    }294    let mean = |count: usize| excesses[..count].iter().sum::<f64>() / count as f64;295    let scaled = mean(layers) * layers as f64;296    let near = mean(layers / 2) * (layers / 2) as f64;297    let logged = scaled + (layers as f64).ln() / 4.0;298    let miss = logged - Branch::of(layers).constant();299    Ok(Decay {300        layers,301        excess: mean(layers),302        scaled,303        logged,304        miss,305        linear: miss * layers as f64,306        residual: miss * layers as f64 * layers as f64,307        slope: layers308            .is_multiple_of(4)309            .then(|| (scaled - near) / 2f64.ln()),310    })311}312313#[cfg(test)]314mod tests {315    use super::*;316317    #[test]318    fn the_arm_ink_is_the_closed_form() {319        let star = Star::new(23).unwrap();320        let read: Vec<(i64, i64)> = (0..8)321            .map(|step| star.arm(2 * step + 1, 0).unwrap().reduced())322            .collect();323        assert_eq!(324            read,325            vec![326                (1, 1),327                (1, 3),328                (2, 5),329                (4, 7),330                (5, 9),331                (5, 11),332                (6, 13),333                (8, 15)334            ]335        );336        for step in 0..8 {337            let number = 2 * step + 1;338            assert_eq!(star.arm(number, 0).unwrap(), arm_law(number).unwrap());339        }340        assert!((0..1001).all(|step| {341            let number = 2 * step + 1;342            star.arm(number, 0).unwrap() == arm_law(number).unwrap()343        }));344    }345346    #[test]347    fn the_hexagon_ink_is_the_cut_ink_law() {348        let star = Star::new(23).unwrap();349        for step in 0..28 {350            let number = 2 * step + 1;351            let n = number as i64;352            let chi = if (3 * n - 1) / 2 % 2 == 0 { 1 } else { -1 };353            let want =354                0.5 + chi as f64 / 8.0 + 1.0 / (2 * n) as f64 - chi as f64 / (8 * n * n) as f64;355            assert!((star.hexagon(number).unwrap().value() - want).abs() < 1e-12);356        }357    }358359    #[test]360    fn the_cell_frame_decay_is_minus_a_quarter() {361        let star = Star::new(23).unwrap();362        let excesses = star.excesses(400, 0).unwrap();363        let read = decay(&excesses, 28).unwrap();364        assert_eq!(format!("{:.6}", read.scaled), "-1.126964");365        assert_eq!(format!("{:.10}", read.logged), "-0.2939128437");366        assert_eq!(format!("{:.8}", read.residual), "-0.11937029");367        let deep = decay(&excesses, 400).unwrap();368        assert_eq!(format!("{:.6}", deep.scaled), "-1.791627");369        assert_eq!(format!("{:.10}", deep.logged), "-0.2937613344");370        assert_eq!(format!("{:.8}", deep.residual), "-0.11978958");371        assert!((deep.slope.expect("a window at L = 0 mod 4") + 0.25).abs() < 1e-4);372        assert_eq!(decay(&excesses, 102).unwrap().slope, None);373        assert_eq!(format!("{:.10}", Branch::Zero.constant()), "-0.2937605857");374        assert_eq!(Branch::of(400).residual(), Some(-23.0 / 192.0));375        assert_eq!(Branch::of(402).residual(), Some(25.0 / 192.0));376        assert_eq!(Branch::of(401).name(), "odd");377    }378379    #[test]380    fn a_wider_band_walks_off_the_quarter() {381        let star = Star::new(23).unwrap();382        assert_eq!(width_law(1), width_law(0));383        assert_eq!(width_law(3), width_law(2));384        assert_eq!(format!("{:.6}", width_law(3)), "-0.166667");385        for half in [1usize, 2, 3, 4, 6] {386            let excesses = star.excesses(400, half).unwrap();387            let slope = decay(&excesses, 400)388                .unwrap()389                .slope390                .expect("a window at L = 0 mod 4");391            assert!((slope - width_law(half)).abs() < 2e-4);392        }393    }394395    #[test]396    fn even_layers_are_refused_and_the_code_is_checked() {397        let star = Star::new(23).unwrap();398        assert!(star.hexagon(4).is_err());399        assert!(star.arm(0, 0).is_err());400        assert!(Star::new(1 << 9).is_err());401    }402}