#!/usr/bin/env python3
"""Exact two-variable finite invariant certificates (standard library only).
Characteristic zero uses Q(omega), omega^2+omega+1=0; prime-field tests use
ordinary scalars. This is a finite regression suite, not a general invariant-ring
termination algorithm. Matrices act on coefficient columns of degree-one forms.
"""
from fractions import Fraction as F
from itertools import product
from collections import Counter
from math import comb
from pathlib import Path
import argparse,json
COUNT=Counter()
def require(ok,msg):
 if not ok:raise ValueError(msg)
def check(ok,msg):
 COUNT[msg]+=1
 if not ok:raise RuntimeError(msg)
class Field:
 def __init__(self,p=0):
  require(type(p)is int and (p==0 or (p>=2 and all(p%d for d in range(2,int(p**.5)+1)))),'zero or prime characteristic')
  self.p=p;self.zero=self.c(0);self.one=self.c(1)
 def base(self,a):
  require(type(a)is int or isinstance(a,F),'exact scalar required')
  if not self.p:return F(a)
  require(a.denominator%self.p if isinstance(a,F) else True,'denominator vanishes')
  return (a.numerator*pow(a.denominator,-1,self.p))%self.p if isinstance(a,F)else a%self.p
 def c(self,a,b=0):
  a,b=self.base(a),self.base(b);require(not self.p or b==0,'prime field has no omega coordinate');return(a,b)
 def add(self,x,y):return self.c(x[0]+y[0],x[1]+y[1])
 def neg(self,x):return self.c(-x[0],-x[1])
 def sub(self,x,y):return self.add(x,self.neg(y))
 def mul(self,x,y):return self.c(x[0]*y[0]-x[1]*y[1],x[0]*y[1]+x[1]*y[0]-x[1]*y[1])
 def inv(self,x):
  require(x!=self.zero,'zero scalar inverse')
  if self.p:return self.c(pow(x[0],-1,self.p))
  n=x[0]*x[0]-x[0]*x[1]+x[1]*x[1];return self.c((x[0]-x[1])/n,-x[1]/n)
 def div(self,x,y):return self.mul(x,self.inv(y))
 def pow(self,x,n):
  require(type(n)is int and n>=0,'nonnegative exponent');out=self.one
  while n:
   if n&1:out=self.mul(out,x)
   x=self.mul(x,x);n//=2
  return out
 def sum(self,x):
  out=self.zero
  for a in x:out=self.add(out,a)
  return out
 def show(self,a):return str(a[0]) if a[1]==0 else str(a[0])+' + ('+str(a[1])+')*omega'
def eye(f,n):return [[f.one if i==j else f.zero for j in range(n)]for i in range(n)]
def mm(f,A,B):return [[f.sum(f.mul(a,b)for a,b in zip(row,col))for col in zip(*B)]for row in A]
def submat(f,A,B):return [[f.sub(a,b)for a,b in zip(x,y)]for x,y in zip(A,B)]
def rref(f,A):
 A=[row[:]for row in A];piv=[];r=0
 for j in range(len(A[0])if A else 0):
  i=next((i for i in range(r,len(A))if A[i][j]!=f.zero),None)
  if i is None:continue
  A[i],A[r]=A[r],A[i];s=f.inv(A[r][j]);A[r]=[f.mul(s,x)for x in A[r]]
  for i in range(len(A)):
   if i!=r:
    s=A[i][j];A[i]=[f.sub(x,f.mul(s,y))for x,y in zip(A[i],A[r])]
  piv.append(j);r+=1
  if r==len(A):break
 return A,piv
def kernel(f,A):
 A,piv=rref(f,A);n=len(A[0])if A else 0;out=[]
 for j in range(n):
  if j not in piv:
   v=[f.zero]*n;v[j]=f.one
   for i,k in enumerate(piv):v[k]=f.neg(A[i][j])
   out.append(v)
 return out
def mkey(A):return tuple(map(tuple,A))
def invert2(f,A):
 d=f.sub(f.mul(A[0][0],A[1][1]),f.mul(A[0][1],A[1][0]));s=f.inv(d)
 return [[f.mul(s,A[1][1]),f.mul(s,f.neg(A[0][1]))],[f.mul(s,f.neg(A[1][0])),f.mul(s,A[0][0])]]
