#!/usr/bin/env python3
"""Exact conic recovery, parametrization, polarity and incidence certificates.

Standard library only. Uses Fraction over Q and modular inverses over the
listed odd prime fields. Polynomial identities are checked coefficientwise.
The accompanying prose proves the general statements; finite tests are separate.
"""
from fractions import Fraction as F
from itertools import product, combinations, permutations
from collections import Counter
from math import gcd, lcm
from pathlib import Path
import argparse, json
CHECKS=Counter()
def check(v,label):
    if not v: raise RuntimeError(label)
    CHECKS[label]+=1

def scalar(x,p=0):
    x=F(x)
    return (x.numerator*pow(x.denominator,-1,p))%p if p else x

def inverse(x,p=0):
    if not x: raise ValueError('division by zero')
    return pow(int(x)%p,-1,p) if p else 1/F(x)

def vector(x,p=0): return [scalar(a,p) for a in x]
def transpose(A): return list(map(list,zip(*A)))
def dot(x,y,p=0): return scalar(sum(a*b for a,b in zip(x,y)),p)
def mv(A,x,p=0): return [dot(row,x,p) for row in A]
def mm(A,B,p=0): return [[dot(row,col,p) for col in transpose(B)] for row in A]
def cross(a,b,p=0): return vector([a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]],p)
def det(A,p=0):
    n=len(A)
    if n==0:return scalar(1,p)
    return scalar(sum((-1)**j*A[0][j]*det([r[:j]+r[j+1:] for r in A[1:]],p) for j in range(n)),p)
def rref(A,p=0):
    B=[vector(r,p) for r in A];pivots=[];i=0
    for j in range(len(B[0])):
        pivot=next((k for k in range(i,len(B)) if B[k][j]),None)
        if pivot is None:continue
        B[i],B[pivot]=B[pivot],B[i];a=inverse(B[i][j],p);B[i]=vector([a*x for x in B[i]],p)
        for k in range(len(B)):
            if k!=i:
                a=B[k][j];B[k]=vector([x-a*y for x,y in zip(B[k],B[i])],p)
        pivots.append(j);i+=1
        if i==len(B):break
    return B,pivots

def nullspace(A,p=0):
    B,piv=rref(A,p);n=len(B[0]);out=[]
    for j in range(n):
        if j in piv:continue
        x=[scalar(0,p)]*n;x[j]=scalar(1,p)
        for i,k in enumerate(piv):x[k]=scalar(-B[i][j],p)
        out.append(x)
    return out

def invmat(A,p=0):
    n=len(A);B,piv=rref([vector(r,p)+vector([i==j for j in range(n)],p) for i,r in enumerate(A)],p)
    if piv[:n]!=list(range(n)):raise ValueError('singular matrix')
    return [r[n:] for r in B]

