"""Exact finite verification for the manuscript on Bergman's third unit.

This script uses only Python's standard library. The cocycle models use F_2
coordinates; the independent generalized-dihedral model uses Z/4 coordinates.
It verifies P, H, the equivariant section, the model E=P x C_2^3,
the quotient E/K[2], and the explicit third unit. The universal properties
and the identification of E are proved in the accompanying manuscript.
"""

from itertools import product


# Coordinates of P are (v,w) in F_2^4 x F_2^3.
# The W-basis a,b,c is encoded by the bit values 1,2,4.
P_COMM = {
    (0, 1): 1, (0, 2): 2, (0, 3): 3,
    (1, 2): 4, (1, 3): 5, (2, 3): 6,
}


def p_cocycle(v, u):
    out = 0
    for (i, j), value in P_COMM.items():
        if (v >> j) & 1 and (u >> i) & 1:
            out ^= value
    return out


def p_mul(g, h):
    v, w = g
    u, z = h
    return v ^ u, w ^ z ^ p_cocycle(v, u)


def p_inv(g):
    v, w = g
    return v, w ^ p_cocycle(v, v)


def p_conj(g, h):
    return p_mul(p_mul(g, h), p_inv(g))


# H has coordinates F_2^3 x F_2^2; [h1,h2]=a,
# [h1,h3]=b and [h2,h3]=a+b.
H_COMM = {(0, 1): 1, (0, 2): 2, (1, 2): 3}


def h_cocycle(v, u):
    out = 0
    for (i, j), value in H_COMM.items():
        if (v >> j) & 1 and (u >> i) & 1:
            out ^= value
    return out


def h_mul(g, h):
    v, w = g
    u, z = h
    return v ^ u, w ^ z ^ h_cocycle(v, u)


def h_inv(g):
    v, w = g
    return v, w ^ h_cocycle(v, v)


def h_conj(g, h):
    return h_mul(h_mul(g, h), h_inv(g))


def projection(g):
    """P -> H, xi -> hi (i <= 3), x4 -> h1*h2*h3."""
    v, w = g
    fourth = (v >> 3) & 1
    first_three = v & 7
    projected_w = (w & 3) ^ (3 if (w & 4) else 0)
    projected_w ^= fourth * (((v >> 1) ^ (v >> 2)) & 1)
    return first_three ^ (7 if fourth else 0), projected_w


def section(t):
    """A conjugation-equivariant section Omega(H) -> Omega(P)."""
    v, w = t
    a, b = w & 1, (w >> 1) & 1
    if v in (0, 1):
        return v, a | (b << 1)
    if v == 2:
        return 2, (a ^ b) | (b << 2)
    if v == 4:
        return 4, ((a ^ b) << 1) | (a << 2)
    if v == 7:
        return 8, b | (a << 1) | ((a ^ b) << 2)
    raise ValueError("The argument does not belong to Omega(H)")


def e_mul(g, h):
    return p_mul(g[0], h[0]), g[1] ^ h[1]


def e_inv(g):
    return p_inv(g[0]), g[1]


def e_label(t):
    v, w = t
    tag = {1: 1, 2: 2, 3: 4}.get(w, 0) if v == 0 else 0
    return section(t), tag


def q_mul(g, h):
    return h_mul(g[0], h[0]), g[1] ^ h[1]


def quotient(g):
    p, _ = g
    return projection(p), (p[0] >> 3) & 1


def d_mul(g, h):
    """Independent generalized-dihedral coordinates (i,j,epsilon)."""
    i, j, epsilon = g
    k, ell, eta = h
    sign = -1 if epsilon else 1
    return (i + sign * k) % 4, (j + sign * ell) % 4, epsilon ^ eta


def dihedral_coordinates(g):
    v, w = g
    out = (0, 0, 0)
    generators = [(0, 0, 1), (1, 0, 1), (0, 1, 1)]
    for i, generator in enumerate(generators):
        if (v >> i) & 1:
            out = d_mul(out, generator)
    if w & 1:
        out = d_mul(out, (2, 0, 0))
    if w & 2:
        out = d_mul(out, (0, 2, 0))
    return out


def generated(elements, mul, identity):
    subgroup = {identity}
    while True:
        new = subgroup | {mul(g, s) for g in subgroup for s in elements}
        if new == subgroup:
            return subgroup
        subgroup = new


def verify_cocycle(cocycle, vector_order):
    # This identity proves associativity for all elements, since the second
    # coordinates enter the multiplication additively and are central.
    for v, u, z in product(range(vector_order), repeat=3):
        assert cocycle(v, u) ^ cocycle(v ^ u, z) == (
            cocycle(u, z) ^ cocycle(v, u ^ z)
        )


