#!/usr/bin/env python3
"""Exact finite certificates for local curves. Python standard library only.
No filesystem writes unless --output is supplied; checks remain active under -O.
Finite checks supplement the articles' general proofs, not normalization of an
arbitrary input curve. Coefficients are rational or reduced modulo a stated prime.
"""
from fractions import Fraction as Q
from itertools import product
from math import gcd
from collections import Counter
from pathlib import Path
import argparse,json
CHECKS=0;GROUPS=Counter()
def check(ok,label):
 global CHECKS
 CHECKS+=1;GROUPS[label]+=1
 if not ok:raise RuntimeError(label)
def add(p,q):
 r=p.copy()
 for m,c in q.items():r[m]=r.get(m,Q(0))+c
 return {m:c for m,c in r.items() if c}
def scale(p,a):return {m:c*a for m,c in p.items() if c*a}
def mul(p,q):
 r={}
 for (i,j),c in p.items():
  for (k,l),d in q.items():r[i+k,j+l]=r.get((i+k,j+l),Q(0))+c*d
 return {m:c for m,c in r.items() if c}
def mon(i,j,c=1):return {(i,j):Q(c)} if c else {}
def power(p,n):
 r=mon(0,0)
 for _ in range(n):r=mul(r,p)
 return r
def derivative(p,axis):
 r={}
 for m,c in p.items():
  if m[axis]:
   n=list(m);n[axis]-=1;r[tuple(n)]=c*m[axis]
 return r
def normal(p,G):
 p=p.copy();r={}
 while p:
  i,j=max(p);c=p[i,j]
  for g in G:
   a,b=max(g)
   if i>=a and j>=b:
    p=add(p,scale(mul(mon(i-a,j-b),g),-c/g[a,b]));break
  else:r[i,j]=c;del p[i,j]
 return r
def rank(columns,rows,prime=0):
 if not columns:return 0
 A=[[columns[j][i] for j in range(len(columns))] for i in range(rows)]
 if prime:A=[[int(v)%prime for v in row] for row in A]
 else:A=[[Q(v) for v in row] for row in A]
 r=0
 for j in range(len(columns)):
  k=next((k for k in range(r,rows) if A[k][j]),None)
  if k is None:continue
  A[r],A[k]=A[k],A[r];inv=pow(A[r][j],-1,prime) if prime else 1/A[r][j]
  A[r]=[(v*inv)%prime if prime else v*inv for v in A[r]]
  for k in range(rows):
   if k==r or not A[k][j]:continue
   c=A[k][j];A[k]=[(v-c*w)%prime if prime else v-c*w for v,w in zip(A[k],A[r])]
  r+=1
  if r==rows:break
 return r
def matvec(M,v):return [sum(a*b for a,b in zip(row,v)) for row in M]
def matmul(A,B):return [[sum(A[i][k]*B[k][j] for k in range(len(B))) for j in range(len(B[0]))] for i in range(len(A))]
def ident(n):return [[Q(i==j) for j in range(n)] for i in range(n)]
def matpow(M,n):
 R=ident(len(M))
 for _ in range(n):R=matmul(R,M)
 return R

# The extra lowest relation and an actual finite-dimensional representation.
f=add(mon(2,0),mon(0,3));g=mon(1,1);h=mon(0,4);G=[f,g,h]
check(add(mul(mon(0,1),f),scale(mul(mon(1,0),g),-1))==h,'extra_y4_relation')
for i in range(3):
 for j in range(i):
  a=max(G[i]);b=max(G[j]);l=(max(a[0],b[0]),max(a[1],b[1]))
  s=add(mul(mon(l[0]-a[0],l[1]-a[1]),G[i]),scale(mul(mon(l[0]-b[0],l[1]-b[1]),G[j]),-1))
  check(not normal(s,G),'Buchberger_S_pair')
basis=[(0,0),(1,0),(0,1),(0,2),(0,3)];pos={v:i for i,v in enumerate(basis)}
Mx=[[Q(0)]*5 for _ in range(5)];My=[[Q(0)]*5 for _ in range(5)]
Mx[1][0]=1;Mx[4][1]=-1;My[2][0]=1;My[3][2]=1;My[4][3]=1
zero=[[Q(0)]*5 for _ in range(5)]
check(matmul(Mx,My)==matmul(My,Mx)==zero,'commuting_five_dimensional_representation')
x2=matpow(Mx,2);y3=matpow(My,3)
check(all(x2[i][j]+y3[i][j]==0 for i in range(5) for j in range(5)),'representation_defining_relation')
e0=[Q(1),Q(0),Q(0),Q(0),Q(0)]
for i in range(9):
 for j in range(9):
  rem=normal(mon(i,j),G);v=[Q(0)]*5
  for b,c in rem.items():v[pos[b]]=c
  check(v==matvec(matmul(matpow(Mx,i),matpow(My,j)),e0),'normal_form_vs_matrix')
