#!/usr/bin/env python3
"""Exact model certificates for planar cycles; no numerical ODE integration.
Python standard library only. Use --output to choose the JSON destination.
Analytic topology/compactness hypotheses are proved in the accompanying pages.
"""
from fractions import Fraction as Q
from collections import Counter
from itertools import product
from pathlib import Path
import argparse,json,copy
checks=Counter()
def check(ok,label):
 checks[label]+=1
 if not ok:raise ArithmeticError(label)
# Sparse bivariate polynomials in x,y, with exact rational coefficients.
def poly(items):return {p:Q(v) for p,v in items.items() if v}
def add(*args):
 d={}
 for a in args:
  for p,v in a.items():d[p]=d.get(p,Q(0))+v
 return poly(d)
def scale(c,a):return poly({p:c*v for p,v in a.items()})
def mul(a,b):
 d={}
 for p,v in a.items():
  for q,w in b.items():
   z=(p[0]+q[0],p[1]+q[1]);d[z]=d.get(z,Q(0))+v*w
 return poly(d)
def derivative(a,j):
 d={}
 for p,v in a.items():
  if p[j]:z=list(p);z[j]-=1;d[tuple(z)]=p[j]*v
 return poly(d)
def value(a,x,y):return sum(v*x**i*y**j for (i,j),v in a.items())
ONE={(0,0):Q(1)};X={(1,0):Q(1)};Y={(0,1):Q(1)};S=add(mul(X,X),mul(Y,Y))

def vector_field(kappa,delta,omega=Q(1)):
 radial=add(ONE,scale(-1,S))
 P=add(scale(-omega,Y),mul(X,radial),scale(kappa,X),scale(delta,ONE))
 R=add(scale(omega,X),mul(Y,radial))
 return P,R

def identities(kappa,delta,omega):
 P,R=vector_field(kappa,delta,omega)
 radial=add(mul(X,P),mul(Y,R))
 angular=add(mul(X,R),scale(-1,mul(Y,P)))
 div=add(derivative(P,0),derivative(R,1))
 check(radial==add(mul(S,add(ONE,scale(-1,S))),scale(kappa,mul(X,X)),scale(delta,X)),'symbolic_radial_numerator')
 check(angular==add(scale(omega,S),scale(-kappa,mul(X,Y)),scale(-delta,Y)),'symbolic_angular_numerator')
 check(div==add(scale(2+kappa,ONE),scale(-4,S)),'symbolic_divergence')
 return P,R

def certificate(kappa,delta,omega=Q(1),a=Q(3,4),b=Q(5,4)):
 if not (0<a<b):raise ValueError('positive ordered radii required')
 low=a*(1-a*a)+min(Q(0),kappa)*a-abs(delta)
 high=b*(1-b*b)+max(Q(0),kappa)*b+abs(delta)
 speed_min=omega-abs(kappa)/2-abs(delta)/a
 speed_max=omega+abs(kappa)/2+abs(delta)/a
 div_max=2+kappa-4*a*a
 exists=low>0 and high<0 and speed_min>0
 unique=exists and div_max<0
 ans={'kappa':kappa,'delta':delta,'omega':omega,'inner_radius':a,'outer_radius':b,
 'radial_inner_lower':low,'radial_outer_upper':high,'angular_lower':speed_min,'angular_upper':speed_max,'divergence_upper':div_max,
 'existence':'CERTIFIED' if exists else 'UNVERIFIED',
 'unique_cycle_and_all_annulus_distance_convergence':'CERTIFIED' if unique else 'UNVERIFIED',
 'period_pi_coefficients':None,'negative_log_multiplier_pi_coefficient':None}
 if exists:
  ans['period_pi_coefficients']=[2/speed_max,2/speed_min]
 if unique:
  ans['negative_log_multiplier_pi_coefficient']=(-div_max)*2/speed_max
 return ans

def verify(rec):
 fields=['kappa','delta','omega','inner_radius','outer_radius']
 if any(k not in rec or not isinstance(rec[k],Q) for k in fields):return False
 try:wanted=certificate(*(rec[k] for k in fields))
 except (ValueError,ZeroDivisionError):return False
 return rec==wanted

