# Independent check (A9) of the GL_2/T orbit list and Hecke action, using the
# Moebius action of explicit coset representatives on ends in P^1(F_p(t)),
# exact polynomial arithmetic over F_p, and the disc model of the tree.
# No tree is built; distances come from t-adic valuations of differences.
import sys, random
p=int(sys.argv[1]); trials=int(sys.argv[2]); seed=int(sys.argv[3]) if len(sys.argv)>3 else 1
random.seed(seed)
def norm(a):
    a=[c%p for c in a]
    while a and a[-1]==0: a.pop()
    return a
def add(a,b):
    n=max(len(a),len(b)); return norm([(a[i] if i<len(a) else 0)+(b[i] if i<len(b) else 0) for i in range(n)])
def neg(a): return norm([-c for c in a])
def sub(a,b): return add(a,neg(b))
def mul(a,b):
    if not a or not b: return []
    r=[0]*(len(a)+len(b)-1)
    for i,x in enumerate(a):
        for j,y in enumerate(b): r[i+j]+=x*y
    return norm(r)
def ordp(a):
    a=norm(a)
    if not a: return None
    for i,c in enumerate(a):
        if c: return i
def tp(k): return [0]*k+[1]
# an end is a pair (N,D) of polynomials in t (projective point [N:D]); infinity = ([1],[])
INF=([1],[])
def isinf(z): return not norm(z[1])
def ordv(z):  # t-adic valuation of z=N/D (finite z)
    o=ordp(z[0]); return None if o is None else o-ordp(z[1])
def diff_ord(a,c): # ord(a-c), None means a=c (infinite)
    N=sub(mul(a[0],c[1]),mul(c[0],a[1])); o=ordp(N)
    return None if o is None else o-ordp(a[1])-ordp(c[1])
BIG=10**6
def dvert(c,r,a,b):
    # distance from vertex D(c,r) (c finite end, r integer) to geodesic (a,b)
    if isinf(a): a,b=b,a
    if isinf(b):
        al=diff_ord(a,c); al=BIG if al is None else al
        return max(0,r-al)
    al=diff_ord(a,c); be=diff_ord(b,c); e=diff_ord(a,b)
    al=BIG if al is None else al; be=BIG if be is None else be
    a1=min(al,r); b1=min(be,r)
    return (r-max(a1,b1))+(e-min(a1,b1))
ZERO=([],[1])
def inO(z): return (not isinf(z)) and (ordv(z) is None or ordv(z)>=0)
def cls(g):
    a,b=g
    d0=dvert(ZERO,0,a,b); d1=dvert(ZERO,-1,a,b)   # v0=D(0,0)=O, v1=D(0,-1)=t^{-1}O
    if d0==0 and d1==0:
        return 'A+' if (inO(a) and not inO(b)) else 'A-'
    assert abs(d0-d1)==1, (d0,d1)
    return ('0',d0) if d0<d1 else ('1',d1)
def mob(M,z):  # M=[[al,be],[ga,de]] polynomial entries (scalar multiples irrelevant)
    (al,be),(ga,de)=M; N,D=z
    return (add(mul(al,N),mul(be,D)), add(mul(ga,N),mul(de,D)))
# inverses (up to scalar and power of t) of coset reps of IsI/I
def s1inv(c): return [[[],[1]],[[1],[(-c)%p]]]            # (u(c)s1)^{-1} ~ [[0,1],[1,-c]]
def s0inv(c): return [[[(-c)%p,0][:1] if c%p else [], [1]],[tp(2),[]]]  # t*[[-c,t^-1],[t,0]] = [[-ct,1],[t^2,0]]
def s0inv(c):
    return [[norm([0,(-c)%p]),[1]],[tp(2),[]]]
OMinv=[[[],[1]],[tp(1),[]]]   # omega^{-1} ~ omega=[[0,1],[t,0]] up to the central t
def rand_end(V,P):
    if random.random()<0.08: return INF
    lo=random.randint(0,V); # numerator poly of degree<V+P, denominator t^lo  -> valuation >= -lo
    N=norm([random.randrange(p) for _ in range(random.randint(1,V+P))])
    if not N: N=[0]*random.randint(0,3)+[1]
    if random.random()<0.3:   # occasionally a genuinely rational (non-Laurent) end
        D=norm([random.randrange(p) for _ in range(3)]); D=D if D and D[0] else [1]
        D=mul(D,tp(lo)); return (N,D)
    return (N,tp(lo))
q=p
def pred(s,O):  # predicted coefficients of T_s 1_O from A6-2, as dict class->coeff
    def sw(X): return X if X in ('A+','A-') else (('1' if X[0]=='0' else '0'),X[1])
    if s=='s0': return {sw(k):v for k,v in pred('s1',sw(O)).items()}
    if O=='A+': return {'A-':1,('0',0):1}
    if O=='A-': return {'A+':1,('0',0):1}
    i,m=O
    if O==('0',0): return {'A+':q-1,'A-':q-1,('0',0):q-2}
    if i=='1': return {('0',m+1):1}
    return {('1',m-1):q,('0',m):q-1}
def coeff_at(s,O,Oprime): return pred(s,O).get(Oprime,0)
seen={}; bad=0; nchk=0
for _ in range(trials):
    a=rand_end(3,3); b=rand_end(3,3)
    if isinf(a) and isinf(b): continue
    if not isinf(a) and not isinf(b) and diff_ord(a,b) is None: continue
    g=(a,b); Op=cls(g)
    if Op not in ('A+','A-') and Op[1]>3: continue
    seen[Op]=seen.get(Op,0)+1
    for s,inv in (('s1',s1inv),('s0',s0inv)):
        counts={}
        for c in range(q):
            h=(mob(inv(c),a),mob(inv(c),b)); X=cls(h); counts[X]=counts.get(X,0)+1
        # (T_s 1_O)(g) = counts[O]; predicted = coefficient of 1_{Op} in T_s 1_O
        allO=set(counts)|{'A+','A-'}|{(i,m) for i in '01' for m in range(6)}
        for O in allO:
            nchk+=1
            if counts.get(O,0)!=coeff_at(s,O,Op): bad+=1; print('MISMATCH',s,O,Op,counts)
    h=(mob(OMinv,a),mob(OMinv,b)); X=cls(h)
    sw={'A+':'A-','A-':'A+'}.get(Op, None) or (('1' if Op[0]=='0' else '0'),Op[1])
    nchk+=1
    if X!=sw: bad+=1; print('OMEGA MISMATCH',Op,X)
print('p=',p,'classes seen:',sorted(seen.items(),key=str))
print('checks',nchk,'mismatches',bad)