import sys; sys.path.insert(0, "/private/tmp/claude-501/-Users-miltonlin-Projects-Trees/75c2c3d9-ef06-4644-b59b-d12583740b0c/scratchpad/solve/A6")
from e6_setup import *
from fractions import Fraction as Fr
from math import gcd
import pickle
# Weyl group as tuples of images of basis coroots (columns): w = (w(alpha_1^vee),...,w(alpha_6^vee))
I6=tuple(tuple(1 if i==j else 0 for j in range(n)) for i in range(n))
def apply(w, lam): return tuple(sum(w[j][i]*lam[j] for j in range(n)) for i in range(n))
def compose_refl(j, w): return tuple(refl(j,col) for col in w)   # s_j o w
W={I6}; fr=[I6]
while fr:
    nf=[]
    for w in fr:
        for j in range(n):
            w2=compose_refl(j,w)
            if w2 not in W: W.add(w2); nf.append(w2)
    fr=nf
print('|W|=',len(W))
# shifted torus point x(u,v): x(alpha_i^vee) = t^{<2rho_L,alpha_i^vee>} * y(phi alpha_i^vee), y(e1)=u, y(e6)=v, e1=phi a1, e6=phi a6.
# x(lam) for lam in coroot coords: t-exponent E(lam)=<2rho_L,lam>, u-exponent = lam_1, v-exponent = lam_6  (phi(lam) = lam_1 e1 + lam_6 e6 since phi kills black)
black=[1,2,3,4]
# inverse of w: find by search in set? compute via transposition w.r.t. form: w^{-1} = A^{-1} w^T A ; easier: build dict
def inv(w):
    # solve w*x = e_i for each i: since w is a permutation-free integer matrix, use the fact that W is a group: w^{-1} = unique element with w^{-1}(w e_i)=e_i; brute force by applying words is costly; use adjugate via fractions
    M=[[Fr(w[j][i]) for j in range(n)] for i in range(n)]  # M[i][j] = (w(alpha_j^vee))_i
    N=[row[:]+[Fr(1 if i==k else 0) for k in range(n)] for i,row in enumerate(M)]
    for c in range(n):
        p=next(r for r in range(c,n) if N[r][c]!=0); N[c],N[p]=N[p],N[c]
        pv=N[c][c]; N[c]=[x/pv for x in N[c]]
        for r in range(n):
            if r!=c and N[r][c]!=0:
                f=N[r][c]; N[r]=[x-f*y for x,y in zip(N[r],N[c])]
    Minv=[[N[i][n+j] for j in range(n)] for i in range(n)]
    return tuple(tuple(int(Minv[i][j]) for i in range(n)) for j in range(n))  # columns
def wx_images(u,v):
    # (u,v) as (exp_u, ang_u, exp_v, ang_v) with u = t^exp * e^{2 pi i ang}; W_X=S3 generated by s1:(u,v)->(u^-1, u v), s6:(u,v)->(u v, v^-1)
    def s1(p): eu,au,ev,av=p; return (-eu,(-au)%1,eu+ev,(au+av)%1)
    def s6(p): eu,au,ev,av=p; return (eu+ev,(au+av)%1,-ev,(-av)%1)
    orb={(u[0],u[1],v[0],v[1])}; fr=list(orb)
    while fr:
        nf=[]
        for p in fr:
            for q in (s1(p),s6(p)):
                if q not in orb: orb.add(q); nf.append(q)
        fr=nf
    return orb
glued=[]; curves=[]; count_norm=0
for w in W:
    wi=inv(w)
    cols=[apply(wi,I6[i]) for i in range(n)]   # c_i = w^{-1} alpha_i^vee
    # condition for black i: x(c_i) = t^2  ->  u^{c_i1} v^{c_i6} = t^{2 - E(c_i)}
    eqs=[(c[0],c[5],2-two_rhoL_pair(c)) for c in [cols[i] for i in black]]
    rank=0
    M=[(a,b) for a,b,_ in eqs]
    minors=[M[i][0]*M[j][1]-M[i][1]*M[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 any(a!=0 or b!=0 for a,b in M): rank=1
        else: rank=0
    else: rank=2
    if rank==0:
        count_norm+=1; continue   # w normalises the plane: W_X-type identification
    if rank==1:
        curves.append((w,eqs)); continue
    # rank 2: solve real exponents (a,b) from two independent equations, check all
    i,j=next((i,j) for i in range(4) for j in range(i+1,4) if M[i][0]*M[j][1]-M[i][1]*M[j][0]!=0)
    det=Fr(M[i][0]*M[j][1]-M[i][1]*M[j][0])
    ea=( M[j][1]*eqs[i][2]-M[i][1]*eqs[j][2])/det
    eb=(-M[j][0]*eqs[i][2]+M[i][0]*eqs[j][2])/det
    if any(Fr(a)*ea+Fr(b)*eb!=c for a,b,c in eqs): continue   # inconsistent: no solution
    # roots of unity: (zu,zv) in mu_N^2 with zu^a zv^b = 1 for all rows
    sols=[]
    for ku in range(N):
        for kv in range(N):
            if all(((a*ku+b*kv)%N)==0 for a,b,_ in eqs):
                sols.append((Fr(ku,N),Fr(kv,N)))
    for au,av in sols:
        u=(ea,au); v=(eb,av)
        # image point: u' = x(w^{-1} alpha_1^vee) * t^{6}, v' = x(w^{-1} alpha_6^vee) * t^6
        def xval(c):   # returns (t-exponent, angle)
            return (two_rhoL_pair(c)+c[0]*ea+c[5]*eb, (c[0]*au+c[5]*av)%1)
        e1,a1=xval(cols[0]); e6,a6=xval(cols[5])
        up=(e1+6,a1); vp=(e6+6,a6)
        if (up[0],up[1],vp[0],vp[1]) not in wx_images(u,v):
            glued.append(((ea,au,eb,av),(up[0],up[1],vp[0],vp[1]),w))
print('elements normalising the plane:',count_norm,' rank-1 (curve) cases:',len(curves))
print('glued pairs (y, y\') with y\' not W_X-conjugate to y:',len(glued))
pairs={}
for y,yp,w in glued:
    key=(frozenset(wx_images((y[0],y[1]),(y[2],y[3]))), frozenset(wx_images((yp[0],yp[1]),(yp[2],yp[3]))))
    pairs.setdefault(key,[]).append(w)
print('distinct pairs of W_X-orbits glued:',len(pairs))
for (o1,o2),ws in pairs.items():
    r1=min(o1); r2=min(o2)
    print(' orbit1 rep (log_t|u|, arg u, log_t|v|, arg v)=',r1,'  orbit2 rep=',r2,'  #w=',len(ws))
pickle.dump({'glued':glued,'curves':curves},open('e6_points.pkl','wb'))