models=[certificate(Q(0),Q(1,10)),certificate(Q(1,8),Q(1,10)),certificate(Q(1,3),Q(1,10))]
expected=[(Q(73,320),Q(-193,320),Q(13,15),Q(17,15),Q(-1,4)),
 (Q(73,320),Q(-143,320),Q(193,240),Q(287,240),Q(-1,8)),
 (Q(73,320),Q(-179,960),Q(7,10),Q(13,10),Q(1,12))]
for r,ans in zip(models,expected):
 check(verify(r),'valid_model_certificate')
 check(tuple(r[k] for k in ['radial_inner_lower','radial_outer_upper','angular_lower','angular_upper','divergence_upper'])==ans,'displayed_bound_values')
check(models[0]['period_pi_coefficients']==[Q(30,17),Q(30,13)],'base_period_bounds')
check(models[1]['period_pi_coefficients']==[Q(480,287),Q(480,193)],'anisotropic_period_bounds')
check(models[2]['period_pi_coefficients']==[Q(20,13),Q(20,7)] and models[2]['negative_log_multiplier_pi_coefficient'] is None,'period_bounds_remain_without_uniqueness')
check(models[0]['negative_log_multiplier_pi_coefficient']==Q(15,34),'base_multiplier_exponent')
check(models[1]['negative_log_multiplier_pi_coefficient']==Q(60,287),'anisotropic_multiplier_exponent')
check(models[2]['existence']=='CERTIFIED' and models[2]['unique_cycle_and_all_annulus_distance_convergence']=='UNVERIFIED','existence_does_not_silently_upgrade_to_uniqueness')
# pi > 3 and exp(z) > 1+z+z^2/2 for z>0 are proved/used analytically.
for coefficient,bound in [(Q(15,34),Q(1,3)),(Q(60,287),Q(3,5))]:
 z=3*coefficient;partial=1+z+z*z/2
 check(partial>1/bound,'rational_strict_multiplier_bound')
# Deliberately altered claims or bounds must not pass the same verifier.
for key in ['radial_inner_lower','radial_outer_upper','angular_lower','angular_upper','divergence_upper']:
 bad=copy.deepcopy(models[1]);bad[key]+=Q(1,1000);check(not verify(bad),'reject_tampered_rational_bound')
bad=copy.deepcopy(models[2]);bad['unique_cycle_and_all_annulus_distance_convergence']='CERTIFIED';check(not verify(bad),'reject_unproved_uniqueness')
bad=copy.deepcopy(models[1]);bad['period_pi_coefficients'][0]=Q(30,17);check(not verify(bad),'reject_stale_period_after_model_change')
bad=copy.deepcopy(models[1]);bad['inner_radius']=Q(0);check(not verify(bad),'reject_zero_radius')
bad=copy.deepcopy(models[1]);bad['dimension']=3;check(not verify(bad),'reject_extra_model_field')
failed=certificate(Q(0),Q(1,2));check(failed['existence']=='UNVERIFIED','failure_of_chosen_annulus_is_not_nonexistence')
# Finite parameter tests certify polynomial implementation, not a universal
# inference from sampling. The trigonometric envelopes are proved in the text.
for kappa,delta,omega in product([Q(-1,8),Q(0),Q(1,8),Q(1,3)],[Q(-1,10),Q(0),Q(1,10)],[Q(1),Q(3,2)]):
 P,R=identities(kappa,delta,omega);r=certificate(kappa,delta,omega)
 for j in range(-12,13):
  t=Q(j,7);c=(1-t*t)/(1+t*t);s=2*t/(1+t*t)
  check(c*c+s*s==1,'rational_circle_point')
  check(abs(c*s)<=Q(1,2),'sampled_product_envelope')
  for radius in [Q(3,4),Q(1),Q(5,4)]:
   x=radius*c;y=radius*s;p=value(P,x,y);v=value(R,x,y)
   rd=(x*p+y*v)/radius;ang=(x*v-y*p)/(radius*radius)
   check(rd==radius*(1-radius*radius)+kappa*radius*c*c+delta*c,'sampled_radial_identity')
   check(ang==omega-kappa*s*c-delta*s/radius,'sampled_angular_identity')
   check(r['angular_lower']<=ang<=r['angular_upper'],'sampled_angular_envelope')
   if radius==Q(3,4):check(rd>=r['radial_inner_lower'],'sampled_inner_envelope')
   if radius==Q(5,4):check(rd<=r['radial_outer_upper'],'sampled_outer_envelope')