filtered=[]
for r in range(7):
 cols=[]
 for i in range(r+1):
  M=matmul(matpow(Mx,i),matpow(My,r-i));cols.extend([list(v) for v in zip(*M)])
 filtered.append(rank(cols,5))
check(filtered==[5,4,2,1,0,0,0],'maximal_ideal_filtration')
check([filtered[i]-filtered[i+1] for i in range(6)]==[1,2,1,1,0,0],'graded_dimensions')
for i,j in product(range(12),repeat=2):
 survives=not(i>=2 or (i>=1 and j>=1) or j>=4)
 check(survives==((i,j) in basis),'complete_initial_ideal_standard_basis')

# Local jets calculated by a finite ideal-span matrix, not by selecting the
# lowest polynomial alone. Higher terms remain in every truncated relation.
curves={'parabola':add(mon(0,1),mon(2,0,-1)),
        'cusp':add(mon(0,2),mon(3,0,-1)),
        'node':add(add(mon(0,2),mon(2,0,-1)),mon(3,0,-1)),
        'three_four':add(mon(0,3),mon(4,0,-1))}
jets={}
for name,F in curves.items():
 order=min(sum(m) for m in F);dims=[]
 for r in range(13):
  mons=[(i,j) for i in range(r+1) for j in range(r+1-i)];index={m:i for i,m in enumerate(mons)};cols=[]
  for i in range(max(0,r-order+1)):
   for j in range(r-order-i+1):
    v=[Q(0)]*len(mons)
    for m,c in mul(mon(i,j),F).items():
     if sum(m)<=r:v[index[m]]=c
    cols.append(v)
  d=len(mons)-rank(cols,len(mons));dims.append(d)
  check(d==sum(min(order,n+1) for n in range(r+1)),'local_jet_dimension')
 jets[name]=dims
 dx=derivative(F,0).get((0,0),0);dy=derivative(F,1).get((0,0),0)
 check((dx,dy)==((0,1) if name=='parabola' else (0,0)),'origin_Jacobian')

