unsolved.rs
33.0 kB · rust · 976 lines
1use crate::cocycle::{morse, value, COLUMN, DIAGONAL, FULL, GASKET, ROW, UNIT};2use crate::series::{Frac, Rep};3use crate::word::{render, CODES};4use mrlycore::rng::Rng;56pub const DEEP: usize = 13;7pub const SEEN: usize = 7;8pub const RANGE: usize = 1 << 14;9pub const FAR: usize = 1 << 15;10pub const EXACT: usize = 60;11const SEED: u64 = 1618033988;1213#[derive(Clone, Copy, PartialEq, Eq)]14pub enum Shape {15 ZeroContact,16 GasketDomino,17}1819pub struct Open {20 pub name: &'static str,21 pub rule: &'static str,22 pub shape: Shape,23 pub pairs: Vec<(u8, u8)>,24}2526fn fill_of(word: &[u8]) -> i128 {27 word.iter().map(|code| code.count_ones() as i128).product()28}2930pub fn count(shape: Shape, word: &[u8], heavy: u8, light: u8) -> i128 {31 let last = word.iter().rposition(|code| *code == light);32 match (shape, last) {33 (_, None) => 1,34 (Shape::ZeroContact, Some(at)) => fill_of(&word[..=at]),35 (Shape::GasketDomino, Some(at)) => {36 let mut total = 1i128;37 for (i, code) in word[..=at].iter().enumerate() {38 if *code == heavy {39 total += fill_of(&word[..i]);40 }41 }42 total43 }44 }45}4647fn pair_word(heavy: u8, light: u8, mask: usize, length: usize) -> Vec<u8> {48 (0..length)49 .map(|i| if (mask >> i) & 1 == 1 { light } else { heavy })50 .collect()51}5253fn cross(left: &[u8], right: &[u8]) -> Vec<(u8, u8)> {54 let mut out: Vec<(u8, u8)> = Vec::new();55 for &a in left.iter() {56 for &b in right.iter() {57 out.push((a, b));58 }59 }60 out61}6263pub fn open_families() -> Vec<Open> {64 let dominoes: Vec<u8> = ROW.iter().chain(COLUMN.iter()).copied().collect();65 vec![66 Open {67 name: "gasket and unit",68 rule: "3^(g - j), g the gasket count, j the terminal gasket run",69 shape: Shape::ZeroContact,70 pairs: cross(&GASKET, &UNIT),71 },72 Open {73 name: "gasket and diagonal",74 rule: "2^d 3^(g - j), d the diagonal count",75 shape: Shape::ZeroContact,76 pairs: cross(&GASKET, &DIAGONAL),77 },78 Open {79 name: "full and unit",80 rule: "4^(F - j), F the full count, j the terminal full run",81 shape: Shape::ZeroContact,82 pairs: cross(&[FULL], &UNIT),83 },84 Open {85 name: "full and diagonal",86 rule: "2^d 4^(F - j)",87 shape: Shape::ZeroContact,88 pairs: cross(&[FULL], &DIAGONAL),89 },90 Open {91 name: "gasket and domino",92 rule:93 "1 + sum over the gasket places i <= m of fill(w_1..i-1), m the last domino place",94 shape: Shape::GasketDomino,95 pairs: cross(&GASKET, &dominoes),96 },97 ]98}99100fn forms(rep: &Rep) {101 println!("CLOSED FORMS ON THE 46");102 let mut pairs = 0usize;103 let mut words = 0usize;104 let mut drawn = 0usize;105 for family in open_families().iter() {106 let mut checked = 0usize;107 let mut bad = 0usize;108 let mut seen = 0usize;109 let mut wrong = 0usize;110 for (heavy, light) in family.pairs.iter() {111 for length in 1..=DEEP {112 for mask in 0..(1usize << length) {113 let word = pair_word(*heavy, *light, mask, length);114 let want = count(family.shape, &word, *heavy, *light);115 if value(rep, &word) != want {116 bad += 1;117 }118 checked += 1;119 if length <= SEEN {120 if render(&word).components() as i128 != want {121 wrong += 1;122 }123 seen += 1;124 }125 }126 }127 }128 println!(129 "{}: {} pairs, comp = {}, {checked} words to L = {DEEP}, mismatches {bad}, {seen} words drawn to L = {SEEN}, mismatches {wrong}",130 family.name,131 family.pairs.len(),132 family.rule133 );134 assert_eq!(135 bad, 0,136 "the closed form is exact against the representation"137 );138 assert_eq!(wrong, 0, "the closed form is exact against the drawn cells");139 pairs += family.pairs.len();140 words += checked;141 drawn += seen;142 }143 println!("{pairs} pairs, {words} words against the representation and {drawn} against the drawn cells, mismatches 0; with the 59 above every one of the 105 letter pairs now carries an exact closed form");144 assert_eq!(pairs, 46, "the five families cover the 46 open pairs");145}146147fn class_of(code: u8) -> usize {148 if UNIT.contains(&code) {149 0150 } else if ROW.contains(&code) {151 1152 } else if COLUMN.contains(&code) {153 2154 } else if DIAGONAL.contains(&code) {155 3156 } else if GASKET.contains(&code) {157 4158 } else {159 5160 }161}162163fn class_name(which: usize) -> &'static str {164 [165 "unit",166 "row domino",167 "column domino",168 "diagonal",169 "gasket",170 "full",171 ][which]172}173174fn weight(code: u8) -> (i64, i64) {175 match code.count_ones() {176 1 => (0, 0),177 2 => (1, 0),178 3 => (0, 1),179 _ => (2, 0),180 }181}182183fn phi_weight(code: u8) -> (i64, i64) {184 if DIAGONAL.contains(&code) {185 (1, 0)186 } else {187 (0, 0)188 }189}190191fn exponent_weights(a: u8, b: u8) -> ((i64, i64), (i64, i64)) {192 let (ca, cb) = (class_of(a), class_of(b));193 let light = (0i64, 0i64);194 let two = (1i64, 0i64);195 let three = (0i64, 1i64);196 let four = (2i64, 0i64);197 let (low, high) = if ca <= cb { (ca, cb) } else { (cb, ca) };198 let (first, second) = match (low, high) {199 (0, 0) => (light, light),200 (0, 1) | (0, 2) => (light, two),201 (0, 3) => (light, two),202 (0, 4) => (light, three),203 (0, 5) => (light, four),204 (1, 1) | (2, 2) => (light, light),205 (1, 2) => (two, two),206 (1, 3) | (2, 3) => (two, two),207 (1, 4) | (2, 4) => (two, three),208 (1, 5) | (2, 5) => (light, two),209 (3, 3) => (two, two),210 (3, 4) => (two, three),211 (3, 5) => (two, four),212 (4, 4) => (light, light),213 _ => (light, light),214 };215 if ca <= cb {216 (first, second)217 } else {218 (second, first)219 }220}221222fn nats(pair: (i64, i64)) -> f64 {223 pair.0 as f64 * 2f64.ln() + pair.1 as f64 * 3f64.ln()224}225226fn ledger() {227 println!();228 println!("THE FREQUENCY LEDGER");229 let mut saturating = 0usize;230 let mut short: Vec<(usize, usize)> = Vec::new();231 let mut exact = 0usize;232 let mut refuted: Vec<(usize, usize)> = Vec::new();233 let mut between: Vec<(usize, usize)> = Vec::new();234 for i in 0..CODES.len() {235 for j in i + 1..CODES.len() {236 let (a, b) = (CODES[i], CODES[j]);237 let (wa, wb) = exponent_weights(a, b);238 let key = (class_of(a).min(class_of(b)), class_of(a).max(class_of(b)));239 if wa == weight(a) && wb == weight(b) {240 saturating += 1;241 } else {242 short.push(key);243 }244 if wa == phi_weight(a) && wb == phi_weight(b) {245 exact += 1;246 } else {247 refuted.push(key);248 }249 let middle = (nats(wa) + nats(wb)) / 2.0;250 let prediction = (nats(phi_weight(a)) + nats(phi_weight(b))) / 2.0;251 let ceiling = (nats(weight(a)) + nats(weight(b))) / 2.0;252 if middle > prediction + 1e-12 && middle < ceiling - 1e-12 {253 between.push(key);254 }255 }256 }257 let tally = |list: &[(usize, usize)]| {258 let mut seen: Vec<((usize, usize), usize)> = Vec::new();259 for key in list.iter() {260 match seen.iter_mut().find(|(at, _)| at == key) {261 Some((_, n)) => *n += 1,262 None => seen.push((*key, 1)),263 }264 }265 seen.iter()266 .map(|((low, high), n)| {267 if low == high {268 format!("{} {n}", class_name(*low))269 } else {270 format!("{} and {} {n}", class_name(*low), class_name(*high))271 }272 })273 .collect::<Vec<_>>()274 .join(", ")275 };276 println!(277 "at interior frequency the exponent saturates the fill ceiling on {saturating} of the 105 pairs and falls short on {}: {}",278 short.len(),279 tally(&short)280 );281 println!(282 "Phi(f) = (f_6 + f_9) log 2 is exact on {exact} pairs and refuted on {}: {}",283 refuted.len(),284 tally(&refuted)285 );286 println!(287 "at equal frequencies exactly {} pairs sit strictly between the prediction and the fill ceiling: {}",288 between.len(),289 tally(&between)290 );291 assert_eq!(saturating + short.len(), 105, "every pair is counted once");292 assert_eq!(exact + refuted.len(), 105, "every pair is counted once");293 assert_eq!(between.len(), 4, "domino against the full tile is alone");294}295296struct Track {297 two: Vec<i64>,298 three: Vec<i64>,299 gaskets: Vec<usize>,300 reached: Vec<usize>,301 last: Vec<usize>,302}303304impl Track {305 fn new(word: &[u8], gasket: u8) -> Track {306 let n = word.len();307 let mut two = vec![0i64; n + 1];308 let mut three = vec![0i64; n + 1];309 let mut gaskets: Vec<usize> = Vec::new();310 let mut reached = vec![0usize; n + 1];311 let mut last = vec![0usize; n + 1];312 for i in 1..=n {313 let heavy = word[i - 1] == gasket;314 two[i] = two[i - 1] + if heavy { 0 } else { 1 };315 three[i] = three[i - 1] + if heavy { 1 } else { 0 };316 if heavy {317 gaskets.push(i);318 last[i] = last[i - 1];319 } else {320 last[i] = i;321 }322 reached[i] = gaskets.len();323 }324 Track {325 two,326 three,327 gaskets,328 reached,329 last,330 }331 }332333 fn log2_fill(&self, at: usize, scale: f64) -> f64 {334 self.two[at] as f64 + self.three[at] as f64 * scale335 }336337 fn top(&self, length: usize) -> Option<usize> {338 let cut = self.last[length];339 if cut == 0 {340 return None;341 }342 let reach = self.reached[cut];343 if reach == 0 {344 None345 } else {346 Some(self.gaskets[reach - 1])347 }348 }349350 fn ratio(&self, length: usize, scale: f64) -> Option<f64> {351 let top = self.top(length)?;352 let base = self.log2_fill(top - 1, scale);353 let reach = self.reached[self.last[length]];354 let mut total = (-base).exp2();355 for index in (0..reach).rev() {356 let place = self.gaskets[index];357 let term = (self.log2_fill(place - 1, scale) - base).exp2();358 if term < 1e-25 {359 break;360 }361 total += term;362 }363 Some(total)364 }365366 fn log2_comp(&self, length: usize, scale: f64) -> f64 {367 match self.top(length) {368 None => 0.0,369 Some(top) => {370 self.log2_fill(top - 1, scale)371 + self372 .ratio(length, scale)373 .expect("a top place has a ratio")374 .log2()375 }376 }377 }378}379380fn gcd(a: u128, b: u128) -> u128 {381 let (mut a, mut b) = (a, b);382 while b != 0 {383 let t = a % b;384 a = b;385 b = t;386 }387 if a == 0 {388 1389 } else {390 a391 }392}393394fn exact_saturation(word: &[u8], heavy: u8, light: u8, length: usize) -> (u128, u128) {395 let comp = count(Shape::GasketDomino, &word[..length], heavy, light) as u128;396 let fill = fill_of(&word[..length]) as u128;397 let g = gcd(comp, fill);398 (comp / g, fill / g)399}400401fn zero_contact_log2(word: &[u8], light: u8, scale: f64) -> f64 {402 match word.iter().rposition(|code| *code == light) {403 None => 0.0,404 Some(at) => {405 let mut two = 0i64;406 let mut three = 0i64;407 let mut four = 0i64;408 for code in word[..=at].iter() {409 match code.count_ones() {410 2 => two += 1,411 3 => three += 1,412 4 => four += 1,413 _ => (),414 }415 }416 two as f64 + three as f64 * scale + four as f64 * 2.0417 }418 }419}420421fn morse_word(gasket: u8, domino: u8, swap: usize, length: usize) -> Vec<u8> {422 morse(length)423 .iter()424 .map(|bit| {425 if (*bit as usize) == swap {426 gasket427 } else {428 domino429 }430 })431 .collect()432}433434fn morse_value(rep: &Rep) {435 println!();436 println!("THE THUE-MORSE VALUE ON THE GASKET-DOMINO PAIRS");437 let scale = 3f64.log2();438 let target = (1.0 + scale) / 2.0;439 let bound = 108f64.ln() + 1.5f64.ln() / 2.0;440 let dominoes: Vec<u8> = ROW.iter().chain(COLUMN.iter()).copied().collect();441 let mut worst = 0.0f64;442 let mut short = 0usize;443 let mut empty = 0usize;444 for &gasket in GASKET.iter() {445 for &domino in dominoes.iter() {446 for swap in 0..2usize {447 let word = morse_word(gasket, domino, swap, RANGE);448 let track = Track::new(&word, gasket);449 for length in 1..=DEEP {450 if value(rep, &word[..length])451 != count(Shape::GasketDomino, &word[..length], gasket, domino)452 {453 short += 1;454 }455 }456 for length in 4..=RANGE {457 if track.top(length).is_none() {458 empty += 1;459 }460 let gap = (track.log2_comp(length, scale)461 - length as f64 * (1.0 + scale) / 2.0)462 .abs()463 * 2f64.ln();464 worst = worst.max(gap);465 }466 }467 }468 }469 assert_eq!(short, 0, "the closed form reads every Thue-Morse prefix");470 assert_eq!(empty, 0, "a last gasket place exists from L = 4");471 println!("all 16 gasket-domino pairs, both letter readings, closed form against the representation to L = {DEEP}, mismatches 0, and the prefixes of length 4 or more with no gasket place at or before the last domino place number {empty}");472 let word = morse_word(7, 3, 0, RANGE);473 let track = Track::new(&word, 7);474 let other = morse_word(7, 3, 1, RANGE);475 let flip = Track::new(&other, 7);476 let marks = [256usize, 1024, 4096, RANGE];477 let line = |run: &Track| {478 marks479 .iter()480 .map(|length| {481 format!(482 "{:.12} ({length})",483 run.log2_comp(*length, scale) / *length as f64484 )485 })486 .collect::<Vec<_>>()487 .join(", ")488 };489 println!(490 "prefix rate log2 comp / L over (3,7), gasket at t = 0: {}",491 line(&track)492 );493 println!(494 "prefix rate log2 comp / L over (3,7), gasket at t = 1: {}",495 line(&flip)496 );497 println!("the limit is (1/2) log 6, which is {target:.15} in log 2 units and 0.895879734614027... in nats, both floats");498 println!(499 "over 4 <= L <= {RANGE}, all 16 pairs and both readings, the largest |log comp - (L/2) log 6| is {worst:.6} nats against the certificate log 108 + (1/2) log(3/2) = {bound:.6} nats, itself below 4.885"500 );501 assert!(worst < bound, "the certificate holds at every length");502 let mut deficit = 0.0f64;503 let mut ties = 0usize;504 for (name, run, letters) in [505 ("gasket at t = 0", &track, &word),506 ("gasket at t = 1", &flip, &other),507 ] {508 let sat =509 |length: usize| (run.log2_comp(length, scale) - run.log2_fill(length, scale)).exp2();510 let mut floor = (1.0f64, 0usize);511 let mut roof = (0.0f64, 0usize);512 let mut past = (0.0f64, 0usize);513 let mut over = 0usize;514 for length in 1..=RANGE {515 let here = sat(length);516 if here < floor.0 {517 floor = (here, length);518 }519 if length >= 4 && here > roof.0 {520 roof = (here, length);521 }522 if length >= 5 && here > past.0 {523 past = (here, length);524 }525 if here > 5.0 / 12.0 {526 over += 1;527 }528 deficit = deficit.max(-here.log2());529 }530 let show = |at: usize| {531 if at > EXACT {532 return format!("{:.10} at L = {at}, a float", sat(at));533 }534 let (n, d) = exact_saturation(letters, 7, 3, at);535 format!("{n}/{d} at L = {at}")536 };537 println!(538 "saturation comp/fill, {name}: over 1 <= L <= {RANGE} the minimum is {:.10} at L = {}, the value at L = 4096 is {:.10}, both floats, against the proved floor 1/108 = 0.009259259",539 floor.0,540 floor.1,541 sat(4096)542 );543 println!(544 "saturation comp/fill, {name}: the largest value at L >= 4 is {} and at L >= 5 is {}, exact rationals, and the number of lengths in the whole sweep above 5/12 is {over}, all of them below L = 4, which the certificate excludes",545 show(roof.1),546 show(past.1)547 );548 }549 for length in 1..=RANGE {550 let here = (track.log2_comp(length, scale) - track.log2_fill(length, scale)).exp2();551 if (-here.log2() - deficit).abs() < 1e-9 {552 ties += 1;553 }554 }555 println!(556 "the largest fill deficit log2 fill - log2 comp over both readings is {deficit:.4}, a float, and it is attained on a tie set, {ties} lengths in the reading gasket at t = 0 alone, so no single length may be named for it; the proved ceiling is log2 108 = {:.4}",557 108f64.log2()558 );559 assert!(560 deficit < 108f64.log2(),561 "the deficit stays under its ceiling"562 );563 let mut rates: Vec<String> = Vec::new();564 for big in [7u8, 15] {565 for swap in 0..2usize {566 let word: Vec<u8> = morse(RANGE)567 .iter()568 .map(|bit| if (*bit as usize) == swap { big } else { 6 })569 .collect();570 let exponent = zero_contact_log2(&word, 6, scale);571 let first = if swap == 0 { big } else { 6 };572 let second = if swap == 0 { 6 } else { big };573 rates.push(format!(574 "({first},{second}) {:.12}",575 exponent / RANGE as f64576 ));577 }578 }579 println!(580 "the same word over four zero-contact pairs at L = {RANGE}, log 2 units, floats: {}, against (1/2) log2 6 = {target:.12} and (1/2) log2 8 = 1.5",581 rates.join(", ")582 );583}584585fn named_words() -> Vec<(String, Vec<u8>)> {586 let mut out: Vec<(String, Vec<u8>)> = Vec::new();587 out.push(("Thue-Morse (3,7)".into(), morse_word(7, 3, 0, 4096)));588 out.push(("Thue-Morse swapped".into(), morse_word(7, 3, 1, 4096)));589 out.push((590 "(3,7)^2000".into(),591 (0..4000).map(|i| if i % 2 == 0 { 3 } else { 7 }).collect(),592 ));593 out.push((594 "(7,3)^2000".into(),595 (0..4000).map(|i| if i % 2 == 0 { 7 } else { 3 }).collect(),596 ));597 let mut rng = Rng::new(SEED);598 out.push((599 format!("Bernoulli(1/2) at seed {SEED}"),600 (0..4000)601 .map(|_| if rng.boolean() { 7 } else { 3 })602 .collect(),603 ));604 let mut block: Vec<u8> = vec![3; 2000];605 block.extend(vec![7; 2000]);606 block.extend(vec![3; 2000]);607 out.push(("3^2000 7^2000 3^2000".into(), block));608 out609}610611fn sandwich() {612 println!();613 println!("THE SANDWICH");614 let scale = 3f64.log2();615 println!("T < comp <= 1 + (3/2) T at T = fill(w_1..i*-1), i* the last gasket place at or before the last domino place");616 for (name, word) in named_words().iter() {617 let track = Track::new(word, 7);618 let mut low = f64::MAX;619 let mut high = 0.0f64;620 let mut bad = 0usize;621 for length in 1..=word.len() {622 let ratio = match track.ratio(length, scale) {623 None => continue,624 Some(value) => value,625 };626 low = low.min(ratio);627 high = high.max(ratio);628 let top = track.top(length).expect("a ratio needs a top place");629 let size = track.log2_fill(top - 1, scale);630 if ratio > 1.5 + (-size).exp2() + 1e-12 || ratio <= 1.0 {631 bad += 1;632 }633 }634 println!("{name}: comp / T in [{low:.4}, {high:.4}], floats, violations {bad}");635 assert_eq!(bad, 0, "the sandwich holds at every length");636 }637}638639fn text(value: &Frac) -> String {640 let (num, den) = value.parts();641 if den == 1 {642 format!("{num}")643 } else {644 format!("{num}/{den}")645 }646}647648fn chart_step(state: &[Frac; 4], matrix: &[Vec<i64>], divisor: i64) -> [Frac; 4] {649 let mut out = [Frac::zero(); 4];650 for i in 0..4 {651 for j in 0..4 {652 out[j] = out[j].add(&state[i].mul(&Frac::int(matrix[i][j])));653 }654 }655 for slot in out.iter_mut() {656 *slot = slot.div(&Frac::int(divisor));657 }658 out659}660661fn image(point: (Frac, Frac), gasket: bool) -> (Frac, Frac) {662 let one = Frac::int(1);663 let (b, c) = point;664 if gasket {665 (666 one.add(&b).div(&Frac::int(3)),667 one.add(&c).div(&Frac::int(3)),668 )669 } else {670 (one.add(&b).div(&Frac::int(2)), Frac::zero())671 }672}673674fn inside(point: &(Frac, Frac)) -> bool {675 let (b, c) = point;676 let zero = Frac::zero();677 let one = Frac::int(1);678 let half = Frac::new(1, 2);679 !b.below(&zero) && !one.below(b) && !c.below(&zero) && !half.below(c) && !one.below(&b.add(c))680}681682fn cone(rep: &Rep) {683 println!();684 println!("THE INVARIANT CONE AND WHY IT IS NOT THE MECHANISM");685 println!("phi = (1,2,2,4)^T is a common RIGHT eigenvector, M_c phi = popcount(c) phi, so it normalises the ROW orbit lambda M_(c_1) ... M_(c_k) alone");686 let phi = [1i64, 2, 2, 4];687 let mut bad = 0usize;688 for &code in CODES.iter() {689 let matrix = &rep.matrices[&code];690 for i in 0..4 {691 let row: i64 = (0..4).map(|j| matrix[i][j] * phi[j]).sum();692 if row != code.count_ones() as i64 * phi[i] {693 bad += 1;694 }695 }696 }697 assert_eq!(bad, 0, "phi is a right eigenvector at every code");698 println!("checked on all 15 code matrices, mismatches {bad}; M_gasket gamma = M_full gamma = gamma, M_unit gamma = phi, M_diagonal gamma = 2 phi, and M_diagonal = 2 phi lambda");699 let word: [u8; 7] = [7, 3, 7, 7, 3, 3, 7];700 let mut state = [Frac::int(1), Frac::zero(), Frac::zero(), Frac::zero()];701 let mut point = (Frac::zero(), Frac::zero());702 let mut orbit: Vec<String> = Vec::new();703 let mut drift = 0usize;704 for code in word.iter() {705 state = chart_step(&state, &rep.matrices[code], code.count_ones() as i64);706 point = image(point, *code == 7);707 if state[1] != point.0 || state[2] != point.1 || state[3] != Frac::zero() {708 drift += 1;709 }710 orbit.push(format!("({},{})", text(&point.0), text(&point.1)));711 }712 assert_eq!(drift, 0, "the chart is the normalised row orbit");713 println!("in the chart n_4 = 0, (b,c) = (n_2,n_3) and comp/fill = 1 - b - c, the two maps are N_gasket(b,c) = ((1+b)/3, (1+c)/3) and N_domino(b,c) = ((1+b)/2, 0)");714 println!("along the word 7,3,7,7,3,3,7 the chart reads {}, matching the raw matrices at every step, mismatches {drift}", orbit.join(" "));715 let vertices = [716 (Frac::zero(), Frac::zero()),717 (Frac::int(1), Frac::zero()),718 (Frac::zero(), Frac::new(1, 2)),719 (Frac::new(1, 2), Frac::new(1, 2)),720 ];721 let mut escapes = 0usize;722 for vertex in vertices.iter() {723 for gasket in [true, false] {724 if !inside(&image(vertex.clone(), gasket)) {725 escapes += 1;726 }727 }728 }729 assert_eq!(escapes, 0, "S is invariant under both maps");730 println!("S = {{0 <= b <= 1, 0 <= c <= 1/2, b + c <= 1}} is invariant under both maps, checked at all four vertices, escapes {escapes}");731 let mut images: Vec<String> = Vec::new();732 let mut widest = Frac::zero();733 for vertex in vertices.iter() {734 let mut at = vertex.clone();735 for gasket in [true, false, true] {736 at = image(at, gasket);737 }738 let sum = at.0.add(&at.1);739 if widest.below(&sum) {740 widest = sum;741 }742 images.push(format!("({},{})", text(&at.0), text(&at.1)));743 }744 println!("gasket-domino-gasket sends the vertices to {} with largest b + c = {}, strictly inside S, so the pair semigroup is primitive in this chart", images.join(" "), text(&widest));745 let edge = image(image((Frac::int(1), Frac::zero()), false), true);746 println!(747 "domino-gasket sends (1,0) to ({},{}) with b + c = 1, on the face, so length 3 is minimal",748 text(&edge.0),749 text(&edge.1)750 );751 let mut positive = 0usize;752 let mut products = 0usize;753 for length in 1..=12usize {754 for mask in 0..(1usize << length) {755 let mut matrix: Vec<Vec<i128>> = (0..4)756 .map(|i| (0..4).map(|j| if i == j { 1i128 } else { 0 }).collect())757 .collect();758 for i in 0..length {759 let code = if (mask >> i) & 1 == 1 { 7u8 } else { 3u8 };760 let next = &rep.matrices[&code];761 let mut out = vec![vec![0i128; 4]; 4];762 for r in 0..4 {763 for k in 0..4 {764 if matrix[r][k] == 0 {765 continue;766 }767 for c in 0..4 {768 out[r][c] += matrix[r][k] * next[k][c] as i128;769 }770 }771 }772 matrix = out;773 }774 products += 1;775 if matrix.iter().all(|row| row.iter().all(|entry| *entry > 0)) {776 positive += 1;777 }778 }779 }780 let negatives = |code: u8| {781 rep.matrices[&code]782 .iter()783 .flat_map(|row| row.iter())784 .filter(|entry| **entry < 0)785 .count()786 };787 println!(788 "entrywise positivity in the standard basis is the wrong test and fails: M_3 carries {} negative entries and M_7 carries {}, and {positive} of the {products} products of {{M_3, M_7}} of length 1 to 12 are entrywise positive",789 negatives(3),790 negatives(7)791 );792 assert_eq!(positive, 0, "no product is entrywise positive");793 let mut power: Vec<Vec<i128>> = (0..4)794 .map(|i| (0..4).map(|j| if i == j { 1i128 } else { 0 }).collect())795 .collect();796 let mut norms: Vec<String> = Vec::new();797 let mut wrong = 0usize;798 for length in 1..=32usize {799 let next = &rep.matrices[&3];800 let mut out = vec![vec![0i128; 4]; 4];801 for r in 0..4 {802 for k in 0..4 {803 if power[r][k] == 0 {804 continue;805 }806 for c in 0..4 {807 out[r][c] += power[r][k] * next[k][c] as i128;808 }809 }810 }811 power = out;812 let peak = power813 .iter()814 .flat_map(|row| row.iter())815 .map(|entry| entry.abs())816 .max()817 .expect("a matrix has entries");818 if peak != (1i128 << (length + 2)) - 2 {819 wrong += 1;820 }821 if [1usize, 2, 4, 8, 16, 32].contains(&length) {822 norms.push(format!("{peak}"));823 }824 if value(rep, &vec![3u8; length]) != 1 {825 wrong += 1;826 }827 }828 println!("the norm exponent is a different number from the component exponent: max |entry(M_3^L)| = 2^(L+2) - 2 reads {} at L = 1, 2, 4, 8, 16, 32, so the norm rate along 3^inf is log 2, while comp(A_(3^L)) = 1 at every L to 32 and the component rate is 0; mismatches {wrong}", norms.join(", "));829 assert_eq!(wrong, 0, "the norm witness is exact");830}831832fn boundary(rep: &Rep) {833 println!();834 println!("THE BOUNDARY FREQUENCY (1,0) OVER (3,7)");835 let scale = 3f64.log2();836 let squares = |i: usize| {837 let root = (i as f64).sqrt() as usize;838 (root.saturating_sub(1)..root + 2).any(|r| r * r == i)839 };840 let powers: Vec<u8> = (1..=FAR)841 .map(|i| if i.is_power_of_two() { 7 } else { 3 })842 .collect();843 let boxes: Vec<u8> = (1..=FAR).map(|i| if squares(i) { 7 } else { 3 }).collect();844 let flats = vec![3u8; FAR];845 let mut bad = 0usize;846 for word in [&powers, &boxes, &flats] {847 for length in 1..=DEEP {848 if value(rep, &word[..length]) != count(Shape::GasketDomino, &word[..length], 7, 3) {849 bad += 1;850 }851 }852 }853 assert_eq!(bad, 0, "the closed form covers the three boundary words");854 println!("three words whose gasket frequency is 0, closed form against the representation to L = {DEEP}, mismatches {bad}");855 let run = Track::new(&powers, 7);856 let marks = [2048usize, 2049, 4096, 4097, 8192, 8193, 16384, 16385, 32768];857 let rates = marks858 .iter()859 .map(|length| {860 format!(861 "{:.9} ({length})",862 run.log2_comp(*length, scale) / *length as f64863 )864 })865 .collect::<Vec<_>>()866 .join(", ");867 println!("the gasket at the powers of 2, prefix rate in log 2 units, floats: {rates}");868 let mut low = (f64::MAX, 0usize);869 let mut high = (0.0f64, 0usize);870 for length in 4097..=8192usize {871 let rate = run.log2_comp(length, scale) / length as f64;872 if rate < low.0 {873 low = (rate, length);874 }875 if rate > high.0 {876 high = (rate, length);877 }878 }879 println!(880 "over the single block 4097 <= L <= 8192 the rate sweeps from {:.5} at L = {} down to {:.5} at L = {}, floats, so the accumulation set is the whole interval and not two points",881 high.0, high.1, low.0, low.1882 );883 let boxed = Track::new(&boxes, 7);884 let squared = [4096usize, 8192, 16384, 32768]885 .iter()886 .map(|length| {887 format!(888 "{:.9} ({length})",889 boxed.log2_comp(*length, scale) / *length as f64890 )891 })892 .collect::<Vec<_>>()893 .join(", ");894 println!(895 "the gasket at the squares, prefix rate in log 2 units, floats: {squared}, closing on 1"896 );897 let flat = Track::new(&flats, 7);898 let steady = (1..=FAR)899 .map(|length| flat.log2_comp(length, scale))900 .fold(0.0f64, f64::max);901 println!("the constant word 3^L has comp 1 and rate {steady} at every L to {FAR}");902 assert_eq!(steady, 0.0, "the constant word never grows");903}904905fn by_products(rep: &Rep) {906 println!();907 println!("BY-PRODUCTS");908 let mut terms: Vec<String> = Vec::new();909 let mut bad = 0usize;910 for k in 1..=8usize {911 let word: Vec<u8> = (0..2 * k)912 .map(|i| if i % 2 == 0 { 7u8 } else { 3u8 })913 .collect();914 let want = (6i128.pow(k as u32) + 4) / 5;915 let seen = count(Shape::GasketDomino, &word, 7, 3);916 if seen != want || value(rep, &word) != want {917 bad += 1;918 }919 if k <= 6 {920 terms.push(format!("{seen}"));921 }922 }923 println!("comp(A_((7,3)^k)) = (6^k + 4)/5, reading {} at k = 1 to 6, checked to k = 8 against the closed form and the representation, mismatches {bad}", terms.join(", "));924 assert_eq!(bad, 0, "the periodic word has that closed form");925 let mut best = 0i128;926 let mut witness: Vec<u8> = Vec::new();927 for mask in 0..(1usize << 8) {928 let word = pair_word(7, 3, mask, 8);929 let seen = count(Shape::GasketDomino, &word, 7, 3);930 if seen > best {931 best = seen;932 witness = word.clone();933 }934 }935 println!(936 "the largest component count at L = 8 over (3,7) is {best} at the word {}, and 1094 = 2 x 547, while every closed form on the other 89 pairs gives a count of the form 2^a 3^b",937 witness938 .iter()939 .map(|code| format!("{code}"))940 .collect::<Vec<_>>()941 .join(",")942 );943 let scale = 3f64.log2();944 let mut block: Vec<u8> = vec![3, 7];945 while block.len() < 4096 {946 let size = block.len();947 let mut next = block.clone();948 next.extend(vec![7u8; size]);949 next.extend(vec![3u8; size]);950 block = next;951 }952 let track = Track::new(&block, 7);953 let mut low = f64::MAX;954 let mut high = 0.0f64;955 let mut thin = 1.0f64;956 for length in 1024..=4096usize {957 let rate = track.log2_comp(length, scale) / length as f64;958 low = low.min(rate);959 high = high.max(rate);960 let heavy = track.three[length] as f64 / length as f64;961 thin = thin.min(heavy.min(1.0 - heavy));962 }963 println!("the tripling word W_(k+1) = W_k 7^|W_k| 3^|W_k| at seed 3,7 has both letters at density at least {thin:.4} over 1024 <= L <= 4096 yet its prefix rate ranges over [{low:.4}, {high:.4}] in log 2 units with no narrowing, all floats");964}965966pub fn study(rep: &Rep) {967 println!();968 println!("THE FORTY-SIX");969 forms(rep);970 ledger();971 morse_value(rep);972 sandwich();973 cone(rep);974 boundary(rep);975 by_products(rep);976}