"""JAM 3, line 3 instrument (the other flash): CONWAY'S MATRIX, EXACT.
No sampling anywhere. Atoms found by empirical decay-splitting (a
boundary is a split iff the halves evolve independently for T steps:
las^t(L) + las^t(R) == las^t(L+R), t = 1..T), leftmost-first
recursion (sound: if s = L+R splits and L = A+B splits, then s splits
at A too - associativity of non-interaction). Closure from seed "1":
decompose every atom's decay until no new atoms appear. Then the
integer transition matrix M (row = atom, entries = child counts) and
its exact eigenvalues. AUDIT GATE built in: the Perron root must
reproduce Conway's constant 1.303577269 or the whole object is
rejected. lambda2 answers the sheet's question directly."""
from itertools import groupby
from functools import lru_cache
import numpy as np

T = 20


@lru_cache(maxsize=None)
def las(s: str) -> str:
    return "".join(str(len(list(g))) + k for k, g in groupby(s))


@lru_cache(maxsize=None)
def splits_at(L: str, R: str) -> bool:
    l, r, w = L, R, L + R
    for _ in range(T):
        l, r, w = las(l), las(r), las(w)
        if l + r != w:
            return False
    return True


@lru_cache(maxsize=None)
def decompose(s: str) -> tuple:
    for i in range(1, len(s)):
        if s[i - 1] != s[i] and splits_at(s[:i], s[i:]):
            return (s[:i],) + decompose(s[i:])
    return (s,)


def main():
    # grow a seed string long enough to contain several atoms, then
    # close the atom set under decay
    s = "1"
    for _ in range(24):
        s = las(s)
    atoms = set(decompose(s))
    frontier = list(atoms)
    while frontier:
        a = frontier.pop()
        for child in decompose(las(a)):
            if child not in atoms:
                atoms.add(child)
                frontier.append(child)
    atoms = sorted(atoms, key=lambda a: (len(a), a))
    n = len(atoms)
    idx = {a: i for i, a in enumerate(atoms)}
    print("atoms found by decay-splitting closure:", n)
    M = np.zeros((n, n))
    for a in atoms:
        for child in decompose(las(a)):
            M[idx[a], idx[child]] += 1.0
    ev = np.linalg.eigvals(M)
    ev = sorted(ev, key=lambda z: -abs(z))
    l1, l2 = ev[0], ev[1]
    print("lambda1 = %.9f  (Conway's constant = 1.303577269; gate %s)"
          % (abs(l1), "PASS" if abs(abs(l1) - 1.303577269) < 1e-6 else "FAIL"))
    print("lambda2 = %.6f%+.6fj  |lambda2| = %.6f  angle = %.4f rad"
          % (l2.real, l2.imag, abs(l2), abs(np.angle(l2))))
    print("|lambda2|/lambda1 = %.6f   complex pair: %s"
          % (abs(l2) / abs(l1), abs(l2.imag) > 1e-9))
    print("next moduli:", "  ".join("%.4f" % abs(z) for z in ev[2:8]))


if __name__ == "__main__":
    main()