#!/usr/bin/env python3
"""Exact certificates for positive binary forms and fundamental discriminants.

Only Python's standard library. Reduce accepts nonprimitive positive forms;
compose accepts primitive positive forms of the same negative fundamental D.
Ideal products are checked independently using two-dimensional integer lattices.
The finite tests audit implementations; the accompanying pages prove the claims.
"""
from pathlib import Path
from fractions import Fraction as F
from math import gcd,isqrt
from itertools import product
from collections import Counter
import argparse,json,random
I=(1,0,0,1);S=(0,-1,1,0)
def ints(v,length):
 if not isinstance(v,(list,tuple)) or len(v)!=length or any(type(x) is not int for x in v):raise ValueError('integer tuple required')
 return tuple(v)
def disc(q):a,b,c=q;return b*b-4*a*c
def valid(q,primitive=False):
 q=ints(q,3)
 if q[0]<=0 or disc(q)>=0:raise ValueError('positive definite form required')
 if primitive and gcd(*q)!=1:raise ValueError('primitive form required')
 return q
def det(u):p,q,r,s=u;return p*s-q*r
def mm(u,v):p,q,r,s=u;a,b,c,d=v;return(p*a+q*c,p*b+q*d,r*a+s*c,r*b+s*d)
def inv(u):
 if det(u)!=1:raise ValueError('SL2 matrix required')
 p,q,r,s=u;return(s,-q,-r,p)
def change(q,u):
 a,b,c=ints(q,3);p,t,r,s=ints(u,4)
 return(a*p*p+b*p*r+c*r*r,2*a*p*t+b*(p*s+t*r)+2*c*r*s,a*t*t+b*t*s+c*s*s)
def val(q,x,y):a,b,c=q;return a*x*x+b*x*y+c*y*y
def reduced(q):
 a,b,c=valid(q)
 return abs(b)<=a<=c and (b>=0 or (abs(b)!=a and a!=c))
def reduce_form(q):
 start=valid(q);q=start;u=I;steps=[]
 while True:
  a,b,c=q;k=(a-b)//(2*a)
  if k:
   v=(1,k,0,1);new=change(q,v);steps.append({'op':'T','k':k,'before':q,'after':new,'matrix':v});q=new;u=mm(u,v)
  a,b,c=q
  if a>c:
   new=change(q,S);steps.append({'op':'S','before':q,'after':new,'matrix':S});q=new;u=mm(u,S)
  else:
   if a==c and b<0:
    new=change(q,S);steps.append({'op':'boundary-S','before':q,'after':new,'matrix':S});q=new;u=mm(u,S)
   break
 return {'input':start,'reduced':q,'matrix':u,'steps':steps}
def verify_reduction(q,certificate):
 try:
  q=valid(q);u=ints(certificate['matrix'],4);out=valid(certificate['reduced'])
  if ints(certificate['input'],3)!=q:return False
  if det(u)!=1 or change(q,u)!=out or not reduced(out):return False
  acc=I;cur=q
  for step in certificate['steps']:
   v=ints(step['matrix'],4)
   if det(v)!=1 or tuple(step['before'])!=cur or change(cur,v)!=tuple(step['after']):return False
   if step['op']=='T' and v!=(1,step['k'],0,1):return False
   if step['op'] in ('S','boundary-S') and v!=S:return False
   if step['op'] not in ('T','S','boundary-S'):return False
   cur=change(cur,v);acc=mm(acc,v)
  return acc==u and cur==out
 except (ValueError,KeyError,TypeError):return False
def squarefree(n):
 n=abs(n)
 if n==0:return False
 p=2
 while p*p<=n:
  if n%(p*p)==0:return False
  p+=1
 return True
