radix.rs

1use crate::gauss::Ring;2use std::collections::HashSet;3use std::f64::consts::TAU;45fn divides(ring: Ring, z: (i64, i64), w: (i64, i64)) -> bool {6    let n = ring.norm(w.0, w.1) as i64;7    let (x, y) = ring.mul(z, ring.conjugate(w.0, w.1));8    x % n == 0 && y % n == 09}1011fn argument(ring: Ring, a: i64, b: i64) -> f64 {12    let (x, y) = ring.place(a, b);13    let t = y.atan2(x);14    if t < 0.0 {15        t + TAU16    } else {17        t18    }19}2021/// The base of a radix design: a ring and an element of norm at least two, the scale every word is read against.22#[derive(Clone, Copy, Debug, PartialEq, Eq)]23pub struct Base {24    ring: Ring,25    value: (i64, i64),26}2728impl Base {29    /// Fixes a base in a ring, panicking below norm two.30    pub fn new(ring: Ring, value: (i64, i64)) -> Base {31        assert!(32            ring.norm(value.0, value.1) >= 2,33            "a base needs norm two or more"34        );35        Base { ring, value }36    }37    /// Returns the ring.38    pub fn ring(self) -> Ring {39        self.ring40    }41    /// Returns the base element.42    pub fn value(self) -> (i64, i64) {43        self.value44    }45    /// Returns the norm `q` of the base: the count of residue classes and the square of the scale.46    pub fn norm(self) -> u64 {47        self.ring.norm(self.value.0, self.value.1)48    }49    /// Returns the base raised to a level.50    pub fn power(self, level: usize) -> (i64, i64) {51        let mut out = (1, 0);52        for _ in 0..level {53            out = self.ring.mul(out, self.value);54        }55        out56    }57    /// Returns whether two points are congruent modulo the base.58    pub fn congruent(self, z: (i64, i64), w: (i64, i64)) -> bool {59        divides(self.ring, (z.0 - w.0, z.1 - w.1), self.value)60    }61    /// Returns the canonical complete residue system modulo the base: the `q` representatives of least norm, ties broken by argument in `[0, 2 pi)`.62    ///63    /// The system is built greedily over every point of norm at most `q` sorted by norm and then by argument, which reaches every class because the covering radius of the lattice `b R` is `|b| / sqrt(2)` on the square lattice and `|b| / sqrt(3)` on the hexagonal.64    ///65    /// ```66    /// use mrlynum::gauss::Ring;67    /// use mrlynum::radix::Base;68    /// assert_eq!(Base::new(Ring::Gaussian, (1, 1)).residues(), vec![(0, 0), (1, 0)]);69    /// ```70    pub fn residues(self) -> Vec<(i64, i64)> {71        let q = self.norm();72        let reach = (2 * q).isqrt() as i64 + 1;73        let mut pool = Vec::new();74        for a in -reach..=reach {75            for b in -reach..=reach {76                if self.ring.norm(a, b) <= q {77                    pool.push((a, b));78                }79            }80        }81        pool.sort_by(|&(a, b), &(c, d)| {82            let left = (self.ring.norm(a, b), argument(self.ring, a, b));83            let right = (self.ring.norm(c, d), argument(self.ring, c, d));84            left.085                .cmp(&right.0)86                .then_with(|| left.1.total_cmp(&right.1))87        });88        let mut out: Vec<(i64, i64)> = Vec::with_capacity(q as usize);89        for z in pool {90            if out.len() == q as usize {91                break;92            }93            if !out.iter().any(|&w| self.congruent(z, w)) {94                out.push(z);95            }96        }97        assert_eq!(out.len(), q as usize, "the residue system is incomplete");98        out99    }100    /// Returns the index in the canonical residue system of the class of a point.101    pub fn class(self, z: (i64, i64)) -> usize {102        self.residues()103            .iter()104            .position(|&w| self.congruent(z, w))105            .expect("every point lies in a residue class")106    }107    /// Returns whether the conjugate of the base is an associate of the base, which is when the mirror joins the symmetry group.108    pub fn mirrored(self) -> bool {109        let c = self.ring.conjugate(self.value.0, self.value.1);110        self.ring111            .associates(self.value.0, self.value.1)112            .contains(&c)113    }114    /// Returns the symmetry group of the base as permutations of the canonical residue indices: every unit multiplication, and every unit times conjugation when the conjugate of the base is an associate of the base.115    ///116    /// Multiplying every digit by a unit `v` carries the attractor of a design to `v` times that attractor over the same base, so the unit group acts; conjugation carries the base to its conjugate and joins the group exactly when that is an associate.117    ///118    /// The list is the Burnside multiset, one entry per abstract group element, and not the order of the permutation group it induces: the action need not be faithful, so entries repeat, and the acting image is the deduplicated list.119    pub fn group(self) -> Vec<Vec<usize>> {120        let residues = self.residues();121        let mirror = self.mirrored();122        let mut out = Vec::new();123        for unit in self.ring.associates(1, 0) {124            for flip in [false, true] {125                if flip && !mirror {126                    continue;127                }128                let map: Vec<usize> = residues129                    .iter()130                    .map(|&z| {131                        let w = if flip {132                            self.ring.conjugate(z.0, z.1)133                        } else {134                            z135                        };136                        let image = self.ring.mul(unit, w);137                        residues138                            .iter()139                            .position(|&r| self.congruent(image, r))140                            .expect("a unit multiple lands in a class")141                    })142                    .collect();143                out.push(map);144            }145        }146        out147    }148}149150/// A radix design: a digit set inside one ring, placed by a base with a unit twist per digit.151#[derive(Clone, Debug, PartialEq, Eq)]152pub struct Radix {153    base: Base,154    digits: Vec<(i64, i64)>,155    twists: Vec<(i64, i64)>,156}157158impl Radix {159    /// Builds a design from a base, a digit list and a unit twist per digit, panicking on a length mismatch, a twist whose norm is not one, or two digits congruent modulo the base.160    ///161    /// Pairwise incongruent digits are the hypothesis of the untwisted fill law: they are what recovers the last digit from the word read modulo the base, so a repeated class is refused here rather than silently gluing words.162    pub fn new(base: Base, digits: Vec<(i64, i64)>, twists: Vec<(i64, i64)>) -> Radix {163        assert_eq!(digits.len(), twists.len(), "one twist per digit");164        for &(a, b) in &twists {165            assert_eq!(base.ring().norm(a, b), 1, "a twist is a unit");166        }167        for (i, &z) in digits.iter().enumerate() {168            for &w in &digits[..i] {169                assert!(170                    !base.congruent(z, w),171                    "the digits {w:?} and {z:?} are congruent modulo the base"172                );173            }174        }175        Radix {176            base,177            digits,178            twists,179        }180    }181    /// Builds an untwisted design from a code over the canonical residue system: bit `i` of the code selects residue `i`.182    pub fn from_code(base: Base, code: u128) -> Radix {183        let residues = base.residues();184        assert!(185            code >> residues.len() == 0,186            "code out of range for the base"187        );188        let digits: Vec<(i64, i64)> = residues189            .into_iter()190            .enumerate()191            .filter(|(i, _)| (code >> i) & 1 == 1)192            .map(|(_, z)| z)193            .collect();194        let twists = vec![(1, 0); digits.len()];195        Radix::new(base, digits, twists)196    }197    /// Returns the design with the twists named by their index in the unit list, the units in turning order from one.198    pub fn with_twists(self, units: &[usize]) -> Radix {199        let list = self.base.ring().associates(1, 0);200        let twists = units.iter().map(|&i| list[i]).collect();201        Radix::new(self.base, self.digits, twists)202    }203    /// Returns the base.204    pub fn base(&self) -> Base {205        self.base206    }207    /// Returns the ring.208    pub fn ring(&self) -> Ring {209        self.base.ring()210    }211    /// Returns the digits.212    pub fn digits(&self) -> &[(i64, i64)] {213        &self.digits214    }215    /// Returns the twists.216    pub fn twists(&self) -> &[(i64, i64)] {217        &self.twists218    }219    /// Returns the digit count `|F|`.220    pub fn size(&self) -> usize {221        self.digits.len()222    }223    /// Returns the code of the classes the digits occupy, which names the design only when the digits are the canonical representatives.224    pub fn code(&self) -> u128 {225        let mut out = 0u128;226        for &d in &self.digits {227            out |= 1 << self.base.class(d);228        }229        out230    }231    /// Returns whether every digit is the canonical representative of its class.232    pub fn canonical(&self) -> bool {233        let residues = self.base.residues();234        self.digits.iter().all(|d| residues.contains(d))235    }236    /// Returns the count of words of a level, `|F|^L`.237    pub fn fill(&self, level: usize) -> u128 {238        (self.size() as u128).pow(level as u32)239    }240    /// Returns the similarity dimension `log |F| / log sqrt(q)`, the ratio of the digit count to the scale of the base.241    pub fn dimension(&self) -> f64 {242        (self.size() as f64).ln() / (self.base.norm() as f64).ln() * 2.0243    }244    /// Returns the level-`L` points in exact ring coordinates scaled by `b^L`.245    ///246    /// The place map of digit `d` is `phi_d(x) = (u_d x + d) / b`, a word `d_1 ... d_L` lands on `phi_(d_1)(... phi_(d_L)(0))`, and that point times `b^L` is `sum_(j=1..L) (prod_(i<j) u_(d_i)) d_j b^(L-j)`, which stays in the ring. Words come in digit-lexicographic order, the first digit slowest.247    pub fn words(&self, level: usize) -> Vec<(i64, i64)> {248        let ring = self.ring();249        let mut out = vec![(0i64, 0i64)];250        for step in 1..=level {251            let shift = self.base.power(step - 1);252            let mut next = Vec::with_capacity(out.len() * self.size());253            for (&d, &u) in self.digits.iter().zip(self.twists.iter()) {254                let head = ring.mul(d, shift);255                for &p in &out {256                    let t = ring.mul(u, p);257                    next.push((head.0 + t.0, head.1 + t.1));258                }259            }260            out = next;261        }262        out263    }264    /// Returns the level-`L` points in the plane, the scaled words divided by `b^L`.265    pub fn plane(&self, level: usize) -> Vec<(f64, f64)> {266        let ring = self.ring();267        let p = self.base.power(level);268        let (px, py) = ring.place(p.0, p.1);269        let den = px * px + py * py;270        self.words(level)271            .into_iter()272            .map(|(a, b)| {273                let (x, y) = ring.place(a, b);274                ((x * px + y * py) / den, (y * px - x * py) / den)275            })276            .collect()277    }278    /// Returns the count of distinct level-`L` points: the glue count, which is the fill exactly when no two words name one point.279    pub fn distinct(&self, level: usize) -> usize {280        self.words(level).into_iter().collect::<HashSet<_>>().len()281    }282}283284/// Returns the Koch curve as a radix design: base `3` on the hexagonal lattice, digits `0, 1, 2 + omega, 2`, twists `1, e^(i pi/3), e^(-i pi/3), 1`.285///286/// The digits are not the canonical residues: `2` and `-1` share a class and the canonical system holds `-1`, so the code alone does not name this design.287pub fn koch() -> Radix {288    let base = Base::new(Ring::Eisenstein, (3, 0));289    let digits = vec![(0, 0), (1, 0), (2, 1), (2, 0)];290    Radix::new(base, digits, vec![(1, 0); 4]).with_twists(&[0, 1, 5, 0])291}292293/// Returns the Sierpinski gasket as a radix design: base `2` on the hexagonal lattice, three of the four residues, code `7`.294pub fn gasket() -> Radix {295    Radix::from_code(Base::new(Ring::Eisenstein, (2, 0)), 7)296}297298/// Returns the twindragon as a radix design: base `1 + i` on the square lattice, the full residue system, code `3`.299pub fn twindragon() -> Radix {300    Radix::from_code(Base::new(Ring::Gaussian, (1, 1)), 3)301}302303/// Returns the terdragon as a radix design: base `2 + omega` on the hexagonal lattice, the full residue system, code `7`, twisted by `1, omega, 1`.304///305/// The untwisted code misses the curve: reading the L-system `F -> F + F - F` at `120` degrees as a turtle, three segments to a level, and normalising by the endpoint gives the segment starts word for word only under this twist.306pub fn terdragon() -> Radix {307    Radix::from_code(Base::new(Ring::Eisenstein, (2, 1)), 7).with_twists(&[0, 2, 0])308}309310/// Returns the flowsnake as a radix design: base `3 + omega` of norm seven on the hexagonal lattice, the full residue system, code `127`.311pub fn flowsnake() -> Radix {312    Radix::from_code(Base::new(Ring::Eisenstein, (3, 1)), 127)313}314315/// Returns the plane design of a cell code as a radix design: base the rational integer `m`, of norm `m^2`, on the square lattice, no twist, digits the box residues `{x + y i : 0 <= x, y < m}`.316///317/// Bit `y m + x` of the code is the cell at row `y` and column `x` of the `m` by `m` tile, the column the real part and the row the imaginary part, so the level-`L` words scaled by `m^L` are exactly the filled cells of the level-`L` tile read as `(column, row)`.318pub fn tile(m: u64, code: u128) -> Radix {319    let base = Base::new(Ring::Gaussian, (m as i64, 0));320    let cells = (m * m) as usize;321    assert!(code >> cells == 0, "code out of range for the base");322    let digits: Vec<(i64, i64)> = (0..cells)323        .filter(|i| (code >> i) & 1 == 1)324        .map(|i| ((i as u64 % m) as i64, (i as u64 / m) as i64))325        .collect();326    let twists = vec![(1, 0); digits.len()];327    Radix::new(base, digits, twists)328}329330#[cfg(test)]331mod tests {332    use super::*;333334    fn mul(z: (f64, f64), w: (f64, f64)) -> (f64, f64) {335        (z.0 * w.0 - z.1 * w.1, z.0 * w.1 + z.1 * w.0)336    }337338    #[test]339    fn the_canonical_residues_are_a_complete_system() {340        let three = Base::new(Ring::Eisenstein, (3, 0));341        assert_eq!(342            three.residues(),343            vec![344                (0, 0),345                (1, 0),346                (1, 1),347                (0, 1),348                (-1, 0),349                (-1, -1),350                (0, -1),351                (2, 1),352                (1, 2)353            ]354        );355        assert_eq!(356            Base::new(Ring::Gaussian, (1, 1)).residues(),357            vec![(0, 0), (1, 0)]358        );359        assert_eq!(360            Base::new(Ring::Gaussian, (2, 1)).residues(),361            vec![(0, 0), (1, 0), (0, 1), (-1, 0), (0, -1)]362        );363        let bases = [364            Base::new(Ring::Gaussian, (2, 0)),365            Base::new(Ring::Gaussian, (1, 1)),366            Base::new(Ring::Gaussian, (2, 1)),367            Base::new(Ring::Gaussian, (3, 0)),368            Base::new(Ring::Eisenstein, (2, 0)),369            Base::new(Ring::Eisenstein, (2, 1)),370            Base::new(Ring::Eisenstein, (3, 0)),371            Base::new(Ring::Eisenstein, (3, 1)),372        ];373        for base in bases {374            let residues = base.residues();375            assert_eq!(residues.len(), base.norm() as usize);376            for (i, &z) in residues.iter().enumerate() {377                for &w in &residues[..i] {378                    assert!(!base.congruent(z, w), "{base:?} {z:?} {w:?}");379                }380            }381            for a in -9..=9 {382                for b in -9..=9 {383                    let hits = residues384                        .iter()385                        .filter(|&&r| base.congruent((a, b), r))386                        .count();387                    assert_eq!(hits, 1, "{base:?} {a} {b}");388                }389            }390            for map in base.group() {391                let seen: HashSet<usize> = map.iter().copied().collect();392                assert_eq!(seen.len(), residues.len(), "{base:?}");393            }394            assert_eq!(395                base.group().len(),396                base.ring().units() * if base.mirrored() { 2 } else { 1 }397            );398        }399        let multiset = Base::new(Ring::Gaussian, (1, 1)).group();400        let image: HashSet<&Vec<usize>> = multiset.iter().collect();401        assert_eq!(multiset.len(), 8);402        assert_eq!(image.len(), 1);403    }404405    #[test]406    #[should_panic(expected = "congruent modulo the base")]407    fn a_congruent_digit_pair_is_refused() {408        let base = Base::new(Ring::Eisenstein, (3, 0));409        assert!(base.congruent((2, 0), (-1, 0)));410        Radix::new(base, vec![(0, 0), (2, 0), (-1, 0)], vec![(1, 0); 3]);411    }412413    #[test]414    fn the_canonical_residues_build_a_design() {415        let bases = [416            Base::new(Ring::Gaussian, (1, 1)),417            Base::new(Ring::Gaussian, (2, 1)),418            Base::new(Ring::Eisenstein, (2, 0)),419            Base::new(Ring::Eisenstein, (3, 0)),420            Base::new(Ring::Eisenstein, (3, 1)),421        ];422        for base in bases {423            let digits = base.residues();424            let twists = vec![(1, 0); digits.len()];425            let design = Radix::new(base, digits, twists);426            assert_eq!(design.size(), base.norm() as usize);427            assert!(design.canonical());428        }429        assert_eq!(koch().size(), 4);430    }431432    #[test]433    fn the_terdragon_is_the_twist_the_l_system_reads() {434        let design = terdragon();435        assert_eq!(design.twists().to_vec(), vec![(1, 0), (0, 1), (1, 0)]);436        let level = 2;437        let mut word = b"F".to_vec();438        for _ in 0..level {439            let mut next = Vec::with_capacity(word.len() * 5);440            for &c in &word {441                if c == b'F' {442                    next.extend_from_slice(b"F+F-F");443                } else {444                    next.push(c);445                }446            }447            word = next;448        }449        let root = 3f64.sqrt();450        let step = [(1.0, 0.0), (-0.5, root / 2.0), (-0.5, -root / 2.0)];451        let mut heading = 0usize;452        let mut at = (0.0f64, 0.0f64);453        let mut starts = Vec::new();454        for c in word {455            match c {456                b'F' => {457                    starts.push(at);458                    at = (at.0 + step[heading].0, at.1 + step[heading].1);459                }460                b'+' => heading = (heading + 1) % 3,461                _ => heading = (heading + 2) % 3,462            }463        }464        let den = at.0 * at.0 + at.1 * at.1;465        let want: Vec<(f64, f64)> = starts466            .into_iter()467            .map(|p| {468                (469                    (p.0 * at.0 + p.1 * at.1) / den,470                    (p.1 * at.0 - p.0 * at.1) / den,471                )472            })473            .collect();474        let got = design.plane(level);475        assert_eq!(got.len(), 9);476        assert_eq!(want.len(), 9);477        let worst = got478            .iter()479            .zip(want.iter())480            .map(|(a, b)| ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt())481            .fold(0.0f64, f64::max);482        assert!(worst < 1e-12, "deviation {worst}");483    }484485    #[test]486    fn the_twisted_koch_words_are_the_textbook_maps() {487        let koch = koch();488        assert_eq!(koch.size(), 4);489        assert!(!koch.canonical());490        let root = 3f64.sqrt();491        let maps: [((f64, f64), (f64, f64)); 4] = [492            ((1.0 / 3.0, 0.0), (0.0, 0.0)),493            ((1.0 / 6.0, root / 6.0), (1.0 / 3.0, 0.0)),494            ((1.0 / 6.0, -root / 6.0), (0.5, root / 6.0)),495            ((1.0 / 3.0, 0.0), (2.0 / 3.0, 0.0)),496        ];497        let mut truth = vec![(0.0f64, 0.0f64)];498        for level in 1..=2 {499            let mut next = Vec::new();500            for &(scale, shift) in &maps {501                for &p in &truth {502                    let t = mul(scale, p);503                    next.push((t.0 + shift.0, t.1 + shift.1));504                }505            }506            truth = next;507            let got = koch.plane(level);508            assert_eq!(got.len(), truth.len());509            assert_eq!(got.len(), 4usize.pow(level as u32));510            let worst = got511                .iter()512                .zip(truth.iter())513                .map(|(a, b)| ((a.0 - b.0).powi(2) + (a.1 - b.1).powi(2)).sqrt())514                .fold(0.0f64, f64::max);515            assert!(worst < 1e-12, "level {level} deviation {worst}");516        }517        assert_eq!(koch.distinct(2), 16);518    }519520    #[test]521    fn the_real_base_with_no_twist_is_the_plane_cell() {522        let carpet = tile(3, 0b111101111);523        assert_eq!(carpet.size(), 8);524        assert!((carpet.dimension() - 8f64.ln() / 3f64.ln()).abs() < 1e-12);525        assert_eq!(carpet.fill(2), 64);526        let got: HashSet<(i64, i64)> = carpet.words(2).into_iter().collect();527        let mut want = HashSet::new();528        for row in 0..9i64 {529            for col in 0..9i64 {530                let hole = |r: i64, c: i64| r == 1 && c == 1;531                if !hole(row / 3, col / 3) && !hole(row % 3, col % 3) {532                    want.insert((col, row));533                }534            }535        }536        assert_eq!(got.len(), 64);537        assert_eq!(got, want);538    }539540    #[test]541    fn the_fill_law_counts_every_word_and_the_twist_only_moves_it() {542        for design in [gasket(), twindragon(), terdragon(), flowsnake(), koch()] {543            for level in 0..=4 {544                assert_eq!(design.words(level).len() as u128, design.fill(level));545                assert!(design.distinct(level) <= design.words(level).len());546            }547            assert!(548                (design.dimension() * (design.base().norm() as f64).ln() / 2.0549                    - (design.size() as f64).ln())550                .abs()551                    < 1e-12552            );553        }554        let plain = twindragon();555        let turned = plain.clone().with_twists(&[0, 1]);556        assert_eq!(plain.words(4).len(), turned.words(4).len());557        assert_ne!(plain.words(4), turned.words(4));558        let glue = Radix::from_code(Base::new(Ring::Gaussian, (2, 0)), 3).with_twists(&[0, 2]);559        assert_eq!(glue.words(2), vec![(0, 0), (1, 0), (2, 0), (1, 0)]);560        assert_eq!(glue.fill(2), 4);561        assert_eq!(glue.distinct(2), 3);562        assert_eq!(gasket().code(), 7);563        assert_eq!(flowsnake().code(), 127);564        assert!((gasket().dimension() - 3f64.ln() / 2f64.ln()).abs() < 1e-12);565        assert!((twindragon().dimension() - 2.0).abs() < 1e-12);566        assert!((terdragon().dimension() - 2.0).abs() < 1e-12);567        assert!((flowsnake().dimension() - 2.0).abs() < 1e-12);568    }569}