#!/usr/bin/env python3
"""Verify the 2-generator presentation and clockwork lifts for *632 (p6m, g268..g271).

Checks:
  1. In the render frame of g268 (the plain 1-colouring of *632 / p6m), the two
     generators
       alpha = P Q  (counterclockwise 1/6-turn about the 6-centre at (0, 0))
       R            (vertical mirror across x = -1/2, halfway to the next 6-centre)
     generate all three Coxeter triangle mirrors P, Q, R and both primitive
     translations t_1, t_2 via explicit short words:
       R   = R
       Q   = R alpha^-1 R alpha R
       P   = alpha R alpha^-1 R alpha R
       t_1 = alpha R alpha R
       t_2 = alpha^-1 R alpha^-1 R
     and therefore generate all 12 point-group cosets and the full lattice Z^2.

  2. The 2-generator presentation
       Gamma = < alpha, R | alpha^6 = R^2 = (alpha R alpha^-1 R)^3
                          = (alpha R alpha^-1 R alpha R)^2 = 1 >
     is Tietze-equivalent to the Coxeter presentation
       < P, Q, R | P^2 = Q^2 = R^2 = (P Q)^6 = (Q R)^3 = (R P)^2 = 1 >.
     In particular, (R P)^2 = 1 follows identically in the free group from
     (alpha R alpha^-1 R)^3 = 1 and P^2 = R^2 = 1 because
       P R P R = (alpha R alpha^-1 R)^3  when  P = alpha R alpha^-1 R alpha R.

  3. Every forward clockwork homomorphism tau : Gamma -> R/Z on the 2-generator
     presentation satisfies 2*tau(alpha) = 0 and 2*tau(R) = 0 (mod 1), giving
     an unconstrained 2 x 2 table (tau(alpha), tau(R)) in {0, 1/2}^2 matching
     the four forward clockwork groups of p6m:
       (0,   0  ) -> g268 (*632,       P6mm,  No. 183)
       (0,   1/2) -> g270 (*~6~3~2,    P6cc,  No. 184)
       (1/2, 1/2) -> g269 (*6_3~3~2_1, P6_3cm, No. 185)
       (1/2, 0  ) -> g271 (*~6_332_1,  P6_3mc, No. 186)
"""

from __future__ import annotations
from collections import deque
from fractions import Fraction

Mat = tuple[tuple[int, int], tuple[int, int]]
Vec = tuple[Fraction, Fraction]
Aff = tuple[Mat, Vec]

I2: Mat = ((1, 0), (0, 1))
ZERO: Vec = (Fraction(0), Fraction(0))
IDENTITY: Aff = (I2, ZERO)


def mat_mul(a: Mat, b: Mat) -> Mat:
    return (
        (a[0][0] * b[0][0] + a[0][1] * b[1][0], a[0][0] * b[0][1] + a[0][1] * b[1][1]),
        (a[1][0] * b[0][0] + a[1][1] * b[1][0], a[1][0] * b[0][1] + a[1][1] * b[1][1]),
    )


def mat_vec(a: Mat, v: Vec) -> Vec:
    return (
        a[0][0] * v[0] + a[0][1] * v[1],
        a[1][0] * v[0] + a[1][1] * v[1],
    )


def aff_mul(g: Aff, h: Aff) -> Aff:
    m1, v1 = g
    m2, v2 = h
    mv = mat_vec(m1, v2)
    return (mat_mul(m1, m2), (mv[0] + v1[0], mv[1] + v1[1]))


def aff_inv(g: Aff) -> Aff:
    m, v = g
    det = m[0][0] * m[1][1] - m[0][1] * m[1][0]
    assert det in (1, -1)
    mi: Mat = (
        (m[1][1] * det, -m[0][1] * det),
        (-m[1][0] * det, m[0][0] * det),
    )
    mv = mat_vec(mi, v)
    return (mi, (-mv[0], -mv[1]))


def aff_pow(g: Aff, n: int) -> Aff:
    if n < 0:
        return aff_pow(aff_inv(g), -n)
    out = IDENTITY
    for _ in range(n):
        out = aff_mul(out, g)
    return out


def eval_word(gens: dict[str, Aff], word: list[str]) -> Aff:
    out = IDENTITY
    for symbol in word:
        if symbol.endswith("^-1"):
            out = aff_mul(out, aff_inv(gens[symbol[:-3]]))
        else:
            out = aff_mul(out, gens[symbol])
    return out


def free_reduce(word: list[tuple[str, int]]) -> list[tuple[str, int]]:
    stack: list[tuple[str, int]] = []
    for sym, exp in word:
        if exp == 0:
            continue
        if stack and stack[-1][0] == sym:
            new_exp = stack[-1][1] + exp
            if new_exp == 0:
                stack.pop()
            else:
                stack[-1] = (sym, new_exp)
        else:
            stack.append((sym, exp))
    return stack