def canon(x,p=0):
    x=vector(x,p)
    if not any(x):raise ValueError('zero vector has no projective class')
    if p:
        k=next(a for a in x if a);return tuple(scalar(a*inverse(k,p),p) for a in x)
    den=lcm(*(a.denominator for a in x));nums=[int(a*den) for a in x];g=gcd(*nums);nums=[a//g for a in nums]
    if next(a for a in nums if a)<0:nums=[-a for a in nums]
    return tuple(nums)

def quadratic(H,x,p=0):return dot(x,mv(H,x,p),p)
def require_conic(H,p=0):
    if p==2 or len(H)!=3 or any(len(r)!=3 for r in H):raise ValueError('three coordinates and characteristic not two required')
    if H!=transpose(H) or not det(H,p):raise ValueError('invertible symmetric matrix required')

def recover(points,p=0):
    if len(points)!=5:raise ValueError('five points required')
    points=[canon(x,p) for x in points]
    if len(set(points))!=5 or any(not det(transpose(list(t)),p) for t in combinations(points,3)):
        raise ValueError('five distinct points with no three collinear required')
    rows=[vector([x*x,y*y,z*z,x*y,x*z,y*z],p) for x,y,z in points];basis=nullspace(rows,p)
    if len(basis)!=1:raise RuntimeError('five-point rank certificate failed')
    a,b,c,d,e,f=basis[0];half=inverse(scalar(2,p),p)
    H=[vector([a,d*half,e*half],p),vector([d*half,b,f*half],p),vector([e*half,f*half,c],p)]
    require_conic(H,p);return H,rows,basis[0]

def parameter_matrix(H,point,p=0):
    require_conic(H,p);point=vector(canon(point,p),p)
    if quadratic(H,point,p):raise ValueError('given point is not on the conic')
    h=mv(H,point,p);i=next(i for i,a in enumerate(h) if a);v=[scalar(0,p)]*3;v[i]=inverse(h[i],p)
    half=inverse(scalar(2,p),p);qv=quadratic(H,v,p);v=vector([a-qv*half*b for a,b in zip(v,point)],p)
    r=nullspace([mv(H,point,p),mv(H,v,p)],p)[0];d=quadratic(H,r,p)
    T=transpose([point,r,vector([-d*half*a for a in v],p)])
    expected=[[0,0,-d*half],[0,d,0],[-d*half,0,0]]
    check(mm(transpose(T),mm(H,T,p),p)==[vector(row,p) for row in expected],'constructed_congruence_coefficient_identity')
    check(bool(det(T,p)) and bool(d),'constructed_basis_invertible')
    return T,d

def nu(pair,p=0):
    s,t=pair;return vector([s*s,s*t,t*t],p)
def reverse_parameter(T,x,p=0):
    X,Y,Z=mv(invmat(T,p),x,p)
    if scalar(X*Z-Y*Y,p):raise ValueError('point not in parametrized conic')
    return canon((X,Y) if X else (Y,Z),p)
def pascal(points,p=0):
    if len(points)!=6 or len({canon(x,p) for x in points})!=6:raise ValueError('six distinct points required')
    edges=[canon(cross(points[i],points[(i+1)%6],p),p) for i in range(6)]
    meets=[canon(cross(edges[i],edges[i+3],p),p) for i in range(3)]
    L=canon(cross(meets[0],meets[1],p),p)
    check(all(dot(L,x,p)==0 for x in meets),'Pascal_all_three_incidences')
    check(len(set(meets))==3,'Pascal_three_distinct_meets')
    return edges,meets,L

def brianchon(H,points,p=0):
    tangents=[canon(mv(H,x,p),p) for x in points];vertices=[canon(cross(tangents[i],tangents[(i+1)%6],p),p) for i in range(6)]
    diagonals=[canon(cross(vertices[i],vertices[i+3],p),p) for i in range(3)];pole=canon(cross(diagonals[0],diagonals[1],p),p)
    check(all(dot(l,pole,p)==0 for l in diagonals),'Brianchon_all_three_incidences')
    check(len(set(diagonals))==3,'Brianchon_three_distinct_diagonals')
    return tangents,vertices,diagonals,pole

def crossratio(A,B,C,D,p=0):
    if len({canon(x,p) for x in [A,B,C,D]})!=4:raise ValueError('four distinct ordered points required')
    bracket=lambda x,y:scalar(x[0]*y[1]-x[1]*y[0],p)
    return scalar(bracket(A,C)*bracket(B,D)*inverse(bracket(A,D)*bracket(B,C),p),p)

class Poly:
    """Sparse integer polynomials for coefficientwise identities, not evaluations."""
    def __init__(self,x=0):self.c=dict(x) if isinstance(x,dict) else ({(0,)*5:x} if x else {})
    def __add__(self,b):
        if not isinstance(b,Poly):b=Poly(b)
        c=self.c.copy()
        for k,v in b.c.items():c[k]=c.get(k,0)+v
        return Poly({k:v for k,v in c.items() if v})
    __radd__=__add__
    def __neg__(self):return Poly({k:-v for k,v in self.c.items()})
    def __sub__(self,b):return self+-aspoly(b)
    def __rsub__(self,b):return aspoly(b)+-self
    def __mul__(self,b):
        b=aspoly(b);c=Counter()
        for k,v in self.c.items():
            for l,w in b.c.items():c[tuple(x+y for x,y in zip(k,l))]+=v*w
        return Poly({k:v for k,v in c.items() if v})
    __rmul__=__mul__
    def __pow__(self,n):
        out=Poly(1)
        for _ in range(n):out=out*self
        return out

def aspoly(x):return x if isinstance(x,Poly) else Poly(x)
def pdet(A):
    if len(A)==1:return aspoly(A[0][0])
    return sum(((-1)**j)*aspoly(A[0][j])*pdet([r[:j]+r[j+1:] for r in A[1:]]) for j in range(len(A)))
def polynomial_identities():
    a,b,c,u,v=[Poly({tuple(int(i==j) for i in range(5)):1}) for j in range(5)]
    P=[1,0,-a*b];Q=[b+c-1,b*c,b*c];R=[1,c,c*(1+a)-a];L=[a*b*c,-a*(b+c-1)-c,c]
    for X in [P,Q,R]:check(not sum(aspoly(x)*y for x,y in zip(X,L)).c,'universal_Pascal_dot_polynomial')
    # Rename c,u,v as the generic sixth coordinates U,V,W.
    U,V,W=c,u,v
    Q=[b*(b-1)*U+V-W,b*(b*V-W),b*(b*V-W)];R=[U,V,(a+1)*V-a*U]
    check(not (pdet(transpose([P,Q,R]))+a*b*(b-1)*(U*W-V*V)).c,'universal_Pascal_converse_determinant_polynomial')
    A=[[a*a,2*a*b,b*b],[a*c,a*u+b*c,b*u],[c*c,2*c*u,u*u]]
    check(not (pdet(A)-(a*u-b*c)**3).c,'symmetric_square_determinant_polynomial')

def models():
    base=[[0,0,1],[0,-2,0],[1,0,0]];fields=[0,3,5,7,11,13];counts={}
    for p in fields:
        model_count=0
        for k in range(1,9):
            M=[[1,k,0],[0,1,k+1],[1,0,1-k*(k+1)]] # determinant one
            check(det(M,p)==1,'nonorthogonal_model_matrix_determinant')
            H=mm(transpose(invmat(M,p)),mm(base,invmat(M,p),p),p);require_conic(H,p)
            params=[(1,t) for t in (range(p) if p else range(-3,5))]+[(0,1)]
            points=[canon(mv(M,nu(t,p),p),p) for t in params]
            check(len(set(points))==len(params),'complete_parameter_points_distinct')
            T,d=parameter_matrix(H,points[(k+1)%len(points)],p)
            for pair,point in zip(params,points):
                check(quadratic(H,point,p)==0,'transformed_parameter_on_conic')
                back=reverse_parameter(M,point,p);check(back==canon(pair,p),'original_parameter_inverse')
                alt=reverse_parameter(T,point,p);check(canon(mv(T,nu(alt,p),p),p)==point,'constructed_parameter_inverse_cover')
                ell=mv(H,point,p);check(dot(ell,point,p)==0 and quadratic(invmat(H,p),ell,p)==0,'tangent_dual_quadric_and_contact')
                check(canon(mv(invmat(H,p),ell,p),p)==point,'pole_polar_inverse')
            for a,b in zip(points,points[1:]):
                check(bool(canon(cross(a,b,p),p)),'distinct_chord_nonzero')
                # A+lambda B and A-lambda B form a harmonic pair.
                lam=scalar(2,p);P=vector([x+lam*y for x,y in zip(a,b)],p);Q=vector([x-lam*y for x,y in zip(a,b)],p)
                check(dot(mv(H,P,p),Q,p)==0,'secant_harmonic_polar_incidence')
                check(crossratio((1,0),(0,1),(1,lam),(1,-lam),p)==scalar(-1,p),'harmonic_crossratio_minus_one')
            if p:
                plane=[(1,x,y) for x,y in product(range(p),repeat=2)]+[(0,1,y) for y in range(p)]+[(0,0,1)]
                zero_set={canon(x,p) for x in plane if quadratic(H,x,p)==0}
                check(zero_set==set(points),'finite_field_entire_plane_zero_locus_equals_parameters')
            if len(points)>=5:
                for ix in [list(range(5)),list(range(len(points)-5,len(points)))]:
                    K,rows,coeff=recover([points[i] for i in ix],p)
                    check(canon(sum(K,[]),p)==canon(sum(H,[]),p),'five_points_recover_original_conic')
                    check(len(rref(rows,p)[1])==5,'five_point_constraint_rank')
            if len(points)>=6:
                chosen=points[:6];_,_,L=pascal(chosen,p);_,_,_,W=brianchon(H,chosen,p)
                check(canon(mv(invmat(H,p),L,p),p)==W,'Brianchon_point_is_Pascal_pole')
            A2=[[2,1],[1,1]];four=[(1,0),(0,1),(1,1),(2,1)]
            if len({canon(x,p) for x in four})==4:
                check(crossratio(*four,p)==crossratio(*[mv(A2,x,p) for x in four],p),'four_point_crossratio_invariant')
            model_count+=1
        counts[str(p or 'Q')]=model_count
    return counts

def all_orders():
    H=[[0,0,1],[0,-2,0],[1,0,0]];points=[nu(t,5) for t in [(0,1),(1,0),(1,1),(1,2),(1,3),(1,4)]]
    for order in permutations(range(6)):
        row=[points[i] for i in order];_,_,L=pascal(row,5);_,_,_,W=brianchon(H,row,5)
        check(canon(mv(invmat(H,5),L,5),5)==W,'all_720_orders_polar_duality_over_F5')
    return 720

def capstone():
    M=[[1,1,0],[0,1,1],[1,0,1]];H=[[-2,0,1],[0,-2,1],[1,1,0]];parameters=[(0,1),(1,0),(1,1),(1,2),(1,3),(1,-1)]
    points=[canon(mv(M,nu(t))) for t in parameters];K,rows,coeff=recover(points[:5])
    check(canon(sum(K,[]))==canon(sum(H,[])),'capstone_recovered_H')
    check(det([row[1:] for row in rows])==6,'capstone_rank_five_minor')
    triples=[det(transpose(list(x))) for x in combinations(points[:5],3)]
    check(triples==[1,4,3,1,1,1,2,3,6,1],'capstone_ten_noncollinearity_determinants')
    edges,meets,L=pascal(points);tangents,vertices,diags,W=brianchon(H,points)
    check(L==(3,-2,3) and W==(1,11,14),'capstone_Pascal_line_Brianchon_point')
    polars=[canon(mv(H,p)) for p in [(1,-1,0),(2,-1,1),(1,1,2)]]
    check(polars==[(1,-1,0),(3,-3,-1),(0,0,1)],'capstone_three_polar_lines')
    # Quadratic-field coordinates are stored as a+b*sqrt(2), with rational a,b.
    mul=lambda x,y:(x[0]*y[0]+2*x[1]*y[1],x[0]*y[1]+x[1]*y[0])
    for sign in [-1,1]:
        r=[(2,sign),(1,sign),(3,0)];ell=[(-1,-2*sign),(1,-2*sign),(3,2*sign)]
        for i in range(3):check(tuple(sum(H[i][j]*r[j][k] for j in range(3)) for k in range(2))==ell[i],'quadratic_field_tangent_coefficients')
        q=[0,0]
        for i,j in product(range(3),repeat=2):
            z=mul(r[i],r[j]);q=[q[k]+H[i][j]*z[k] for k in range(2)]
        check(q==[0,0],'quadratic_field_contact_on_conic')
        check(tuple(sum(ell[i][k]*[2,-1,1][i] for i in range(3)) for k in range(2))==(0,0),'quadratic_field_tangent_through_rational_point')
    N=[[1,1,0],[0,1,1],[6,0,1]];at_infinity=canon(mv(N,[1,0,-6]));NL=canon(mv(transpose(invmat(N)),[6,1,1]))
    check(at_infinity==(1,-6,0) and NL==(6,1,6) and dot(at_infinity,NL)==0,'projective_infinity_migration')
    source=[(0,1),(1,1),(1,0),(2,1)];good=[(1,1),(2,1),(3,1),(7,3)];bad=good[:3]+[(5,2)]
    check(crossratio(*source)==crossratio(*good)==F(1,2),'four_point_compatible_target')
    check(crossratio(*bad)==F(2,3),'four_point_incompatible_target')
    # Local obstruction is checked exhaustively; primitive-integer descent is in prose.
    check(all(x==y==0 for x,y in product(range(3),repeat=2) if (x*x+y*y)%3==0),'mod_three_obstruction_first_stage')
    return {'points':points,'five_point_rows':rows,'kernel':canon(coeff),'triple_determinants':triples,'H':H,'H_inverse':invmat(H),'edges':edges,'Pascal_points':meets,'Pascal_line':L,'tangents':tangents,'vertices':vertices,'diagonals':diags,'Brianchon_point':W,'polars':polars,'infinite_intersection':at_infinity,'crossratios':{'source':F(1,2),'bad_target':F(2,3)}}

def rejected_inputs():
    H=[[0,0,1],[0,-2,0],[1,0,0]]
    bad=[lambda:canon([0,0,0]),lambda:invmat([[1,0],[0,0]]),lambda:require_conic(H,2),lambda:parameter_matrix(H,[1,0,1]),lambda:recover([(1,t,0) for t in range(4)]+[(0,0,1)]),lambda:pascal([(1,0,0)]*6),lambda:crossratio((1,0),(1,0),(0,1),(1,1))]
    for f in bad:
        try:f()
        except ValueError:check(True,'invalid_input_rejected')
        else:raise RuntimeError('invalid input accepted')

def serial(x):
    if isinstance(x,F):return str(x)
    if isinstance(x,dict):return {str(k):serial(v) for k,v in x.items()}
    if isinstance(x,(list,tuple)):return [serial(v) for v in x]
    return x

def main():
    CHECKS.clear()
    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path,required=True);args=ap.parse_args()
    polynomial_identities();counts=models();orders=all_orders();result=capstone();rejected_inputs()
    out={'status':'PASS','checks':sum(CHECKS.values()),'groups':dict(CHECKS),'field_models':counts,'ordered_hexagons_F5':orders,'scope':'Exact field operations; coefficientwise polynomial identities and finite model tests. General proofs and the rational nonexistence descent are in the prose. Arithmetic-operation counts are not bit-cost bounds.','capstone':result}
    args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(json.dumps(serial(out),ensure_ascii=False,indent=2)+'\n');print(json.dumps({'status':'PASS','checks':sum(CHECKS.values())}))
if __name__=='__main__':main()