def main():
    if not __debug__:
        raise RuntimeError("Run this verification program without -O or -OO.")
    one = (0, 0)
    p_elements = list(product(range(16), range(8)))
    h_elements = list(product(range(8), range(4)))
    verify_cocycle(p_cocycle, 16)
    verify_cocycle(h_cocycle, 8)
    for g in p_elements:
        assert p_mul(g, one) == p_mul(one, g) == g
        assert p_mul(g, p_inv(g)) == p_mul(p_inv(g), g) == one
    for g in h_elements:
        assert h_mul(g, one) == h_mul(one, g) == g
        assert h_mul(g, h_inv(g)) == h_mul(h_inv(g), g) == one
    for g, h in product(p_elements, repeat=2):
        assert projection(p_mul(g, h)) == h_mul(projection(g), projection(h))
    assert {projection(g) for g in p_elements} == set(h_elements)

    x = [(1 << i, 0) for i in range(4)]
    h_generators = [(1 << i, 0) for i in range(3)]
    r = one
    for t in x:
        r = p_mul(r, t)
    r2 = p_mul(r, r)
    assert r == (15, 0) and r2 == (0, 7)
    assert r2 != one and p_mul(r2, r2) == one
    assert all(p_conj(g, r) == r for g in p_elements)
    kernel = {g for g in p_elements if projection(g) == one}
    assert kernel == {one, r, r2, p_mul(r2, r)}
    assert generated(x, p_mul, one) == set(p_elements)
    assert generated(h_generators, h_mul, one) == set(h_elements)

    p_omega = [g for g in p_elements if p_mul(g, g) == one]
    h_omega = [g for g in h_elements if h_mul(g, g) == one]
    assert {v for v, w in p_omega} == {0, 1, 2, 4, 8}
    assert {v for v, w in h_omega} == {0, 1, 2, 4, 7}
    assert len(p_omega) == 40 and len(h_omega) == 20
    assert section(one) == one
    for t in h_omega:
        assert projection(section(t)) == t
        assert p_mul(section(t), section(t)) == one
    for g in p_elements:
        for t in h_omega:
            assert section(h_conj(projection(g), t)) == p_conj(g, section(t))
    for t, u in product(h_omega, repeat=2):
        assert section(h_conj(t, u)) == p_conj(section(t), section(u))
    assert all(section(projection(t)) == t for t in x)
    assert generated([section(t) for t in h_omega], p_mul, one) == set(p_elements)

    center_h = {g for g in h_elements if all(h_mul(g, h) == h_mul(h, g) for h in h_elements)}
    derived_h = generated([
        h_mul(h_conj(g, h), h_inv(h)) for g, h in product(h_elements, repeat=2)
    ], h_mul, one)
    assert len(center_h) == len(derived_h) == 4
    assert derived_h == center_h
    assert all(h_mul(g, g) == one for g in derived_h)

    # The explicit model E=P x C_2^3 and its labelled generators.
    e_one = (one, 0)
    e_elements = list(product(p_elements, range(8)))
    labels = [e_label(t) for t in h_omega]
    assert e_label(one) == e_one
    assert all(e_mul(t, t) == e_one for t in labels)
    for s, t in product(h_omega, repeat=2):
        assert e_label(h_conj(s, t)) == e_mul(
            e_mul(e_label(s), e_label(t)), e_label(s)
        )
    assert generated(labels, e_mul, e_one) == set(e_elements)
    e_omega = [g for g in e_elements if e_mul(g, g) == e_one]
    e_kernel = {g for g in e_elements if projection(g[0]) == one}
    k_two = {g for g in e_kernel if e_mul(g, g) == e_one}
    assert len(e_elements) == 1024 and len(e_omega) == 320
    assert len(e_kernel) == 32 and len(k_two) == 16
    assert e_kernel == set(product(kernel, range(8)))
    assert k_two == set(product({one, r2}, range(8)))
    relators = {
        e_mul(g, e_inv(e_label(projection(g[0]))))
        for g in e_omega
    }
    assert relators == k_two

    # Q=H x C_2 and the quotient map E -> Q.
    q_one = (one, 0)
    q_elements = list(product(h_elements, range(2)))
    assert {quotient(g) for g in e_elements} == set(q_elements)
    assert {g for g in e_elements if quotient(g) == q_one} == k_two
    for g in e_elements:
        for t in labels:
            assert quotient(e_mul(g, t)) == q_mul(quotient(g), quotient(t))
    q_omega = [g for g in q_elements if q_mul(g, g) == q_one]
    unit = {t: quotient(e_label(t)) for t in h_omega}
    assert len(q_elements) == 64 and len(q_omega) == 40
    assert len(set(unit.values())) == 20
    assert all(unit[t] == (t, int(t[0] == 7)) for t in h_omega)
    assert (one, 1) not in set(unit.values())
    assert generated(list(unit.values()), q_mul, q_one) == set(q_elements)

    # Independent identification H = (Z/4)^2 semidirect Z/2.
    d_elements = set(product(range(4), range(4), range(2)))
    assert {dihedral_coordinates(g) for g in h_elements} == d_elements
    for g, h in product(h_elements, repeat=2):
        assert dihedral_coordinates(h_mul(g, h)) == d_mul(
            dihedral_coordinates(g), dihedral_coordinates(h)
        )
    for t in h_omega:
        i, j, _ = dihedral_coordinates(t)
        assert unit[t][1] == (i % 2) * (j % 2)

    print("All exact checks passed.")
    print("|P| = 128; |H| = 32; |ker(P -> H)| = 4.")
    print("r = x1*x2*x3*x4 = (15,0); r^2 = (0,7) != 1; r^4 = 1.")
    print("The kernel is central cyclic C4 generated by r.")
    print("|Omega(P)| = 40; |Omega(H)| = 20, including identity.")
    print("|Z(H)| = |[H,H]| = 4; [H,H] is elementary abelian.")
    print("The section Omega(H) -> Omega(P) preserves identity and conjugation.")
    print("Its image generates P and contains x1,x2,x3,x4.")
    print("|E| = 1024; |Omega(E)| = 320; the 20 labels generate E.")
    print("|K| = 32; |K[2]| = 16; the coequalizer relators are exactly K[2].")
    print("|Q| = 64; |Omega(Q)| = 40; ker(E -> Q) = K[2].")
    print("The third unit has 20 distinct values and misses (1,1) in H x C2.")
    print("In dihedral coordinates its last coordinate is (i mod 2)(j mod 2).")
    print("The unit image generates all of Q.")


if __name__ == "__main__":
    main()