#!/usr/bin/env python3
"""Exact, standard-library lattice counting and half-open cone certificates.

The proof pages establish the all-parameter theorems. This program independently
compares cone fundamental-domain counts, direct inequalities, polygon refinement,
and the concrete capstone formulas. No floating-point geometry is used.
"""
from fractions import Fraction as F
from pathlib import Path
from itertools import product,combinations
from collections import Counter
from math import gcd,lcm,comb,factorial
import argparse,json,random

def rational(x):
 if isinstance(x,bool) or not isinstance(x,(int,F,str)):raise ValueError('exact rational required')
 return F(x)
def vertices(v):
 if not v:raise ValueError('nonempty vertex list required')
 d=len(v[0])
 if d==0 or any(len(p)!=d for p in v):raise ValueError('consistent positive ambient dimension required')
 out=tuple(tuple(rational(x) for x in p) for p in v)
 if len(set(out))!=len(out):raise ValueError('repeated vertex')
 return out
def floor(x):return x.numerator//x.denominator
def ceil(x):return -floor(-x)
def determinant(a):
 a=[list(map(F,row)) for row in a];n=len(a)
 if any(len(row)!=n for row in a):raise ValueError('square matrix required')
 out=F(1)
 for j in range(n):
  p=next((i for i in range(j,n) if a[i][j]),None)
  if p is None:return F(0)
  if p!=j:a[p],a[j]=a[j],a[p];out=-out
  pivot=a[j][j];out*=pivot
  for i in range(j+1,n):
   c=a[i][j]/pivot
   for k in range(j+1,n):a[i][k]-=c*a[j][k]
 return out
def inverse(a):
 n=len(a);a=[list(map(F,row))+[F(i==j) for j in range(n)] for i,row in enumerate(a)]
 for j in range(n):
  p=next((i for i in range(j,n) if a[i][j]),None)
  if p is None:raise ValueError('dependent generators')
  a[p],a[j]=a[j],a[p];c=a[j][j];a[j]=[x/c for x in a[j]]
  for i in range(n):
   if i!=j:
    c=a[i][j];a[i]=[x-c*y for x,y in zip(a[i],a[j])]
 return tuple(tuple(row[n:]) for row in a)
def coordinate_map(cols):
 k=len(cols);d=len(cols[0])
 if any(len(v)!=d for v in cols):raise ValueError('inconsistent generators')
 rows=next((r for r in combinations(range(d),k) if determinant([[cols[j][i] for j in range(k)] for i in r])),None)
 if rows is None:raise ValueError('dependent generators')
 inv=inverse([[cols[j][i] for j in range(k)] for i in rows])
 def coords(p):
  if len(p)!=d:return None
  c=tuple(sum(inv[j][i]*p[rows[i]] for i in range(k)) for j in range(k))
  if any(sum(c[j]*cols[j][i] for j in range(k))!=p[i] for i in range(d)):return None
  return c
 return coords

def fundamental_domain(cols,strict=()):
 cols=tuple(tuple(v) for v in cols);k=len(cols);d=len(cols[0]);strict=frozenset(strict)
 if any(type(x) is not int for v in cols for x in v):raise ValueError('integer generators required')
 if any(type(j) is not int or not 0<=j<k for j in strict):raise ValueError('invalid strict coordinate')
 c=coordinate_map(cols)
 limits=[range(sum(min(0,v[i]) for v in cols),sum(max(0,v[i]) for v in cols)+1) for i in range(d)]
 out=[]
 for p in product(*limits):
  a=c(p)
  if a is not None and all((0<x<=1 if j in strict else 0<=x<1) for j,x in enumerate(a)):out.append(p)
 return out