def valid_group(f,G):
 require(bool(G)and all(len(A)==2 and all(len(r)==2 for r in A)for A in G),'nonempty group of two by two matrices')
 keys={mkey(A)for A in G};require(len(keys)==len(G)and mkey(eye(f,2))in keys,'distinct group elements including identity')
 for A in G:
  invert2(f,A)
  for B in G:require(mkey(mm(f,A,B))in keys,'group must be closed')
 return G
def cyclic(f,A,n):
 G=[eye(f,2)]
 for _ in range(1,n):G.append(mm(f,G[-1],A))
 require(mm(f,G[-1],A)==G[0],'declared order');return valid_group(f,G)
def clean(f,p):return {a:c for a,c in p.items()if c!=f.zero}
def psum(f,p,q):
 out=dict(p)
 for a,c in q.items():out[a]=f.add(out.get(a,f.zero),c)
 return clean(f,out)
def pscale(f,p,c):return clean(f,{a:f.mul(v,c)for a,v in p.items()})
def pmul(f,p,q):
 out={}
 for a,x in p.items():
  for b,y in q.items():
   e=(a[0]+b[0],a[1]+b[1]);out[e]=f.add(out.get(e,f.zero),f.mul(x,y))
 return clean(f,out)
def ppow(f,p,n):
 out={(0,0):f.one}
 for _ in range(n):out=pmul(f,out,p)
 return out
def act(f,A,p):
 # Column j gives the image of variable j; x -> A00*x+A10*y.
 x=clean(f,{(1,0):A[0][0],(0,1):A[1][0]});y=clean(f,{(1,0):A[0][1],(0,1):A[1][1]});out={}
 for (i,j),v in p.items():out=psum(f,out,pscale(f,pmul(f,ppow(f,x,i),ppow(f,y,j)),v))
 return out
def average(f,G,p):
 inv=f.inv(f.c(len(G)));out={}
 for A in G:out=psum(f,out,act(f,A,p))
 return pscale(f,out,inv)
def action_matrix(f,A,d):
 mon=[(i,d-i)for i in range(d+1)];cols=[act(f,A,{e:f.one})for e in mon]
 return [[p.get(e,f.zero)for p in cols]for e in mon]
def basis_degree(f,G,d):
 require(type(d)is int and d>=0,'nonnegative degree');D=d+1;M=[action_matrix(f,A,d)for A in G];I=eye(f,D)
 fixed=kernel(f,[row for A in M for row in submat(f,A,I)]);result={'dimension':len(fixed),'fixed_basis':fixed}
 if f.c(len(G))!=f.zero:
  scale=f.inv(f.c(len(G)));P=[[f.mul(scale,f.sum(A[i][j]for A in M))for j in range(D)]for i in range(D)]
  check(mm(f,P,P)==P,'average matrix idempotence');_,piv=rref(f,P);cols=[[P[i][j]for i in range(D)]for j in piv]
  check(len(cols)==len(fixed),'average and joint kernel dimension')
  check(all(mm(f,A,[[x]for x in v])==[[x]for x in v]for A in M for v in cols),'average columns fixed')
  check(all(mm(f,P,[[x]for x in v])==[[x]for x in v]for v in fixed),'all joint kernel vectors retained')
  result.update(projection=P,average_basis=cols)
 return result
