gasket.rs

1use crate::rules::{output, single_seed, Diagram};2use mrlycore::tensor::Tensor;34const CANDIDATES: [usize; 4] = [7, 11, 13, 14];5const DEPTH: usize = 8;67fn tile(code: usize, level: usize) -> Vec<Vec<u8>> {8    let mut side = 1usize;9    let mut cells = vec![vec![1u8]];10    for _ in 0..level {11        let next = side * 2;12        let mut out = vec![vec![0u8; next]; next];13        for row in 0..side {14            for column in 0..side {15                if cells[row][column] == 1 {16                    for i in 0..4usize {17                        if (code >> i) & 1 == 1 {18                            out[2 * row + i / 2][2 * column + i % 2] = 1;19                        }20                    }21                }22            }23        }24        side = next;25        cells = out;26    }27    cells28}2930fn crate_tile(code: usize, level: usize) -> Vec<Vec<u8>> {31    let seed = Tensor::of(32        (0..4).map(|i| ((code >> i) & 1) as u8).collect(),33        vec![2, 2],34    );35    let grid = seed.fractal(level);36    let side = 1usize << level;37    (0..side)38        .map(|row| (0..side).map(|column| grid.get(&[row, column])).collect())39        .collect()40}4142fn matches(diagram: &Diagram, code: usize, reading: &str) -> bool {43    for level in 1..=DEPTH {44        let grid = tile(code, level);45        let side = 1usize << level;46        for t in 0..side {47            for j in 0..side {48                let offset = match reading {49                    "right" => j as i64,50                    "left" => j as i64 - (side as i64 - 1),51                    _ => 2 * j as i64 - t as i64,52                };53                if diagram.signed(t, offset) != grid[t][j] {54                    return false;55                }56            }57        }58    }59    true60}6162fn hits(rule: usize, reading: &str) -> Vec<usize> {63    let diagram = single_seed(rule, 1 << DEPTH);64    CANDIDATES65        .into_iter()66        .filter(|&code| matches(&diagram, code, reading))67        .collect()68}6970fn polynomials(rows: usize) -> Vec<Vec<u8>> {71    let width = 2 * rows + 1;72    let mut out = Vec::with_capacity(rows);73    let mut row = vec![0u8; width];74    row[0] = 1;75    for _ in 0..rows {76        out.push(row.clone());77        let mut next = vec![0u8; width];78        for i in 0..width {79            next[i] = row[i];80            if i >= 1 {81                next[i] ^= row[i - 1];82            }83            if i >= 2 {84                next[i] ^= row[i - 2];85            }86        }87        row = next;88    }89    out90}9192fn population_totals(depth: usize) -> (Vec<u64>, Vec<u64>) {93    let rows = 1usize << depth;94    let width = 2 * rows + 1;95    let mut row = vec![0u8; width];96    row[0] = 1;97    let mut totals = Vec::with_capacity(depth + 1);98    let mut pairs = Vec::with_capacity(depth + 1);99    let mut running = 0u64;100    let mut adjacent = 0u64;101    let mut next_mark = 1usize;102    for t in 0..rows {103        running += row.iter().filter(|&&c| c == 1).count() as u64;104        adjacent += row.windows(2).filter(|w| w[0] == 1 && w[1] == 1).count() as u64;105        if t + 1 == next_mark {106            totals.push(running);107            pairs.push(adjacent);108            next_mark *= 2;109        }110        let mut next = vec![0u8; width];111        for i in 0..width {112            next[i] = row[i];113            if i >= 1 {114                next[i] ^= row[i - 1];115            }116            if i >= 2 {117                next[i] ^= row[i - 2];118            }119        }120        row = next;121    }122    (totals, pairs)123}124125pub fn report() {126    println!("THE GASKET IDENTITY");127    for level in 1..=DEPTH {128        for code in CANDIDATES {129            assert_eq!(130                tile(code, level),131                crate_tile(code, level),132                "the study renderer and Tensor::fractal disagree at code {code} level {level}"133            );134        }135    }136    println!("both renderers agree on codes {CANDIDATES:?} at levels 1..{DEPTH}");137    for (rule, reading, want) in [138        (60usize, "right", 13usize),139        (102, "left", 14),140        (90, "shear", 13),141    ] {142        let found = hits(rule, reading);143        println!("rule {rule} read {reading}: matching codes {found:?}");144        assert_eq!(145            found,146            vec![want],147            "rule {rule} does not match exactly bang dim 2, code {want}"148        );149    }150    println!("rule 60 is bang dim 2, code 13; rule 102 is bang dim 2, code 14; rule 90 sheared by j = (t+i)/2 is bang dim 2, code 13, cell for cell to level {DEPTH}");151    println!("RULE 150");152    let rows = polynomials(129);153    let diagram = single_seed(150, 128);154    for t in 0..=128usize {155        for j in 0..=t {156            assert_eq!(157                diagram.signed(t, j as i64 - t as i64),158                rows[t][j],159                "the GF(2) row polynomial and the evolved diagram disagree at ({t},{j})"160            );161        }162    }163    println!("the row polynomial (1+x+x^2)^t matches the evolved diagram on rows 0..128");164    for t in 0..64usize {165        for j in 0..rows[2 * t].len() {166            let want = if j % 2 == 0 { rows[t][j / 2] } else { 0 };167            assert_eq!(168                rows[2 * t][j],169                want,170                "row {} is not row {t} spread by two",171                2 * t172            );173        }174        for j in 0..rows[2 * t + 1].len() {175            let mut want = rows[2 * t][j];176            if j >= 1 {177                want ^= rows[2 * t][j - 1];178            }179            if j >= 2 {180                want ^= rows[2 * t][j - 2];181            }182            assert_eq!(183                rows[2 * t + 1][j],184                want,185                "row {} is not row {} xor its two unit shifts",186                2 * t + 1,187                2 * t188            );189        }190    }191    println!("(1+x+x^2)^(2n) = (1+x^2+x^4)^n on rows 0..128: row 2t is row t spread by two, row 2t+1 is row 2t xor its two unit shifts");192    let populations: Vec<usize> = (0..=64)193        .map(|t| rows[t].iter().filter(|&&c| c == 1).count())194        .collect();195    println!("row populations t = 0..64");196    println!(197        "{}",198        populations199            .iter()200            .map(|p| p.to_string())201            .collect::<Vec<_>>()202            .join(",")203    );204    println!("this is OEIS A071053, the number of ON cells at generation t of rule 150 from a single ON cell");205    let (totals, pairs) = population_totals(12);206    let root = 5f64.sqrt();207    let constant = (5.0 + 3.0 * root) / 10.0;208    println!("k P(k) log2(P)/k P/(1+sqrt5)^k");209    for (k, total) in totals.iter().enumerate().skip(1) {210        println!(211            "{} {} {:.9} {:.9}",212            k,213            total,214            (*total as f64).log2() / k as f64,215            *total as f64 / (1.0 + root).powi(k as i32)216        );217    }218    println!("k P(k) B(k)");219    for k in 0..totals.len() {220        println!("{k} {} {}", totals[k], pairs[k]);221    }222    for k in 0..totals.len() - 1 {223        assert_eq!(224            totals[k + 1],225            4 * totals[k] - 2 * pairs[k],226            "P(k+1) = 4 P(k) - 2 B(k) breaks at k = {k}"227        );228        assert_eq!(229            pairs[k + 1],230            2 * totals[k] - 2 * pairs[k],231            "B(k+1) = 2 P(k) - 2 B(k) breaks at k = {k}"232        );233    }234    println!("the two-term system P(k+1) = 4 P(k) - 2 B(k), B(k+1) = 2 P(k) - 2 B(k) holds at k = 0..{}, matrix [[4,-2],[2,-2]] of trace 2 and determinant -4", totals.len() - 2);235    for k in 2..totals.len() {236        assert_eq!(237            totals[k],238            2 * totals[k - 1] + 4 * totals[k - 2],239            "the cumulative population breaks P(k+1) = 2 P(k) + 4 P(k-1) at k = {k}"240        );241    }242    let mut exact = vec![1u64, 4];243    while exact.len() <= totals.len() {244        let n = exact.len();245        exact.push(2 * exact[n - 1] + 4 * exact[n - 2]);246    }247    assert_eq!(248        totals[..],249        exact[..totals.len()],250        "the cumulative population leaves the closed form"251    );252    for (k, total) in totals.iter().enumerate() {253        let e = k as i32;254        let closed = constant * (1.0 + root).powi(e) + (1.0 - constant) * (1.0 - root).powi(e);255        assert!(256            (closed - *total as f64).abs() < 1e-3,257            "the closed form misses P({e}) = {total}"258        );259    }260    let mut fib = vec![0u64, 1];261    while fib.len() < totals.len() + 3 {262        let n = fib.len();263        fib.push(fib[n - 1] + fib[n - 2]);264    }265    for (k, total) in totals.iter().enumerate() {266        assert_eq!(*total, (1u64 << k) * fib[k + 2], "P({k}) is not 2^k F(k+2)");267    }268    println!("P(k) = 2 P(k-1) + 4 P(k-2) with P(0) = 1, P(1) = 4, so P(k) = c (1+sqrt5)^k + (1-c) (1-sqrt5)^k with c = (5+3 sqrt5)/10 = {constant:.12}");269    println!("equivalently P(k) = 2^k F(k+2), asserted at k = 0..12");270    println!(271        "the exponent is exactly log2(1+sqrt5) = {:.7}",272        (1.0 + root).log2()273    );274    println!("P(k) over rows 0..2^k - 1 is OEIS A087206, whose %N carries the recurrence and whose %F carries 2^n Fibonacci(n+2) and the (1 +- sqrt5) form; A071053 %F states Sum_{{k = 0..2^n - 1}} a(k) = A087206(n). The study re-derives them.");275    assert_eq!(output(150, 1, 1, 1), 1, "rule 150 is not the xor rule");276}