def fundamental(d):
 return type(d) is int and d<0 and ((d%4==1 and squarefree(d)) or (d%4==0 and d//4%4 in (2,3) and squarefree(d//4)))
def enumerate_reduced(d,primitive=True):
 if type(d) is not int or d>=0 or d%4 not in (0,1):raise ValueError('negative quadratic discriminant required')
 out=[]
 for a in range(1,isqrt(abs(d)//3)+1):
  for b in range(-a,a+1):
   if (b*b-d)%(4*a):continue
   q=(a,b,(b*b-d)//(4*a))
   if reduced(q) and (not primitive or gcd(*q)==1):out.append(q)
 return sorted(out)
def egcd(a,b):
 aa,bb=abs(a),abs(b);oldr,r=aa,bb;olds,s=1,0;oldt,t=0,1
 while r:
  k=oldr//r;oldr,r=r,oldr-k*r;olds,s=s,olds-k*s;oldt,t=t,oldt-k*t
 return oldr,olds*(1 if a>=0 else -1),oldt*(1 if b>=0 else -1)
def first_column(x,y):
 g,u,v=egcd(x,y)
 if g!=1:raise ValueError('primitive vector required')
 return (x,-v,y,u)
def coprime_representative(q,N):
 q=valid(q,True)
 if type(N) is not int or N<1:raise ValueError('positive modulus required')
 for x,y in [(1,0),(0,1)]:
  if gcd(val(q,x,y),N)==1:return change(q,first_column(x,y)),first_column(x,y)
 h=1
 while True:
  for x in range(-h,h+1):
   for y in range(1,h+1):
    if max(abs(x),y)!=h or gcd(x,y)!=1:continue
    if gcd(val(q,x,y),N)==1:
     u=first_column(x,y);return change(q,u),u
  h+=1

def compose(q1,q2):
 q1=valid(q1,True);q2=valid(q2,True);d=disc(q1)
 if disc(q2)!=d or not fundamental(d):raise ValueError('same negative fundamental discriminant required')
 q2p,u2=coprime_representative(q2,q1[0]);a,b,_=q1;ap,bp,_=q2p
 # B=b+2*a*t; a*t=(bp-b)/2 mod ap.
 if ap==1:t=0
 else:t=(((bp-b)//2)*pow(a,-1,ap))%ap
 B=b+2*a*t;A=a*ap;k=(A-B)//(2*A);B+=2*A*k # canonical interval (-A,A]
 if (B*B-d)%(4*A):raise AssertionError('CRT integrality')
 raw=(A,B,(B*B-d)//(4*A));red=reduce_form(raw)
 return {'inputs':[q1,q2],'second_input_matrix':u2,'second_representative':q2p,'common_B':B,'common_forms':[change(q1,(1,(B-b)//(2*a),0,1)),change(q2p,(1,(B-bp)//(2*ap),0,1))],'raw_product':raw,'reduction':red,'reduced':red['reduced']}
def verify_composition(q1,q2,cert):
 try:
  q1=valid(q1,True);q2=valid(q2,True);d=disc(q1)
  if disc(q2)!=d or not fundamental(d):return False
  if [ints(x,3) for x in cert['inputs']]!=[q1,q2]:return False
  u=ints(cert['second_input_matrix'],4);qp=ints(cert['second_representative'],3)
  if det(u)!=1 or change(q2,u)!=qp:return False
  a,b,_=q1;ap,bp,_=qp;B=cert['common_B']
  if type(B) is not int or gcd(a,ap)!=1 or (B-b)%(2*a) or (B-bp)%(2*ap):return False
  A=a*ap
  if (B*B-d)%(4*A):return False
  raw=(A,B,(B*B-d)//(4*A))
  common=[change(q1,(1,(B-b)//(2*a),0,1)),change(qp,(1,(B-bp)//(2*ap),0,1))]
  return [tuple(x) for x in cert['common_forms']]==common and tuple(cert['raw_product'])==raw and verify_reduction(raw,cert['reduction']) and tuple(cert['reduced'])==tuple(cert['reduction']['reduced'])
 except (ValueError,KeyError,TypeError):return False

def constants(d):return d%2,((d%2)**2-d)//4
def field_mul(z,w,d):
 delta,t=constants(d);a,b=z;c,e=w
 return(a*c-t*b*e,a*e+b*c+delta*b*e)
def norm(z,d):a,b=z;delta,t=constants(d);return a*a+delta*a*b+t*b*b
def canonical_ideal(q):a,b,c=q;return((a,0),((b-disc(q)%2)//2,1))
def lattice_hnf(vectors):
 # Column lattice: construct positive y gcd and then horizontal kernel gcd.
 g=0;r=0
 for x,y in vectors:
  h,u,v=egcd(g,y);r=u*r+v*x;g=h
 if not g:raise ValueError('rank two lattice required')
 h=0
 for x,y in vectors:h=gcd(h,x-(y//g)*r)
 if not h:raise ValueError('rank two lattice required')
 return((h,0),(r%h,g))
def ideal_form(basis,d):
 if not fundamental(d):raise ValueError('negative fundamental discriminant required')
 (p,r),(q,s)=basis;p,r,q,s=ints((p,r,q,s),4);n=p*s-q*r
 if n<=0:raise ValueError('positive oriented basis required')
 def in_lattice(z):
  x,y=z;return (s*x-q*y)%n==0 and (-r*x+p*y)%n==0
 if not all(in_lattice(field_mul(z,(0,1),d)) for z in basis):raise ValueError('lattice not an O-ideal')
 nums=(norm((p,r),d),norm((p+q,r+s),d)-norm((p,r),d)-norm((q,s),d),norm((q,s),d))
 if any(x%n for x in nums):raise AssertionError('nonintegral normalized norm')
 f=tuple(x//n for x in nums)
 if disc(f)!=d or gcd(*f)!=1:raise AssertionError('form discriminant/content')
 lam=(F(p,f[0]),F(r,f[0]));rebuilt=tuple(field_mul(lam,z,d) for z in canonical_ideal(f))
 if rebuilt!=basis:raise AssertionError('inverse scalar failed')
 return {'form':f,'index':n,'basis':basis,'recovery_scalar':[str(x) for x in lam]}
def ideal_product_form(q1,q2):
 d=disc(q1);vectors=[field_mul(z,w,d) for z in canonical_ideal(q1) for w in canonical_ideal(q2)]
 hnf=lattice_hnf(vectors);out=ideal_form(hnf,d)
 return {'product_generators':vectors,'hnf':hnf,'index':out['index'],'norm_form':out['form'],'reduced':reduce_form(out['form'])['reduced']}
def representations(q,m,primitive=False):
 q=valid(q)
 if type(m) is not int or m<0:raise ValueError('nonnegative represented value required')
 a,b,c=q;D=abs(disc(q));Y=isqrt(4*a*m//D);ans=[]
 for y in range(-Y,Y+1):
  square=4*a*m-D*y*y
  z=isqrt(square)
  if z*z!=square:continue
  for s in sorted(set((z,-z))):
   if (s-b*y)%(2*a)==0:
    x=(s-b*y)//(2*a)
    if not primitive or gcd(x,y)==1:ans.append((x,y))
 return sorted(ans)

def run_checks():
 counts=Counter()
 def ck(k,condition):
  counts[k]+=1
  if not condition:raise AssertionError((k,counts[k]))
 cap={'reduction':reduce_form((24,83,72)),'nonstandard_ideal_basis':ideal_form(((6,1),(10,2)),-23)}
 rng=random.Random(20261009)
 for d in range(-1,-401,-1):
  if d%4 not in (0,1):continue
  forms=enumerate_reduced(d)
  # Explicit representative transformations including exceptional boundaries.
  for f in forms:
   for i in range(16):
    u=I
    for j in range(8):u=mm(u,(1,rng.randrange(-8,9),0,1) if j%2==0 else S)
    q=change(f,u);cert=reduce_form(q)
    ck('transformed_reduction',tuple(cert['reduced'])==f and verify_reduction(q,cert))
   for m in range(1,25):
    # Independent bounded lattice search; square-completion bounds are proved in the page.
    exact=representations(f,m);a,b,c=f;Y=isqrt(4*a*m//abs(d));X=isqrt(4*a*m)+abs(b)*Y+1
    brute=[(x,y) for x in range(-X,X+1) for y in range(-Y,Y+1) if val(f,x,y)==m]
    ck('representation_complete_square',exact==sorted(brute))
  if not fundamental(d):continue
  table={}
  for f,g in product(forms,repeat=2):
   c=compose(f,g);h=tuple(c['reduced']);table[(f,g)]=h
   ck('composition_certificate',verify_composition(f,g,c) and h in forms)
   alt=ideal_product_form(f,g)
   ck('independent_lattice_product',tuple(alt['reduced'])==h and alt['index']==f[0]*g[0])
   # Different input representatives must give the same output.
   v=(1,3,0,1);ck('composition_representative_choice',tuple(compose(change(f,S),change(g,v))['reduced'])==h)
   cf=c['common_forms'];A,B,Cc=c['raw_product'];a1=cf[0][0];a2=cf[1][0]
   for x,y,u,v in [(1,2,-1,3),(0,1,1,0),(2,-3,4,1)]:
    X=x*u-Cc*y*v;Y=a1*x*v+a2*y*u+B*y*v
    ck('bilinear_composition_identity',val(cf[0],x,y)*val(cf[1],u,v)==val((A,B,Cc),X,Y))
  e=next(f for f in forms if f[0]==1)
  for f in forms:
   fi=tuple(reduce_form((f[0],-f[1],f[2]))['reduced'])
   ck('group_identity_inverse',table[(e,f)]==f and table[(f,fi)]==e)
  for f,g,h in product(forms,repeat=3):ck('group_associativity',table[(table[(f,g)],h)]==table[(f,table[(g,h)])])
  if d in (-20,-23,-84):
   cap[str(d)]={'forms':forms,'table':[[forms.index(table[(f,g)]) for g in forms] for f in forms],'indexing':'zero-based indices in forms list'}
 # Very large coordinate changes test exact big-integer certificates.
 for bits in (64,128,256,512,1024,2048):
  for f in [(1,1,1),(1,0,1),(2,1,3),(5,4,5),(2,2,2)]:
   u=mm((1,(1<<bits)+3,0,1),S);u=mm(u,(1,-(1<<(bits//2))-7,0,1))
   q=change(f,u);cert=reduce_form(q);ck('large_integer_reduction',tuple(cert['reduced'])==tuple(reduce_form(f)['reduced']) and verify_reduction(q,cert))
 # Recover ideal bases after integral complex scaling and positive SL2 changes.
 for d in (-3,-4,-7,-8,-20,-23,-84,-163):
  for f in enumerate_reduced(d):
   basis=canonical_ideal(f)
   for z in [(1,0),(-2,1),(3,-2)]:
    old=[field_mul(z,v,d) for v in basis]
    for u in [I,S,(3,5,1,2),(1,-7,0,1)]:
     p,q,r,s=u;new=((p*old[0][0]+r*old[1][0],p*old[0][1]+r*old[1][1]),(q*old[0][0]+s*old[1][0],q*old[0][1]+s*old[1][1]))
     out=ideal_form(new,d);ck('ideal_oriented_inverse',out['form']==change(f,u))
 # Validate genuinely malformed inputs and tampered final matrices/steps.
 bad=[lambda:reduce_form((0,0,0)),lambda:reduce_form((1,0,-2)),lambda:reduce_form((True,0,1)),lambda:compose((1,0,4),(1,0,4)),lambda:compose((1,1,6),(1,0,5)),lambda:compose((2,2,2),(2,2,2)),lambda:ideal_form(((2,0),(0,-1)),-4),lambda:ideal_form(((5,0),(0,1)),-23)]
 for fn in bad:
  try:fn()
  except ValueError:ck('invalid_input_rejected',True)
  else:ck('invalid_input_rejected',False)
 import copy
 for field in ['input','matrix','reduced','steps']:
  cert=copy.deepcopy(cap['reduction'])
  if field=='input':cert[field]=(24,-83,72)
  elif field=='matrix':cert[field]=(2,5,1,-3)
  elif field=='reduced':cert[field]=(2,-1,3)
  else:cert[field][0]['after']=(24,-13,3)
  ck('damaged_reduction_rejected',not verify_reduction((24,83,72),cert))
 cert=compose((2,1,3),(2,1,3));cap['A_squared']=cert;cap['A_times_inverse']=compose((2,1,3),(2,-1,3))
 for field,value in [('inputs',[(2,1,3),(2,-1,3)]),('common_B',1),('raw_product',(6,5,3)),('second_input_matrix',(0,1,1,0)),('reduced',(1,1,6))]:
  damaged=copy.deepcopy(cert);damaged[field]=value
  ck('damaged_composition_rejected',not verify_composition((2,1,3),(2,1,3),damaged))
 cap['representations_13']={str(f):representations(f,13,True) for f in enumerate_reduced(-23)}
 ck('thirteen_principal_obstruction',cap['representations_13'][str((1,1,6))]==[])
 ck('thirteen_nonprincipal_witness',(2,1) in cap['representations_13'][str((2,1,3))])
 cap['A_squared_direct_ideal']=ideal_product_form((2,1,3),(2,1,3))
 cap['prime_ideals_13']=[]
 for r,basis,expected in [(5,((13,0),(-5,1)),(2,1,3)),(9,((13,0),(4,1)),(2,-1,3))]:
  info=ideal_form(basis,-23);cert=reduce_form(info['form'])
  ck('thirteen_prime_ideal',(r*r-r+6)%13==0 and info['index']==13 and tuple(cert['reduced'])==expected and verify_reduction(info['form'],cert))
  cap['prime_ideals_13'].append({'theta_mod_13':r,'positive_basis':basis,'norm_form':info['form'],'reduction':cert})
 return {'status':'PASS','assertions':sum(counts.values()),'groups':dict(counts),'capstone':cap,'limits':['Exact finite tests complement, not replace, the proofs.','Coprime representative search is terminating; no polynomial bit-time claim is made.','Composition excludes nonfundamental discriminants and nonprimitive forms.']}
if __name__=='__main__':
 ap=argparse.ArgumentParser();ap.add_argument('--output',type=Path);args=ap.parse_args();result=run_checks();text=json.dumps(result,ensure_ascii=False,indent=2)+'\n'
 if args.output:args.output.write_text(text)
 print(text,end='')