def present_c3(i,j):
 require(type(i)is int and type(j)is int and i>=0 and j>=0 and (i-j)%3==0,'invariant monomial required')
 r=i%3;return((i-r)//3,(j-r)//3,r)
def expand_c3(a,b,c):return(3*a+c,3*b+c)
def normalize_relation(a,b,c):
 require(all(type(v)is int and v>=0 for v in [a,b,c]),'nonnegative monomial exponents');q,r=divmod(c,3);return a+q,b+q,r
def convolution(a,b):
 c=[F(0)]*(len(a)+len(b)-1)
 for i,x in enumerate(a):
  for j,y in enumerate(b):c[i+j]+=x*y
 while len(c)>1 and c[-1]==0:c.pop()
 return c
def polyadd(a,b):
 c=[(a[i]if i<len(a)else 0)+(b[i]if i<len(b)else 0)for i in range(max(len(a),len(b)))]
 while len(c)>1 and c[-1]==0:c.pop()
 return c
def ratio_coeff(a,b,N):
 require(b[0]!=0,'nonzero series constant');out=[]
 for n in range(N+1):out.append(((a[n]if n<len(a)else 0)-sum(b[j]*out[n-j]for j in range(1,min(n,len(b)-1)+1)))/b[0])
 return out
def showpoly(f,p):return [{'exponents':list(e),'coefficient':f.show(v)}for e,v in sorted(p.items())]
def rejected(fn):
 try:fn()
 except ValueError:return True
 return False

def run():
 COUNT.clear();f=Field();w=f.c(0,1);diag=[[w,f.zero],[f.zero,f.pow(w,2)]];G=cyclic(f,diag,3);records={}
 U={(3,0):f.one};V={(0,3):f.one};W={(1,1):f.one}
 for d in range(19):
  r=basis_degree(f,G,d);expected=[(i,d-i)for i in range(d+1)if(2*i-d)%3==0]
  check(r['dimension']==len(expected),'C3 fixed monomial count')
  for i in range(d+1):
   p={(i,d-i):f.one};check(average(f,G,p)==(p if(i,d-i)in expected else{}),'cyclotomic average selector')
  if d==6:records['C3_degree6_basis']=[[f.show(x)for x in v]for v in r['average_basis']]
 for i,j in product(range(31),repeat=2):
  if(i-j)%3:continue
  a,b,c=present_c3(i,j);check(expand_c3(a,b,c)==(i,j)and c<3,'C3 invariant normal form')
  if i>=j:check(expand_c3((i-j)//3,0,j)==(i,j),'one-sided normal form upper')
  else:check(expand_c3(0,(j-i)//3,i)==(i,j),'one-sided normal form lower')
 for a,b,c in product(range(9),repeat=3):check(expand_c3(*normalize_relation(a,b,c))==expand_c3(a,b,c),'formal relation reduction preserves image')
 for e in product(range(6),repeat=2):
  p={e:f.c(2,1),(e[0]+1,e[1]):f.one}
  for inv in [U,V,W,psum(f,U,W)]:check(average(f,G,pmul(f,inv,p))==pmul(f,inv,average(f,G,p)),'invariant module linearity')
 sixth=ppow(f,{(1,0):f.one,(0,1):f.one},6);avg=average(f,G,sixth)
 check(avg=={(6,0):f.one,(3,3):f.c(20),(0,6):f.one},'sixth binomial average');records['averaged_sixth']=showpoly(f,avg)
 # Nondiagonal and changed-basis rational order-three actions exercise mixed columns.
 A=[[f.zero,f.c(-1)],[f.one,f.c(-1)]];GA=cyclic(f,A,3)
 for T in [[[1,1],[0,1]],[[2,1],[1,1]],[[1,0],[2,1]],[[1,-1],[1,1]]]:
  T=[[f.c(x)for x in r]for r in T];H=valid_group(f,[mm(f,mm(f,T,B),invert2(f,T))for B in GA])
  for d in range(13):check(basis_degree(f,H,d)['dimension']==sum((2*i-d)%3==0 for i in range(d+1)),'nondiagonal basis-covariant dimensions')
 records['rational_order3_degree3_basis']=[[f.show(x)for x in v]for v in basis_degree(f,GA,3)['average_basis']]
 # Two independently derived rational Hilbert expressions, checked by polynomial cross multiplication.
 d1=list(map(F,[1,-2,1]));d2=list(map(F,[1,1,1]));mol_num=[v/3 for v in polyadd(d2,[2*v for v in d1])];mol_den=convolution(d1,d2)
 normal_num=list(map(F,[1,0,1,0,1]));normal_den=convolution(list(map(F,[1,0,0,-1])),list(map(F,[1,0,0,-1])))
 rel_num=list(map(F,[1,0,0,0,0,0,-1]));rel_den=convolution(normal_den,list(map(F,[1,0,-1])))
 check(convolution(mol_num,normal_den)==convolution(normal_num,mol_den),'Molien versus module rational identity')
 check(convolution(rel_num,normal_den)==convolution(normal_num,rel_den),'relation versus module rational identity')
 hs=ratio_coeff(mol_num,mol_den,120);hn=ratio_coeff(normal_num,normal_den,120)
 for d in range(121):
  direct=sum((2*i-d)%3==0 for i in range(d+1));normal=sum(3*(a+b)+2*r==d for r in range(3)for a in range(d//3+1)for b in range(d//3+1));q,r=divmod(d,3)
  check(hs[d]==hn[d]==direct==normal==([q+1,q,q+1][r]),'four-way all tested Hilbert coefficients')
 records['Hilbert_first_19']=[int(x)for x in hs[:19]]
 # Characteristic two shear: no averaging, but independent fixed kernels equal k[y,N].
 f2=Field(2);H=cyclic(f2,[[f2.one,f2.zero],[f2.one,f2.one]],2);N={(2,0):f2.one,(1,1):f2.one};y={(0,1):f2.one}
 check(act(f2,H[1],N)==N,'modular norm polynomial invariant')
 for d in range(25):
  r=basis_degree(f2,H,d);basis=[pmul(f2,ppow(f2,y,d-2*j),ppow(f2,N,j))for j in range(d//2+1)];M=[[p.get((i,d-i),f2.zero)for p in basis]for i in range(d+1)]
  check(len(rref(f2,M)[1])==len(basis)==r['dimension'],'modular generators complete degree piece')
  for p in basis:check(act(f2,H[1],p)==p,'modular basis fixed')
 check(rejected(lambda:average(f2,H,{(0,0):f2.one})),'modular averaging rejected')
 transfer=lambda p:psum(f2,p,act(f2,H[1],p))
 check(transfer({(0,0):f2.one})=={}and transfer({(1,0):f2.one})==y and transfer(y)=={},'transfer not projection')
 check(psum(f2,ppow(f2,{(1,0):f2.one},4),{(2,2):f2.one})==ppow(f2,N,2),'modular fourth-degree example')
 # Nonmodular positive characteristic still has projection; dimensions are computed as integer ranks.
 for p,root in [(7,2),(13,3),(19,7)]:
  k=Field(p);a=k.c(root);J=cyclic(k,[[a,k.zero],[k.zero,k.pow(a,2)]],3)
  for d in range(19):check(basis_degree(k,J,d)['dimension']==sum((2*i-d)%3==0 for i in range(d+1)),'prime-field integer dimension')
  # All p^2 rational points: U,V,W fibres are exactly group orbits for this split action.
  buckets={}
  for x,z in product(range(p),repeat=2):buckets.setdefault((pow(x,3,p),pow(z,3,p),x*z%p),set()).add((x,z))
  for bucket in buckets.values():
   x,z=next(iter(bucket));orbit={(pow(root,j,p)*x%p,pow(root,2*j,p)*z%p)for j in range(3)};check(bucket==orbit,'split prime-field orbit fibres')
 records['modular_ring']='F2[y, x^2+x*y]; averaging denominator rejected; degree d dimension floor(d/2)+1'
 # Failure boundaries and finite coefficients do not imply full equality.
 C2=cyclic(f,[[f.c(-1),f.zero],[f.zero,f.one]],2);x={(1,0):f.one}
 check(average(f,C2,x)=={}and average(f,C2,pmul(f,x,x))==pmul(f,x,x),'average is not multiplicative')
 fake=hs[:];fake[100]+=1;check(fake[:20]==hs[:20]and fake!=hs,'finite prefix is not series equality')
 for fn in [lambda:Field(4),lambda:f.c(.5),lambda:f.inv(f.zero),lambda:present_c3(1,0),lambda:valid_group(f,GA[:-1]),lambda:valid_group(f,[GA[0],GA[0]]),lambda:basis_degree(f,G,-1)]:check(rejected(fn),'invalid input rejected')
 return {'status':'PASS','checks':sum(COUNT.values()),'groups':dict(COUNT),'records':records,'scope':'Exact two-variable finite certificates; arbitrary-characteristic finite generation and full presentations are proved in the articles, not inferred from the tested degree ranges.'}
if __name__=='__main__':
 ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path);args=ap.parse_args();out=json.dumps(run(),ensure_ascii=False,indent=2)+'\n'
 if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(out);print('PASS',sum(COUNT.values()),'checks',args.output)
 else:print(out,end='')