# Independent re-run of the exhaustive point search of the heuristic node (gluings of the shifted torus by W(E6)).
import sys, pickle; sys.path.insert(0, '/private/tmp/claude-501/-Users-miltonlin-Projects-Trees/75c2c3d9-ef06-4644-b59b-d12583740b0c/scratchpad/anim/work/G5/scripts')
from fractions import Fraction as F
from math import gcd
from collections import Counter
from e6model import EXP
A = EXP; n = 6
W = pickle.load(open('/private/tmp/claude-501/-Users-miltonlin-Projects-Trees/75c2c3d9-ef06-4644-b59b-d12583740b0c/scratchpad/anim/work/G5/W_E6.pkl', 'rb'))
rhoL_on = [-3, 1, 1, 1, 1, -3]
def rhoL(c): return sum(rhoL_on[j] * c[j] for j in range(n))
def winv(w):
    # w^{-1} = A^{-1} w^T A  (w preserves the form lam^T A mu); columns of w are w(alpha_j^vee)
    M = [[w[j][i] for j in range(n)] for i in range(n)]   # M[i][j] = i-th coord of w(alpha_j)
    # inverse via Gauss-Jordan
    X = [[F(M[i][j]) for j in range(n)] + [F(int(i == j)) for j in range(n)] for i in range(n)]
    for c in range(n):
        p = next(r for r in range(c, n) if X[r][c] != 0); X[c], X[p] = X[p], X[c]
        pv = X[c][c]; X[c] = [x / pv for x in X[c]]
        for r in range(n):
            if r != c and X[r][c] != 0:
                f = X[r][c]; X[r] = [x - f * y for x, y in zip(X[r], X[c])]
    Mi = [[X[i][n + j] for j in range(n)] for i in range(n)]
    return tuple(tuple(int(Mi[i][j]) for i in range(n)) for j in range(n))
def apply(w, c): return tuple(sum(w[j][i] * c[j] for j in range(n)) for i in range(n))
E = [tuple(int(i == j) for i in range(n)) for j in range(n)]
def wx_orbit_pt(p):
    # p = (su, au, sv, av): u = q^su e^{2 pi i au}, v = q^sv e^{2 pi i av}
    def s1(p): su, au, sv, av = p; return (-su, (-au) % 1, su + sv, (au + av) % 1)
    def s6(p): su, au, sv, av = p; return (su + sv, (au + av) % 1, -sv, (-av) % 1)
    o = {p}; fr = [p]
    while fr:
        nf = []
        for x in fr:
            for y in (s1(x), s6(x)):
                if y not in o: o.add(y); nf.append(y)
        fr = nf
    return frozenset(o)
stats = Counter(); glued = []; curve_ws = []
for w in W:
    wi = winv(w)
    cs = [apply(wi, E[j]) for j in (1, 2, 3, 4)]            # w^{-1} alpha_j^vee, black j
    eqs = [(c[0], c[5], 1 - rhoL(c)) for c in cs]            # u^{c1} v^{c6} = q^{1 - <rho_L, c>}
    Mm = [(a, b) for a, b, _ in eqs]
    minors = [Mm[i][0] * Mm[j][1] - Mm[i][1] * Mm[j][0] for i in range(4) for j in range(i + 1, 4)]
    N = 0
    for m in minors: N = gcd(N, abs(m))
    if N == 0:
        if all(a == 0 and b == 0 for a, b in Mm):
            stats['rank0'] += 1
            if all(k == 0 for _, _, k in eqs): stats['rank0_consistent'] += 1
            continue
        stats['rank1'] += 1
        # consistency of the real system: rows proportional; need consistent right-hand sides
        rows = [(a, b, k) for a, b, k in eqs if (a, b) != (0, 0)]
        if any(k != 0 for a, b, k in eqs if (a, b) == (0, 0)): stats['rank1_inconsistent'] += 1; continue
        a0, b0, _ = rows[0]; g0 = gcd(abs(a0), abs(b0)); p0, q0 = a0 // g0, b0 // g0
        ratios = set()
        for a, b, k in rows:
            t = a // p0 if p0 != 0 else b // q0
            ratios.add(F(k, t))
        if len(ratios) > 1: stats['rank1_inconsistent'] += 1
        else: stats['rank1_curve'] += 1; curve_ws.append(w)
        continue
    stats['rank2'] += 1
    i, j = next((i, j) for i in range(4) for j in range(i + 1, 4) if Mm[i][0] * Mm[j][1] - Mm[i][1] * Mm[j][0] != 0)
    det = F(Mm[i][0] * Mm[j][1] - Mm[i][1] * Mm[j][0])
    su = (Mm[j][1] * eqs[i][2] - Mm[i][1] * eqs[j][2]) / det
    sv = (-Mm[j][0] * eqs[i][2] + Mm[i][0] * eqs[j][2]) / det
    if any(a * su + b * sv != k for a, b, k in eqs): stats['rank2_inconsistent'] += 1; continue
    stats['rank2_consistent'] += 1
    sols = [(F(ku, N), F(kv, N)) for ku in range(N) for kv in range(N) if all((a * ku + b * kv) % N == 0 for a, b, _ in eqs)]
    for au, av in sols:
        def xval(c):  # log_q|x(c)| and arg/2pi of x(c) for x = q^{rho_L} y
            return (rhoL(c) + c[0] * su + c[5] * sv, (c[0] * au + c[5] * av) % 1)
        c1 = apply(wi, E[0]); c6 = apply(wi, E[5])
        l1, a1 = xval(c1); l6, a6 = xval(c6)
        yp = (l1 + 3, a1, l6 + 3, a6)        # y'(e1) = x(w^{-1} alpha_1^vee) q^{3}, likewise e6
        y = (su, au, sv, av)
        stats['rank2_solutions'] += 1
        if yp not in wx_orbit_pt(y):
            glued.append((y, yp, w))
print(dict(stats))
print('non-W_X-conjugate gluings (w, y) found:', len(glued))
pairs = Counter()
for y, yp, w in glued:
    pairs[(min(wx_orbit_pt(y)), min(wx_orbit_pt(yp)))] += 1
for k, v in pairs.items(): print('   orbit pair', k, ' count', v)
ys = Counter(y for y, yp, w in glued); print('distinct glued points y:', len(ys))
pickle.dump(curve_ws, open('/private/tmp/claude-501/-Users-miltonlin-Projects-Trees/75c2c3d9-ef06-4644-b59b-d12583740b0c/scratchpad/anim/work/G5/curve_ws.pkl', 'wb'))