# The two quadratic fiber rings, and the general q=t^2-c family. Each subring
# is k+q k[t]. A polynomial is in it iff its q-remainder has no linear term.
def remainder_q(v,c):
 return (sum(a*c**(i//2) for i,a in enumerate(v) if i%2==0),sum(a*c**(i//2) for i,a in enumerate(v) if i%2))
for c in map(Q,[-1,0,1,2,4]):
 M=[[Q(0),c],[Q(1),Q(0)]]
 check(matmul(M,M)==[[c,Q(0)],[Q(0),c]],'fiber_multiplication')
 for coeff in product([-1,0,1],repeat=7):
  a,b=remainder_q(coeff,c);at,bt=remainder_q([0,*coeff],c)
  check((b==0 and bt==0)==(a==0 and b==0),'conductor_complete_membership')
  check((at,bt)==(b*c,a),'fiber_remainder_times_t')
 for i in range(30):
  # t^i q is an A-combination of x and y, via t^2=x+c.
  rem0,rem1=remainder_q([0]*i+[-c,0,1],c)
  check(rem0==rem1==0,'conductor_generators_all_degrees')
check(matmul([[0,0],[1,0]],[[0,0],[1,0]])==[[0,0],[0,0]],'cusp_nilpotent_fiber')
for sign in [-1,1]:
 e=[Q(1,2),Q(sign,2)]
 prod=[e[0]**2+e[1]**2,2*e[0]*e[1]]
 check(prod==e,'node_fiber_idempotents')

# Numerical semigroup membership is computed from generators. A consecutive
# run of a members certifies the entire remaining tail by adding a.
def semigroup(a,b,limit):
 S={0}
 for n in range(limit+1):
  if n in S:
   if n+a<=limit:S.add(n+a)
   if n+b<=limit:S.add(n+b)
 return S
def conductor(a,b):
 limit=3*a*b+2*a+2*b;S=semigroup(a,b,limit)
 c=next(c for c in range(limit-a+1) if all(c+j in S for j in range(a)))
 return c,S
semigroup_rows=[]
for a in range(2,10):
 for b in range(a+1,16):
  if gcd(a,b)!=1:continue
  c,S=conductor(a,b);gaps=[n for n in range(c) if n not in S]
  check(c-1 not in S and all(c+j in S for j in range(a)),'conductor_tail_certificate')
  for d in sorted(s for s in S if s<c):check(d+(c-1-d) not in S,'conductor_minimality_witness')
  for n in range(c,min(max(S)+1,c+4*a*b)):
   check(n in S,'semigroup_tail_membership')
   check(any(n-c-j in S for j in range(a)),'conductor_A_generators')
  # Actual homogeneous module matrices for b x^{b-1} / a y^{a-1} relation.
  # Free basis coefficient exponents are w-a,w-b; the equation has degree ab.
  degree_kernel={}
  for w in range(a*b+c+a+b+1):
   free=[v for v in [a,b] if w-v in S]
   image=[a if v==a else b for v in free]
   rel=[-b if v==a else a for v in free] if w-a*b in S else []
   rr=rank([rel],len(free)) if rel else 0;ri=rank([image],len(free)) if free else 0
   if rel:check(sum(x*y for x,y in zip(rel,image))==0,'differential_relation_maps_zero')
   k=len(free)-rr-ri;check(k>=0,'homogeneous_differential_kernel_rank')
   if k:degree_kernel[w]=k
   if w>=a*b+c:check(k==0,'differential_finite_tail_certificate')
  # This is a finite comparison for tested two-generator curves, not a claim
  # about arbitrary singular curves.
  check(sum(degree_kernel.values())==2*len(gaps),'tested_two_generator_torsion_dimension')
  semigroup_rows.append({'a':a,'b':b,'conductor':c,'delta':len(gaps),'torsion_dimension':sum(degree_kernel.values())})

# Prove the displayed cyclic generators span every homogeneous kernel in the
# two assigned cases, and compute their annihilator via monomial ideal membership.
assigned=[]
for a,b,ann,expected in [(2,3,[3,4],2),(3,4,[8,9],6)]:
 c,S=conductor(a,b);frobenius=a*b-a-b;weights=[];ann_basis=[];coker=[]
 for e in sorted(S):
  if e>a*b+c:continue
  ideal=any(e-g in S for g in ann);h_exists=e-frobenius in S
  check(ideal==h_exists,'tau_annihilator_ideal')
  if not ideal:ann_basis.append(e)
 for w in range(a*b+c+a+b+1):
  d=int(w-a in S)+int(w-b in S);rel=int(w-a*b in S);image=int(d>0)
  kernel=d-rel-image;tau=int(w-a-b in S and w-a*b not in S)
  check(kernel==tau,'tau_spans_complete_kernel')
  if kernel:weights.append(w)
  if w>=1 and not image:coker.append(w-1)
 check(len(ann_basis)==len(weights)==expected,'assigned_torsion_length')
 check(len(coker)==(1 if a==2 else 3),'assigned_differential_cokernel')
 check(ann_basis==([0,2] if a==2 else [0,3,4,6,7,10]),'assigned_annihilator_basis')
 # Recompute ranks in characteristics dividing a or b. No char-zero division.
 positive_characteristic={}
 for prime in [2,3,5,7]:
  length=0
  for w in range(a*b+c+a+b+1):
   free=[v for v in [a,b] if w-v in S];image=[a if v==a else b for v in free]
   rel=[-b if v==a else a for v in free] if w-a*b in S else []
   rr=rank([rel],len(free),prime) if rel else 0;ri=rank([image],len(free),prime) if free else 0
   length+=len(free)-rr-ri
  positive_characteristic[prime]=length
 check(positive_characteristic==({2:4,3:3,5:2,7:2} if a==2 else {2:8,3:9,5:6,7:6}),'positive_characteristic_torsion_changes')
 assigned.append({'generators':[a,b],'conductor':c,'gaps':[n for n in range(c) if n not in S],'torsion_weights':weights,'annihilator_quotient_exponents':ann_basis,'cokernel_dt_exponents':coker,'torsion_length_by_characteristic':positive_characteristic})
# Distinguish the conductor quotient from the fiber for the migrated curve.
check(conductor(3,4)[0]==6 and min(3,4)==3,'conductor_quotient_not_fiber')
result={'status':'PASS','checks':CHECKS,'groups':dict(GROUPS),'arithmetic':'exact integers/Fraction, plus explicitly named finite characteristics','scope':{'quadratic_subring_polynomial_tests':5*3**7,'two_generator_semigroups':len(semigroup_rows),'local_jet_orders':list(range(13))},'artinian_counterexample':{'basis':['1','x','y','y^2','y^3'],'ideal_power_dimensions':filtered,'associated_graded_relation':['X^2','XY','Y^4']},'local_jet_lengths':jets,'assigned_curves':assigned,'tested_semigroups':semigroup_rows,'limitations':['Finite algebra certificates do not prove normalization of arbitrary input curves.','General universal properties, support equivalence and localization are proved in prose.','The two-generator torsion comparisons are finite tests; no identity for all curve singularities is inferred.']}
ap=argparse.ArgumentParser();ap.add_argument('--output',type=Path);args=ap.parse_args();s=json.dumps(result,ensure_ascii=False,indent=2)+'\n'
if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(s)
print(s,end='')