cocycle.rs
18.7 kB · rust · 597 lines
1use crate::order;2use crate::series::Rep;3use crate::word::{render, CODES};45pub const UNIT: [u8; 4] = [1, 2, 4, 8];6pub const ROW: [u8; 2] = [3, 12];7pub const COLUMN: [u8; 2] = [5, 10];8pub const DIAGONAL: [u8; 2] = [6, 9];9pub const GASKET: [u8; 4] = [7, 11, 13, 14];10pub const FULL: u8 = 15;11pub const REACH: usize = 14;12pub const DRAWN: usize = 7;13pub const HORIZON: usize = 1 << 20;1415#[derive(Clone, Copy, PartialEq, Eq)]16pub enum Form {17 UnitDomino,18 UnitDiagonal,19 RowColumn,20 DominoDiagonal,21 DominoFull,22 Empty,23 All,24}2526pub struct Family {27 pub name: &'static str,28 pub rule: &'static str,29 pub form: Form,30 pub pairs: Vec<(u8, u8)>,31}3233fn gcd(a: u64, b: u64) -> u64 {34 let (mut a, mut b) = (a, b);35 while b != 0 {36 let t = a % b;37 a = b;38 b = t;39 }40 if a == 0 {41 142 } else {43 a44 }45}4647fn ratio(num: u64, den: u64) -> String {48 let g = gcd(num, den);49 if den / g == 1 {50 return format!("{}", num / g);51 }52 format!("{}/{}", num / g, den / g)53}5455pub fn value(rep: &Rep, word: &[u8]) -> i128 {56 let mut state: Vec<i128> = rep.lambda.iter().map(|entry| *entry as i128).collect();57 for code in word {58 let matrix = &rep.matrices[code];59 let mut next = vec![0i128; state.len()];60 for (i, weight) in state.iter().enumerate() {61 if *weight == 0 {62 continue;63 }64 for (j, slot) in next.iter_mut().enumerate() {65 *slot += weight * matrix[i][j] as i128;66 }67 }68 state = next;69 }70 state71 .iter()72 .zip(rep.gamma.iter())73 .map(|(a, b)| a * *b as i128)74 .sum()75}7677pub fn exponent(form: Form, word: &[u8], marked: u8) -> u32 {78 match form {79 Form::UnitDomino => {80 let mut last = 0usize;81 let mut count = 0usize;82 for (i, code) in word.iter().enumerate() {83 if *code == marked {84 last = i + 1;85 count += 1;86 }87 }88 (last - count) as u3289 }90 Form::UnitDiagonal => word.iter().filter(|code| **code == marked).count() as u32,91 Form::RowColumn => {92 let last = *word.last().expect("a word has a last letter");93 let run = word.iter().rev().take_while(|code| **code == last).count();94 (word.len() - run) as u3295 }96 Form::DominoDiagonal => {97 let mut last = 0usize;98 for (i, code) in word.iter().enumerate() {99 if *code == marked {100 last = i + 1;101 }102 }103 last as u32104 }105 Form::DominoFull => {106 let count = word.iter().filter(|code| **code == marked).count();107 let run = word108 .iter()109 .rev()110 .take_while(|code| **code == marked)111 .count();112 (count - run) as u32113 }114 Form::Empty => 0,115 Form::All => word.len() as u32,116 }117}118119fn pair_word(a: u8, b: u8, mask: usize, length: usize) -> Vec<u8> {120 (0..length)121 .map(|i| if (mask >> i) & 1 == 1 { b } else { a })122 .collect()123}124125fn sweep(rep: &Rep, form: Form, a: u8, b: u8, reach: usize) -> (usize, usize) {126 let mut checked = 0usize;127 let mut bad = 0usize;128 for length in 1..=reach {129 for mask in 0..(1usize << length) {130 let word = pair_word(a, b, mask, length);131 let want = 1i128 << exponent(form, &word, b);132 if value(rep, &word) != want {133 bad += 1;134 }135 checked += 1;136 }137 }138 (checked, bad)139}140141fn drawn(form: Form, a: u8, b: u8, reach: usize) -> (usize, usize) {142 let mut checked = 0usize;143 let mut bad = 0usize;144 for length in 1..=reach {145 for mask in 0..(1usize << length) {146 let word = pair_word(a, b, mask, length);147 let want = 1u64 << exponent(form, &word, b);148 if render(&word).components() != want {149 bad += 1;150 }151 checked += 1;152 }153 }154 (checked, bad)155}156157pub fn families() -> Vec<Family> {158 let dominoes: Vec<u8> = ROW.iter().chain(COLUMN.iter()).copied().collect();159 let mut out: Vec<Family> = Vec::new();160 let mut pairs: Vec<(u8, u8)> = Vec::new();161 for &domino in dominoes.iter() {162 for &unit in UNIT.iter() {163 pairs.push((domino, unit));164 }165 }166 out.push(Family {167 name: "unit and domino",168 rule: "2^(k - m), k the last unit place, m the unit count",169 form: Form::UnitDomino,170 pairs,171 });172 let mut pairs: Vec<(u8, u8)> = Vec::new();173 for &unit in UNIT.iter() {174 for &diagonal in DIAGONAL.iter() {175 pairs.push((unit, diagonal));176 }177 }178 out.push(Family {179 name: "unit and diagonal",180 rule: "2^(diagonal count)",181 form: Form::UnitDiagonal,182 pairs,183 });184 let mut pairs: Vec<(u8, u8)> = Vec::new();185 for &row in ROW.iter() {186 for &column in COLUMN.iter() {187 pairs.push((row, column));188 }189 }190 out.push(Family {191 name: "crossed dominoes",192 rule: "2^(L - r), r the terminal run",193 form: Form::RowColumn,194 pairs,195 });196 let mut pairs: Vec<(u8, u8)> = Vec::new();197 for &domino in dominoes.iter() {198 for &diagonal in DIAGONAL.iter() {199 pairs.push((domino, diagonal));200 }201 }202 out.push(Family {203 name: "domino and diagonal",204 rule: "2^k, k the last diagonal place",205 form: Form::DominoDiagonal,206 pairs,207 });208 let pairs: Vec<(u8, u8)> = dominoes.iter().map(|domino| (*domino, FULL)).collect();209 out.push(Family {210 name: "domino and full",211 rule: "2^(n - j), n the full count, j the terminal full run",212 form: Form::DominoFull,213 pairs,214 });215 let pairs: Vec<(u8, u8)> = GASKET.iter().map(|code| (*code, FULL)).collect();216 out.push(Family {217 name: "gasket and full",218 rule: "1",219 form: Form::Empty,220 pairs,221 });222 let mut pairs: Vec<(u8, u8)> = Vec::new();223 for class in [224 UNIT.to_vec(),225 ROW.to_vec(),226 COLUMN.to_vec(),227 GASKET.to_vec(),228 ]229 .iter()230 {231 for i in 0..class.len() {232 for j in i + 1..class.len() {233 pairs.push((class[i], class[j]));234 }235 }236 }237 out.push(Family {238 name: "inside one flat class",239 rule: "1",240 form: Form::Empty,241 pairs,242 });243 out.push(Family {244 name: "inside the diagonal class",245 rule: "2^L",246 form: Form::All,247 pairs: vec![(DIAGONAL[0], DIAGONAL[1])],248 });249 out250}251252fn constants(rep: &Rep) {253 println!("CONSTANT WORDS");254 let mut doubling: Vec<u8> = Vec::new();255 let mut bad = 0usize;256 for &code in CODES.iter() {257 let one = value(rep, &[code]);258 for length in 1..8usize {259 let word = vec![code; length];260 if value(rep, &word) != one.pow(length as u32) {261 bad += 1;262 }263 }264 if one == 2 {265 doubling.push(code);266 }267 }268 println!("comp(c^L) = comp(c)^L for all 15 codes at L = 1..7, mismatches {bad}");269 println!(270 "comp(c) = 2 exactly for the codes {doubling:?} and 1 for the other {}",271 15 - doubling.len()272 );273 println!("so the constant-word rate is log 2 on the diagonal class and 0 elsewhere, and the frequency functional is Phi(f) = (f_6 + f_9) log 2");274 assert_eq!(bad, 0, "a constant word is a power");275 assert_eq!(doubling, vec![6, 9], "only the diagonal class doubles");276}277278fn forms(rep: &Rep) {279 println!();280 println!("CLOSED FORMS");281 for family in families().iter() {282 let mut checked = 0usize;283 let mut bad = 0usize;284 let mut seen = 0usize;285 let mut wrong = 0usize;286 for (a, b) in family.pairs.iter() {287 let (one, two) = sweep(rep, family.form, *a, *b, REACH);288 checked += one;289 bad += two;290 let (three, four) = drawn(family.form, *a, *b, DRAWN);291 seen += three;292 wrong += four;293 }294 println!(295 "{}: {} pairs, comp = {}, {checked} words to L = {REACH}, mismatches {bad}, {seen} words drawn to L = {DRAWN}, mismatches {wrong}",296 family.name,297 family.pairs.len(),298 family.rule299 );300 assert_eq!(301 bad, 0,302 "the closed form is exact against the representation"303 );304 assert_eq!(wrong, 0, "the closed form is exact against the drawn word");305 }306 let named: usize = families().iter().map(|family| family.pairs.len()).sum();307 let all = CODES.len() * (CODES.len() - 1) / 2;308 println!("{named} of the {all} letter pairs carry a closed form, 9 of the 15 pairs of distinct classes and all 6 pairs inside one class; the {} pairs left open are the gasket class against the unit, domino and diagonal classes and the full tile against the unit and diagonal classes", all - named);309}310311fn cut_law(rep: &Rep) {312 println!();313 println!("THE ZERO-CONTACT CUT");314 let mut across = [0u64; 16];315 let mut down = [0u64; 16];316 for &code in CODES.iter() {317 let (rows, columns) = render(&[code]).contacts();318 across[code as usize] = rows;319 down[code as usize] = columns;320 }321 let text = |table: &[u64; 16]| {322 CODES323 .iter()324 .map(|code| format!("{code}:{}", table[*code as usize]))325 .collect::<Vec<_>>()326 .join(" ")327 };328 println!("h {}", text(&across));329 println!("v {}", text(&down));330 let mut total = 0usize;331 let mut applies = 0usize;332 let mut bad = 0usize;333 for length in 1..5usize {334 for word in order::words(&CODES, length) {335 total += 1;336 let mut cut: Option<usize> = None;337 for start in 0..length {338 let rows: u64 = word[start..]339 .iter()340 .map(|code| across[*code as usize])341 .product();342 let columns: u64 = word[start..]343 .iter()344 .map(|code| down[*code as usize])345 .product();346 if rows == 0 && columns == 0 {347 cut = Some(start);348 }349 }350 if let Some(start) = cut {351 applies += 1;352 let ahead: i128 = word[..start]353 .iter()354 .map(|code| code.count_ones() as i128)355 .product();356 if value(rep, &word) != ahead * value(rep, &word[start..]) {357 bad += 1;358 }359 }360 }361 }362 println!("comp = fill(prefix) * comp(suffix) at the last zero-contact suffix: {applies} of {total} words of length 1 to 4 admit the cut, mismatches {bad}");363 assert_eq!(bad, 0, "the cut law is exact where it applies");364}365366pub fn morse(length: usize) -> Vec<u8> {367 (0..length).map(|n| (n.count_ones() & 1) as u8).collect()368}369370fn doubling(length: usize) -> Vec<u8> {371 let mut word: Vec<u8> = vec![1];372 while word.len() < length {373 let mut next: Vec<u8> = Vec::with_capacity(2 * word.len());374 for letter in word.iter() {375 next.push(1);376 next.push(1 - letter);377 }378 word = next;379 }380 word.truncate(length);381 word382}383384fn word_structure() {385 println!();386 println!("THE THUE-MORSE WORD");387 let bits = morse(HORIZON);388 let mut ones = 0usize;389 let mut uneven = 0usize;390 for (i, bit) in bits.iter().enumerate() {391 ones += *bit as usize;392 let length = i + 1;393 if length % 2 == 0 && ones * 2 != length {394 uneven += 1;395 }396 }397 println!("the prefix of every even length L <= {HORIZON} carries exactly L/2 of each letter, exceptions {uneven}");398 let mut longest = 0usize;399 let mut run = 0usize;400 let mut zeros = [0usize; 8];401 let mut zero_run = 0usize;402 for (i, bit) in bits.iter().enumerate() {403 if i > 0 && bits[i - 1] == *bit {404 run += 1;405 } else {406 run = 1;407 }408 zero_run = if *bit == 0 { zero_run + 1 } else { 0 };409 longest = longest.max(run);410 zeros[zero_run.min(7)] += 1;411 }412 println!("the longest terminal run over all L <= {HORIZON} is {longest}, either letter");413 println!(414 "the terminal run of the first letter takes the value 0 on {} prefixes, 1 on {}, 2 on {}, and never more",415 zeros[0], zeros[1], zeros[2]416 );417 let mut complete = [0usize; 8];418 let mut run = 1usize;419 for i in 1..bits.len() {420 if bits[i] == bits[i - 1] {421 run += 1;422 } else {423 complete[run.min(7)] += 1;424 run = 1;425 }426 }427 println!(428 "the first {HORIZON} letters hold {} complete runs of length 1 and {} of length 2, none longer, plus one unfinished run of length {run} at the cut",429 complete[1], complete[2]430 );431 let pd = doubling(bits.len());432 let mut bad = 0usize;433 for i in 0..bits.len() - 1 {434 if (bits[i] ^ bits[i + 1]) != pd[i] {435 bad += 1;436 }437 }438 println!("the run-boundary word t_n xor t_(n+1) is the period-doubling word on the first {} terms, mismatches {bad}", bits.len() - 1);439 assert_eq!(440 uneven, 0,441 "the letter counts are exactly equal at even length"442 );443 assert_eq!(longest, 2, "no three equal letters run together");444 assert_eq!(bad, 0, "the boundary word is the period-doubling word");445}446447struct Reading {448 pair: (u8, u8),449 num: u64,450 den: u64,451 rate: &'static str,452 prediction: &'static str,453 fill: &'static str,454}455456fn morse_rates(rep: &Rep) {457 println!();458 println!("THE THUE-MORSE EXPONENT");459 let bits = morse(HORIZON);460 let readings = [461 Reading {462 pair: (3, 1),463 num: 1,464 den: 2,465 rate: "(1/2) log 2",466 prediction: "0",467 fill: "(1/2) log 2",468 },469 Reading {470 pair: (1, 6),471 num: 1,472 den: 2,473 rate: "(1/2) log 2",474 prediction: "(1/2) log 2",475 fill: "(1/2) log 2",476 },477 Reading {478 pair: (3, 5),479 num: 1,480 den: 1,481 rate: "log 2",482 prediction: "0",483 fill: "log 2",484 },485 Reading {486 pair: (3, 6),487 num: 1,488 den: 1,489 rate: "log 2",490 prediction: "(1/2) log 2",491 fill: "log 2",492 },493 Reading {494 pair: (3, 15),495 num: 1,496 den: 2,497 rate: "(1/2) log 2",498 prediction: "0",499 fill: "(3/2) log 2",500 },501 Reading {502 pair: (7, 15),503 num: 0,504 den: 1,505 rate: "0",506 prediction: "0",507 fill: "(log 3 + log 4)/2",508 },509 ];510 for (family, reading) in families().iter().take(6).zip(readings.iter()) {511 let (a, b) = reading.pair;512 assert!(513 family.pairs.contains(&(a, b)),514 "the pick sits in its family"515 );516 let mut rates: Vec<String> = Vec::new();517 for swap in 0..2usize {518 let word: Vec<u8> = bits519 .iter()520 .map(|bit| if (*bit as usize) == swap { b } else { a })521 .collect();522 let mut short = 0usize;523 for length in 1..=REACH {524 let want = 1i128 << exponent(family.form, &word[..length], b);525 if value(rep, &word[..length]) != want {526 short += 1;527 }528 }529 assert_eq!(short, 0, "the closed form reads the prefix");530 let power = exponent(family.form, &word, b) as u64;531 let gap = (power * reading.den) as i64 - (reading.num * HORIZON as u64) as i64;532 assert!(533 gap.unsigned_abs() <= 2 * reading.den,534 "the prefix rate sits within 2/L of the limit"535 );536 rates.push(ratio(power, HORIZON as u64));537 }538 println!(539 "{}: pair ({a},{b}), rate {}, prediction {}, fill rate {}, prefix rate at L = {HORIZON} is {} and {} under the two letter readings",540 family.name, reading.rate, reading.prediction, reading.fill, rates[0], rates[1]541 );542 }543 println!("the rate is the fill rate on every family above except domino and full, where it lies strictly between the per-letter value 0 and the fill rate, and gasket and full, where the word is connected");544}545546fn boundary(rep: &Rep) {547 println!();548 println!("THE SIMPLEX BOUNDARY");549 let flat = 3u8;550 let mark = 6u8;551 let squares = |i: usize| {552 let root = (i as f64).sqrt() as usize;553 (root.saturating_sub(1)..root + 2).any(|r| r * r == i)554 };555 let mut bad = 0usize;556 for length in 1..=REACH {557 let powers: Vec<u8> = (1..=length)558 .map(|i| if i.is_power_of_two() { mark } else { flat })559 .collect();560 let boxes: Vec<u8> = (1..=length)561 .map(|i| if squares(i) { mark } else { flat })562 .collect();563 let flats = vec![flat; length];564 for word in [powers, boxes, flats] {565 let want = 1i128 << exponent(Form::DominoDiagonal, &word, mark);566 if value(rep, &word) != want {567 bad += 1;568 }569 }570 }571 println!("three words over the pair (3,6) whose marked letter has frequency 0, closed form checked to L = {REACH}, mismatches {bad}");572 let mut low: Vec<String> = Vec::new();573 for k in 2..14usize {574 let length = (1usize << (k + 1)) - 1;575 low.push(ratio(1u64 << k, length as u64));576 }577 println!("the marked letter at the powers of 2: the prefix rate is 1 at every L = 2^k and {} at L = 2^(k+1) - 1, so the rate has no limit, upper 1 and lower 1/2", low.join(", "));578 let mut squares_low: Vec<String> = Vec::new();579 for n in [2usize, 4, 8, 16, 32, 64] {580 squares_low.push(ratio((n * n) as u64, (n * n + 2 * n) as u64));581 }582 println!("the marked letter at the squares: the prefix rate is 1 at every L = n^2 and {} at L = (n+1)^2 - 1, so the rate is 1", squares_low.join(", "));583 println!("the constant word 3^L has rate 0, is periodic and hence minimal, and carries the same letter frequencies as both");584 println!("so at a frequency vector on the boundary of the simplex the exponent is 0, is 1, and fails to exist, and the closed forms give a rate only where both letters have positive frequency");585 assert_eq!(bad, 0, "the closed form covers the three boundary words");586}587588pub fn study(rep: &Rep) {589 println!();590 println!("THE COMPONENT COCYCLE");591 constants(rep);592 forms(rep);593 cut_law(rep);594 word_structure();595 morse_rates(rep);596 boundary(rep);597}