def main() -> None:
    # Coxeter triangle mirrors P, Q, R on the g268 plate (in lattice coordinates
    # with basis b_1 = (1/2, -sqrt(3)/2), b_2 = (1/2, +sqrt(3)/2)):
    #   P: reflection across y = 0       (M = [[0, 1], [1, 0]], v = (0, 0))
    #   Q: reflection across -30 deg ray (M = [[1, 0], [1, -1]], v = (0, 0))
    #   R: reflection across x = -1/2    (M = [[0, -1], [-1, 0]], v = (-1, -1))
    P: Aff = (((0, 1), (1, 0)), (Fraction(0), Fraction(0)))
    Q: Aff = (((1, 0), (1, -1)), (Fraction(0), Fraction(0)))
    R: Aff = (((0, -1), (-1, 0)), (Fraction(-1), Fraction(-1)))

    # Our 2 generators:
    #   alpha = P Q (1/6-turn CCW about (0, 0))
    #   R           (vertical mirror across x = -1/2, halfway to (-1, 0))
    alpha: Aff = aff_mul(P, Q)
    assert alpha == (((1, -1), (1, 0)), (Fraction(0), Fraction(0)))

    gens = {"α": alpha, "R": R}

    # 1. Verify explicit words for R, Q, P, t_1, t_2 in {alpha, R}:
    w_R = ["R"]
    w_Q = ["R", "α^-1", "R", "α", "R"]
    w_P = ["α", "R", "α^-1", "R", "α", "R"]
    w_t1 = ["α", "R", "α", "R"]
    w_t2 = ["α^-1", "R", "α^-1", "R"]

    assert eval_word(gens, w_R) == R
    assert eval_word(gens, w_Q) == Q
    assert eval_word(gens, w_P) == P
    assert eval_word(gens, w_t1) == (I2, (Fraction(1), Fraction(0)))
    assert eval_word(gens, w_t2) == (I2, (Fraction(0), Fraction(1)))

    # 2. Verify all 4 defining relators of Gamma = <alpha, R | ...>:
    rel_1 = aff_pow(alpha, 6)
    rel_2 = aff_pow(R, 2)
    rel_3 = aff_pow(eval_word(gens, ["α", "R", "α^-1", "R"]), 3)
    rel_4 = aff_pow(eval_word(gens, w_P), 2)
    assert rel_1 == IDENTITY
    assert rel_2 == IDENTITY
    assert rel_3 == IDENTITY
    assert rel_4 == IDENTITY

    # Verify the Coxeter relations on P, Q, R and the free-group identity
    # P R P^-1 R == (alpha R alpha^-1 R)^3 when R^2 = 1:
    assert aff_pow(P, 2) == IDENTITY
    assert aff_pow(Q, 2) == IDENTITY
    assert aff_pow(R, 2) == IDENTITY
    assert aff_pow(aff_mul(P, Q), 6) == IDENTITY
    assert aff_pow(aff_mul(Q, R), 3) == IDENTITY
    assert aff_pow(aff_mul(R, P), 2) == IDENTITY

    # Check in the quotient by R^2 = 1 that P R P^-1 R equals (alpha R alpha^-1 R)^3:
    # P = a R a^-1 R a R  =>  P^-1 = R a^-1 R a R a^-1  (using R^-1 = R)
    # P R P^-1 R = (a R a^-1 R a R) R (R a^-1 R a R a^-1) R
    #            = a R a^-1 R a R a^-1 R a R a^-1 R = (a R a^-1 R)^3.
    p_word = [("α", 1), ("R", 1), ("α", -1), ("R", 1), ("α", 1), ("R", 1)]
    p_inv_word = [("R", -1), ("α", -1), ("R", -1), ("α", 1), ("R", -1), ("α", -1)]
    # Substitute R^-1 -> R (since R^2 = 1):
    lhs = [(s, 1 if s == "R" else e) for s, e in (p_word + [("R", 1)] + p_inv_word + [("R", 1)])]
    # Reduce pairs of R*R -> empty and cancel α*α^-1:
    changed = True
    while changed:
        changed = False
        reduced: list[tuple[str, int]] = []
        for sym, exp in lhs:
            if reduced and reduced[-1][0] == sym:
                new_exp = reduced[-1][1] + exp
                if sym == "R":
                    new_exp %= 2
                if new_exp == 0:
                    reduced.pop()
                else:
                    reduced[-1] = (sym, new_exp)
                changed = True
            else:
                reduced.append((sym, exp))
        lhs = reduced
    rhs = [("α", 1), ("R", 1), ("α", -1), ("R", 1)] * 3
    assert lhs == rhs

    # 3. Breadth-first generation check from {alpha, alpha^-1, R}:
    seen: dict[Aff, int] = {IDENTITY: 0}
    queue: deque[Aff] = deque([IDENTITY])
    steps = [alpha, aff_inv(alpha), R]
    while queue:
        cur = queue.popleft()
        d = seen[cur]
        if d >= 6:
            continue
        for s in steps:
            nxt = aff_mul(cur, s)
            if nxt not in seen:
                seen[nxt] = d + 1
                queue.append(nxt)

    point_group = {m for m, _ in seen}
    translations = {v for m, v in seen if m == I2}
    assert len(point_group) == 12
    for req in [(Fraction(1), Fraction(0)), (Fraction(-1), Fraction(0)),
                (Fraction(0), Fraction(1)), (Fraction(0), Fraction(-1))]:
        assert req in translations

    # 4. Classify all forward clockwork homomorphisms tau : Gamma -> R/Z:
    # Relators give:
    #   6 * tau(alpha) in Z
    #   2 * tau(R)     in Z
    #   3 * (tau(alpha) + tau(R) - tau(alpha) + tau(R)) = 6 * tau(R) in Z
    #   2 * (tau(alpha) + tau(R) - tau(alpha) + tau(R) + tau(alpha) + tau(R))
    #     = 2 * tau(alpha) + 6 * tau(R) in Z  ==>  2 * tau(alpha) in Z.
    valid_lifts: list[tuple[Fraction, Fraction, Fraction, Fraction]] = []
    for a_num in range(6):
        for r_num in range(2):
            tau_a = Fraction(a_num, 6)
            tau_r = Fraction(r_num, 2)
            rel1_ok = (6 * tau_a).denominator == 1
            rel2_ok = (2 * tau_r).denominator == 1
            rel3_ok = (6 * tau_r).denominator == 1
            rel4_ok = (2 * tau_a + 6 * tau_r).denominator == 1
            if rel1_ok and rel2_ok and rel3_ok and rel4_ok:
                tau_q = (3 * tau_r) % 1
                tau_p = (tau_a + tau_q) % 1
                valid_lifts.append((tau_a, tau_r, tau_p, tau_q))

    assert len(valid_lifts) == 4
    expected = {
        (Fraction(0), Fraction(0)): ("g268", "*632", "*632", "P6mm", "No. 183"),
        (Fraction(0), Fraction(1, 2)): ("g270", "*²6²3²2", "*~6~3~2", "P6cc", "No. 184"),
        (Fraction(1, 2), Fraction(1, 2)): ("g269", "*¹6²3²2", "*6₃~3~2₁", "P6₃cm", "No. 185"),
        (Fraction(1, 2), Fraction(0)): ("g271", "*²6¹3¹2", "*~6₃32₁", "P6₃mc", "No. 186"),
    }

    print("2-Generator Verification for *632 (p6m)")
    print("========================================")
    print("Generators in g268 plate frame (basis b1=(1/2, -sqrt(3)/2), b2=(1/2, +sqrt(3)/2)):")
    print(f"  α = P Q = {alpha[0]}, v = ({alpha[1][0]}, {alpha[1][1]})  [1/6-turn CCW about (0, 0)]")
    print(f"  R       = {R[0]}, v = ({R[1][0]}, {R[1][1]})  [mirror across x = -1/2]")
    print()
    print("Reconstructed Coxeter triangle mirrors and primitive lattice translations:")
    print(f"  R   = R                 (length 1) -> {eval_word(gens, w_R)}")
    print(f"  Q   = R α^-1 R α R      (length 5) -> {eval_word(gens, w_Q)}")
    print(f"  P   = α R α^-1 R α R    (length 6) -> {eval_word(gens, w_P)}")
    print(f"  t_1 = (α R)²            (length 4) -> v = (+1, 0)")
    print(f"  t_2 = (α^-1 R)²         (length 4) -> v = (0, +1)")
    print(f"  Point-group cosets reached within word length <= 6: {len(point_group)} / 12")
    print()
    print("Presentation:")
    print("  Γ = ⟨α, R | α⁶ = R² = (αRα⁻¹R)³ = (αRα⁻¹RαR)² = 1⟩")
    print("  Free-group identity modulo R² = 1: (PRP⁻¹R) ≡ (αRα⁻¹R)³  [verified]")
    print()
    print("Classification of forward clockwork homomorphisms τ : Γ -> R/Z:")
    print("  τ(α)  τ(R)  | induced τ(P)  τ(Q)  τ(R) | group  Chaim     clockwork   space group")
    print("  ------------+--------------------------+-----------------------------------------")
    for tau_a, tau_r, tau_p, tau_q in valid_lifts:
        gid, chaim, cw, hm, itc = expected[(tau_a, tau_r)]
        print(f"  {str(tau_a):<5} {str(tau_r):<5} |         {str(tau_p):<5} {str(tau_q):<5} {str(tau_r):<5}| "
              f"{gid:<5}  {chaim:<8}  {cw:<10}  {hm} ({itc})")
    print()
    print("All checks passed.")


if __name__ == "__main__":
    main()