from fractions import Fraction as Fr
from itertools import product
# 1. Tree recursion for GL_2/T with q=3,5: C[T1]phi0 spans; C[T1^2]phi0 misses 1_{d=1}.
def T1(vec,q):
    out={}
    for d,c in vec.items():
        terms={0:2,1:1} if d==0 else ({0:q-1,2:1} if d==1 else {d-1:q,d+1:1})
        for e,t in terms.items(): out[e]=out.get(e,0)+c*t
    return out
def rank(rows,D):
    M=[[Fr(r.get(d,0)) for d in range(D)] for r in rows]; rk=0; cols=list(range(D))
    for c in range(D):
        piv=next((i for i in range(rk,len(M)) if M[i][c]!=0),None)
        if piv is None: continue
        M[rk],M[piv]=M[piv],M[rk]
        for i in range(len(M)):
            if i!=rk and M[i][c]!=0:
                f=M[i][c]/M[rk][c]; M[i]=[a-f*b for a,b in zip(M[i],M[rk])]
        rk+=1
    return rk
for q in (3,5):
    v={0:1}; pw=[v]
    for k in range(16): v=T1(v,q); pw.append(v)
    N=8
    print("q",q,"rank C[T1]phi0 on",N+1,"orbits:",rank(pw[:N+1],N+1))
    even=[pw[2*k] for k in range(N+1)]
    r0=rank(even,2*N+1); r1=rank(even+[{1:1}],2*N+1)
    print("   even span rank",r0,"; adding 1_{d=1} ->",r1,"(1_1 in even span:",r0==r1,")")
    t=T1(T1({0:1},q),q); t[0]-= (q+1); print("   T_(1,-1)phi0 =",t)
# 2. Lagrangian 2-planes in (F_q^4, J) for q=3: expect (q+1)(q^2+1)=40
q=3
J=[[0,1,0,0],[-1,0,0,0],[0,0,0,1],[0,0,-1,0]]
def form(u,v): return sum(u[i]*J[i][j]*v[j] for i in range(4) for j in range(4))%q
vecs=[v for v in product(range(q),repeat=4) if any(v)]
def span(vs):
    S=set()
    for a in range(q):
        for b in range(q): S.add(tuple((a*vs[0][i]+b*vs[1][i])%q for i in range(4)))
    return frozenset(S)
planes=set()
for i,u in enumerate(vecs):
    for v in vecs[i+1:]:
        if form(u,v)==0:
            S=span((u,v))
            if len(S)==q*q: planes.add(S)
print("Lagrangian planes q=3:",len(planes),"expected",(q+1)*(q*q+1))
# 3. T_{(2,0,0,0)}*phi0 at x_{(1,1,1,1)}: count M=O^4+O p^{-2}v with M^vee_J = pM : expect 0 (argued); brute force mod p^2 lattice check
# M determined by primitive v mod p^2 up to units; check condition u in M^vee_J  <=> J(u,v) in p^2 O, versus pM = pO^4 + O p^{-1} v not inside O^4 -> never equal.
print("T_(2,0,0,0) count at (1,1,1,1): 0 (pM not contained in O^4)")