def simplex_series(v,q=None,strict=()):
 v=vertices(v);r=len(v)-1
 if q is None:q=lcm(*(x.denominator for p in v for x in p))
 if type(q) is not int or q<=0 or any((x*q).denominator!=1 for p in v for x in p):raise ValueError('invalid common denominator')
 cols=tuple(tuple(int(x*q) for x in p)+(q,) for p in v)
 pts=fundamental_domain(cols,strict);hist=dict(sorted(Counter(p[-1] for p in pts).items()))
 return {'vertices':v,'dimension':r,'q':q,'columns':cols,'strict':sorted(strict),'points':pts,'numerator':hist}

def binomial_polynomial(n,r):
 a=1
 for j in range(r):a*=n-j
 return a//factorial(r)
def count_from_series(cert,t):
 if type(t) is not int or t<0:raise ValueError('nonnegative integer dilation required')
 q=cert['q'];r=cert['dimension']
 return sum(c*comb((t-int(h))//q+r,r) for h,c in cert['numerator'].items() if t>=int(h) and (t-int(h))%q==0)
def extend_quasipolynomial(cert,t):
 if type(t) is not int:raise ValueError('integer evaluation required')
 q=cert['q'];r=cert['dimension']
 if any(not 0<=int(h)<q*(r+1) for h in cert['numerator']):raise ValueError('closed-series degree bound required')
 return sum(c*binomial_polynomial((t-int(h))//q+r,r) for h,c in cert['numerator'].items() if (t-int(h))%q==0)

def cross(a,b,c):return (b[0]-a[0])*(c[1]-a[1])-(b[1]-a[1])*(c[0]-a[0])
def hull(points):
 p=sorted(set(tuple(map(F,x)) for x in points))
 if any(len(v)!=2 for v in p):raise ValueError('planar points required')
 lo=[];hi=[]
 for v in p:
  while len(lo)>=2 and cross(lo[-2],lo[-1],v)<=0:lo.pop()
  lo.append(v)
 for v in reversed(p):
  while len(hi)>=2 and cross(hi[-2],hi[-1],v)<=0:hi.pop()
  hi.append(v)
 h=tuple(lo[:-1]+hi[:-1])
 if len(h)<3:raise ValueError('positive area required')
 return h

def polygon_points(v,t=1,interior=False):
 v=hull(v)
 if type(t) is not int or t<0 or (interior and t==0):raise ValueError('interior requires positive integer t')
 if t==0:return [(0,0)]
 p=[tuple(x*t for x in a) for a in v];out=[]
 for x in range(ceil(min(a[0] for a in p)),floor(max(a[0] for a in p))+1):
  for y in range(ceil(min(a[1] for a in p)),floor(max(a[1] for a in p))+1):
   c=[cross(a,b,(x,y)) for a,b in zip(p,p[1:]+p[:1])]
   if all(z>0 if interior else z>=0 for z in c):out.append((x,y))
 return out

def polygon_certificate(v):
 v=hull(v);q=lcm(*(x.denominator for p in v for x in p));tri=[(v[0],v[i],v[i+1]) for i in range(1,len(v)-1)]
 columns=[tuple(tuple(int(x*q) for x in p)+(q,) for p in t) for t in tri]
 maps=[coordinate_map(c) for c in columns];seed=1
 while True:
  weights=[(seed+1)**i for i in range(len(v))];W=sum(weights)
  w=tuple(sum(weights[i]*v[i][j] for i in range(len(v)))/W for j in range(2))+(F(1),)
  beta=[f(w) for f in maps]
  if all(all(x!=0 for x in b) for b in beta):break
  seed+=1
 closed=Counter();opened=Counter();pieces=[]
 for t,b in zip(tri,beta):
  J=[i for i,x in enumerate(b) if x<0];K=[i for i in range(3) if i not in J]
  a=simplex_series(t,q,J);o=simplex_series(t,q,K);closed.update(a['numerator']);opened.update(o['numerator'])
  pieces.append({'vertices':t,'beta':b,'closed_strict':J,'closed_numerator':a['numerator'],'interior_numerator':o['numerator']})
 return {'vertices':v,'dimension':2,'q':q,'direction':w,'pieces':pieces,'numerator':dict(sorted(closed.items())),'interior_numerator':dict(sorted(opened.items()))}

def verify_polygon_certificate(v,cert):
 try:
  fresh=polygon_certificate(v)
  normalize=lambda obj:json.loads(json.dumps(obj,default=serial))
  return normalize(fresh)==normalize(cert)
 except (ValueError,KeyError,TypeError):return False

def simplex_points(v,t,interior=False):
 v=vertices(v)
 if type(t) is not int or t<0 or (interior and t==0):raise ValueError('invalid dilation')
 if t==0:return [(0,)*len(v[0])]
 cols=tuple(tuple(p)+(F(1),) for p in v);coords=coordinate_map(cols);out=[]
 ranges=[range(ceil(t*min(p[i] for p in v)),floor(t*max(p[i] for p in v))+1) for i in range(len(v[0]))]
 for p in product(*ranges):
  a=coords(p+(t,))
  if a is not None and all(x>0 if interior else x>=0 for x in a):out.append(p)
 return out

def primitive_triangulation(v):
 v=hull(v)
 if any(x.denominator!=1 for p in v for x in p):raise ValueError('integer polygon required')
 v=tuple(tuple(int(x) for x in p) for p in v);tri=[(v[0],v[i],v[i+1]) for i in range(1,len(v)-1)]
 used=set(v);pts=polygon_points(v)
 for p in pts:
  if p in used:continue
  nxt=[];found=False
  for a,b,c in tri:
   signs=[cross(a,b,p),cross(b,c,p),cross(c,a,p)]
   if min(signs)<0:nxt.append((a,b,c));continue
   found=True
   if min(signs)>0:nxt.extend([(a,b,p),(b,c,p),(c,a,p)])
   else:
    j=signs.index(0);u,w,z=[(a,b,c),(b,c,a),(c,a,b)][j];nxt.extend([(u,p,z),(p,w,z)])
  if not found:raise AssertionError('point missed by triangulation')
  tri=nxt;used.add(p)
 edges=Counter(tuple(sorted((a,b))) for t in tri for a,b in zip(t,t[1:]+t[:1]))
 B=sum(gcd(abs(b[0]-a[0]),abs(b[1]-a[1])) for a,b in zip(v,v[1:]+v[:1]));I=len(polygon_points(v,interior=True));A2=sum(cross((0,0),a,b) for a,b in zip(v,v[1:]+v[:1]))
 return {'vertices':v,'area_twice':A2,'I':I,'B':B,'triangles':tri,'triangle_determinants':[cross(*t) for t in tri],'V':len(used),'E':len(edges),'T':len(tri),'boundary_edges':sum(c==1 for c in edges.values()),'edge_multiplicities':dict(Counter(edges.values()))}

def rational_triangle_count(t):
 if t<0:return 0 # actual inequality count, deliberately not quasipolynomial extension
 return sum((t-3*y)//2+1 for y in range(t//3+1))
def rational_triangle_extension(t):
 k,s=divmod(t,6);return 3*k*k+(s+3)*k+[1,1,2,3,4,5][s]
def reeve_formula(t,m):return F(m,6)*t**3+t*t+F(12-m,6)*t+1

def run_checks():
 counts=Counter()
 def ck(label,ok):
  if not ok:raise AssertionError(label)
  counts[label]+=1
 rng=random.Random(20261009)
 triangle=[(0,0),(4,0),(0,3)];p=primitive_triangulation(triangle);rat=[(0,0),(F(1,2),0),(0,F(1,3))]
 cert=polygon_certificate(rat)
 cap={'pick':p,'rational_triangle':cert,'rational_closed_0_to_18':[rational_triangle_count(t) for t in range(19)],'rational_interior_1_to_18':[rational_triangle_count(t-6) for t in range(1,19)],'negative_seven':rational_triangle_extension(-7),'geometric_negative_seven':rational_triangle_count(7),'period_collapse':polygon_certificate([(0,0),(2,0),(0,F(1,2))])}
 polys={hull(triangle)}
 for _ in range(110):
  pts=[(rng.randrange(-3,4),rng.randrange(-3,4)) for _ in range(7)]
  try:polys.add(hull(pts))
  except ValueError:pass
 for v in sorted(polys):
  a=primitive_triangulation(v)
  ck('pick_and_refinement',a['area_twice']==2*a['I']+a['B']-2==a['T'] and all(d==1 for d in a['triangle_determinants']) and a['V']-a['E']+a['T']==1 and a['boundary_edges']==a['B'] and set(a['edge_multiplicities'])<={1,2})
  for t in range(1,5):
   ck('pick_dilation_closed',len(polygon_points(v,t))==F(a['area_twice'],2)*t*t+F(a['B'],2)*t+1)
   ck('pick_dilation_interior',len(polygon_points(v,t,True))==F(a['area_twice'],2)*t*t-F(a['B'],2)*t+1)
 # Non-simplex fans, rational translations and non-unimodular cones.
 rational_polys=[rat,[(0,0),(2,0),(0,F(1,2))],[(0,0),(1,0),(1,1),(0,1)],[(F(1,2),0),(2,0),(F(5,2),1),(1,2),(0,1)]]
 for v in list(sorted(polys))[:24]:
  den=rng.choice([1,2,3]);shift=(F(rng.randrange(2),2),F(rng.randrange(2),2))
  rational_polys.append([tuple(x/den+shift[i] for i,x in enumerate(p)) for p in v])
 for v in rational_polys:
  a=polygon_certificate(v);q=a['q'];c=a['numerator'];o=a['interior_numerator']
  ck('half_open_height_and_reflection',max(c)<3*q and min(o)>0 and o==dict(sorted({3*q-h:n for h,n in c.items()}.items())))
  for t in range(0,13):
   actual=len(polygon_points(v,t));ck('cone_closed_against_inequalities',actual==count_from_series(a,t)==extend_quasipolynomial(a,t))
   if t:
    opened=len(polygon_points(v,t,True));b=dict(a,numerator=o)
    ck('cone_interior_against_inequalities',opened==count_from_series(b,t)==extend_quasipolynomial(a,-t))
  # Independently count closed fan pieces and subtract the shared diagonals.
  vv=list(a['vertices'])
  for t in [1,2,5]:
   s=sum(len(simplex_points([vv[0],vv[i],vv[i+1]],t)) for i in range(1,len(vv)-1))
   s-=sum(len(simplex_points([vv[0],vv[i]],t)) for i in range(2,len(vv)-1))
   ck('fan_closed_inclusion_exclusion',s==len(polygon_points(v,t)))
 for t in range(0,241):
  ck('six_residue_formula',rational_triangle_count(t)==rational_triangle_extension(t)==count_from_series(cert,t))
  if t:ck('six_residue_reciprocity',rational_triangle_extension(-t)==rational_triangle_count(t-6))
 # Relative affine slices: vertices need not span the ambient lattice space.
 slices=[[(F(1,2),0),(F(1,2),1)],[(0,0,0),(2,0,2)],[(F(1,3),F(2,3))],[(F(2,3),)],[(0,0),(F(1,2),F(1,3))]]
 for v in slices:
  a=simplex_series(v);r=a['dimension'];o=simplex_series(v,a['q'],range(r+1))
  for t in range(0,25):
   ck('lower_dimensional_closed',count_from_series(a,t)==len(simplex_points(v,t)))
   if t:ck('relative_interior_reciprocity',(-1)**r*extend_quasipolynomial(a,-t)==len(simplex_points(v,t,True))==count_from_series(o,t))
 # The Reeve family has fixed t=1 counts and variable volume/linear coefficient.
 for m in range(1,17):
  v=[(0,0,0),(1,0,0),(0,1,0),(1,1,m)];a=simplex_series(v)
  expected={0:1};
  if m>1:expected[2]=m-1
  ck('reeve_fundamental_domain',a['numerator']==expected)
  for t in range(0,6):
   ck('reeve_closed',len(simplex_points(v,t))==reeve_formula(t,m)==count_from_series(a,t))
   if t:ck('reeve_interior',len(simplex_points(v,t,True))==-reeve_formula(-t,m))
 for t in range(0,101):
  ck('period_collapse',len(polygon_points([(0,0),(2,0),(0,F(1,2))],t))==(t+1)*(t+2)//2)
 # Directly exercise every open/closed facet assignment for non-unimodular cones.
 for cols in [((0,2),(3,2)),((1,1),(2,3)),((0,0,2),(2,0,2),(0,3,2)),((1,0,2),(1,2,2))]:
  k=len(cols);U=tuple(sum(v[i] for v in cols) for i in range(len(cols[0])))
  for bits in product([0,1],repeat=k):
   J=[j for j,b in enumerate(bits) if b];K=[j for j,b in enumerate(bits) if not b]
   pts=fundamental_domain(cols,J);other=fundamental_domain(cols,K)
   ck('all_half_open_complements',{tuple(U[i]-p[i] for i in range(len(U))) for p in pts}==set(other))
 bad=[lambda:vertices([]),lambda:vertices([(0.1,0),(1,0),(0,1)]),lambda:simplex_series([(0,0),(1,0),(2,0)]),lambda:simplex_series(rat,q=2),lambda:primitive_triangulation(rat),lambda:polygon_points(triangle,0,True),lambda:count_from_series(cert,-1),lambda:count_from_series(cert,True),lambda:fundamental_domain(((1,0),(0,1)),[2])]
 for f in bad:
  try:f()
  except ValueError:ck('invalid_input_rejected',True)
  else:ck('invalid_input_rejected',False)
 cap['complement_demo']={'columns':[(0,2),(3,2)],'closed_points':fundamental_domain(((0,2),(3,2)),()),'open_points':fundamental_domain(((0,2),(3,2)),(0,1))}
 import copy
 ck('certificate_round_trip',verify_polygon_certificate(rat,json.loads(json.dumps(cert,default=serial))))
 for key in ['numerator','interior_numerator','direction','pieces']:
  damaged=copy.deepcopy(cert)
  if key=='numerator':damaged[key][0]+=1
  elif key=='interior_numerator':damaged[key][6]+=1
  elif key=='direction':damaged[key]=(F(0),F(0),F(1))
  else:damaged[key][0]['closed_strict']=[0]
  ck('damaged_certificate_rejected',not verify_polygon_certificate(rat,damaged))
 for t in range(1,11):
  transformed={(x+2*y+t,y-t) for x,y in polygon_points(triangle,t)}
  ck('affine_unimodular_point_bijection',transformed==set(polygon_points([(1,-1),(5,-1),(7,2)],t)))
 cap['large_parameter_2026']={'closed':rational_triangle_extension(2026),'interior':rational_triangle_extension(-2026),'direct_closed':rational_triangle_count(2026),'direct_interior':rational_triangle_count(2020)}
 ck('large_parameter_2026',cap['large_parameter_2026']=={'closed':343070,'interior':341044,'direct_closed':343070,'direct_interior':341044})
 return {'status':'PASS','assertions':sum(counts.values()),'groups':dict(counts),'capstone':cap,'limits':['Finite exact tests do not replace proofs.','Only positive dilations have interior-count semantics.','Negative evaluation follows the quasipolynomial, not geometric reflection.','Enumeration cost may be large in the bit size, denominator and dimension.']}

def serial(obj):
 if isinstance(obj,F):return int(obj) if obj.denominator==1 else str(obj)
 raise TypeError(type(obj).__name__)
if __name__=='__main__':
 p=argparse.ArgumentParser();p.add_argument('--output',type=Path);args=p.parse_args();result=run_checks();s=json.dumps(result,default=serial,ensure_ascii=False,indent=2)+'\n'
 if args.output:args.output.write_text(s)
 print(s,end='')