# Weighted punctured-plane divergence: multiply by the positive s^2.
P,R=vector_field(Q(0),Q(0));div=add(derivative(P,0),derivative(R,1))
weighted_num=add(mul(S,div),scale(-2,add(mul(X,P),mul(Y,R))))
check(weighted_num==scale(-2,mul(S,S)),'symbolic_punctured_Dulac_minus_two')
for a,b in [(Q(1,2),Q(1)),(Q(1),Q(3,2)),(Q(2,3),Q(4,3))]:
 volume_pi=-2*(b*b-a*a);outer_pi=2*(1-b*b);inner_pi=-2*(1-a*a)
 check(volume_pi==outer_pi+inner_pi,'annulus_two_boundary_flux_identity')
# General competitive field, with parameters instantiated independently.
for aa,bb,cc,dd,ee,ff in [(Q(2),Q(1),Q(1,2),Q(3),Q(1,3),Q(1)),(Q(-2),Q(3),Q(-1),Q(4),Q(5),Q(2))]:
 P=mul(X,add(scale(aa,ONE),scale(-bb,X),scale(-cc,Y)))
 R=mul(Y,add(scale(dd,ONE),scale(-ee,X),scale(-ff,Y)))
 numerator=add(mul(mul(X,Y),add(derivative(P,0),derivative(R,1))),scale(-1,mul(Y,P)),scale(-1,mul(X,R)))
 target=add(scale(-bb,mul(mul(X,X),Y)),scale(-ff,mul(X,mul(Y,Y))))
 check(numerator==target,'symbolic_competition_weighted_divergence')
 for x,y in product([Q(1,100),Q(1,3),Q(2),Q(11)],repeat=2):check(-bb/y-ff/x<0,'positive_quadrant_sign_examples')
# Plateau damping is C1 at s=1, but has many periodic circles inside.
chi=lambda s:Q(0) if s<=1 else (s-1)**2
chip=lambda s:Q(0) if s<=1 else 2*(s-1)
check(chi(Q(1))==chip(Q(1))==0,'plateau_C1_seam_values')
for s in [Q(j,8) for j in range(25)]:
 div=-2*chi(s)-2*s*chip(s)
 check(div<=0 and (div==0)==(s<=1),'nonpositive_but_zero_open_disk')
 if s>1:check(div==-2*(s-1)*(3*s-1),'outer_plateau_formula')
for r in [Q(1,4),Q(1,2),Q(3,4),Q(1)]:check(chi(r*r)==0,'explicit_nontrivial_periodic_circle')
# The omega example has exact squared radial distance z'= -2 z^2.
s0=Q(1,4);last=s0*s0
for t in [Q(0),Q(1),Q(10),Q(100),Q(4992),Q(10000)]:
 z=s0*s0/(1+2*s0*s0*t);der=-2*s0**4/(1+2*s0*s0*t)**2
 check(der==-2*z*z and 0<z<=last,'exact_omega_radial_square_identity');last=z
check(s0*s0/(1+2*s0*s0*4992)==Q(1,10000),'omega_radial_tolerance_time')

def clean(x):
 if isinstance(x,Q):return str(x)
 if isinstance(x,dict):return {str(k):clean(v) for k,v in x.items()}
 if isinstance(x,(list,tuple)):return [clean(v) for v in x]
 return x
out={'status':'PASS','checks':sum(checks.values()),'groups':dict(checks),'certificates':models,
 'strict_multiplier_rational_upper_bounds':['1/3','3/5',None],
 'omega_squared_distance_tolerance':{'initial_radius':'3/4','time':4992,'distance':'1/100'},
 'scope':['General existence uses compactness, uniqueness and the planar theorem proved in the pages.','Sampled rational circle points check arithmetic only; global sine envelopes are analytic.','Dulac exclusions require their domain and regularity hypotheses.','UNVERIFIED declines the chosen certificate and is not a nonexistence statement.']}
ap=argparse.ArgumentParser();ap.add_argument('--output',type=Path);args=ap.parse_args();txt=json.dumps(clean(out),ensure_ascii=False,indent=2)+'\n'
if args.output:args.output.parent.mkdir(parents=True,exist_ok=True);args.output.write_text(txt)
print(txt,end='')