"""Independent exact-predicate checks for the ten geometry articles.
Run: python foundation-geometry-capstone.py (Python 3 standard library only).
Fractions are exact. Euclidean path lengths alone use floating square roots.
Reference visibility and exhaustive oracles intentionally favor auditability.
"""
from fractions import Fraction as F
from collections import deque, defaultdict
from itertools import combinations
import math, json, random, pathlib
ROOT=pathlib.Path(__file__).resolve().parents[1]
import sys
OUTPUT=pathlib.Path(sys.argv[1]) if len(sys.argv)>1 else pathlib.Path.cwd()/'foundation-geometry-capstone-results.json'
def pt(x):return tuple(F(z) for z in x)
def sub(a,b):return(a[0]-b[0],a[1]-b[1])
def cross(a,b):return a[0]*b[1]-a[1]*b[0]
def orient(a,b,c):return cross(sub(b,a),sub(c,a))
def area2(P):return sum(cross(a,b) for a,b in zip(P,P[1:]+P[:1]))
def onseg(a,b,p):return orient(a,b,p)==0 and min(a[0],b[0])<=p[0]<=max(a[0],b[0]) and min(a[1],b[1])<=p[1]<=max(a[1],b[1])
def inside(P,p):
 odd=False
 for a,b in zip(P,P[1:]+P[:1]):
  if onseg(a,b,p):return True
  if (a[1]>p[1])!=(b[1]>p[1]):
   x=a[0]+(p[1]-a[1])*(b[0]-a[0])/(b[1]-a[1])
   if x>p[0]:odd=not odd
 return odd

def visible(P,a,b):
 if a==b:return inside(P,a)
 v=sub(b,a); tt={F(0),F(1)}
 for c,d in zip(P,P[1:]+P[:1]):
  w=sub(d,c);den=cross(v,w)
  if den:
   t=cross(sub(c,a),w)/den;u=cross(sub(c,a),v)/den
   if 0<=t<=1 and 0<=u<=1:tt.add(t)
  elif cross(sub(c,a),v)==0:
   k=0 if v[0] else 1
   for z in(c,d):
    t=(z[k]-a[k])/v[k]
    if 0<=t<=1:tt.add(t)
 tt=sorted(tt)
 return all(inside(P,(a[0]+(x+y)/2*v[0],a[1]+(x+y)/2*v[1])) for x,y in zip(tt,tt[1:]))
def simplicity_failure(P):
 n=len(P)
 for i,j in combinations(range(n),2):
  if (i-j)%n in (1,n-1):continue
  a,b=P[i],P[(i+1)%n];c,d=P[j],P[(j+1)%n]
  if (orient(a,b,c)*orient(a,b,d)<0 and orient(c,d,a)*orient(c,d,b)<0) or any((onseg(a,b,c),onseg(a,b,d),onseg(c,d,a),onseg(c,d,b))):return [i,j]
 return None

def in_tri(a,b,c,p):return min(orient(a,b,p),orient(b,c,p),orient(c,a,p))>=0

def ears_slow(P):
 assert area2(P)>0
 ids=list(range(len(P)));T=[];history=[];tests=0
 while len(ids)>3:
  found=False
  for j,b in enumerate(ids):
   a=ids[j-1];c=ids[(j+1)%len(ids)];tests+=1
   if orient(P[a],P[b],P[c])<=0:continue
   if any(in_tri(P[a],P[b],P[c],P[x]) for x in ids if x not in(a,b,c)):continue
   assert visible(P,P[a],P[c]);T.append((a,b,c));history.append({'ear':(a,b,c),'remaining':ids.copy()});ids.pop(j);found=True;break
  assert found,('no ear',ids)
 T.append(tuple(ids));assert sum(orient(P[a],P[b],P[c]) for a,b,c in T)==area2(P)
 return T,history,tests

def ears(P):
 """O(n^2) ear-candidate predicate work with linked neighbors and versioned candidates.
 The extra exact visible() assertion at every removal is a reference audit;
 its work is not included in that ear-test bound or an optimized runtime claim.
 """
 import heapq
 n=len(P);prev=[(i-1)%n for i in range(n)];nxt=[(i+1)%n for i in range(n)];alive=[True]*n;version=[0]*n;heap=[];tests=0;T=[];history=[]
 def evaluate(b):
  nonlocal tests
  tests+=1;version[b]+=1;a=prev[b];c=nxt[b]
  good=orient(P[a],P[b],P[c])>0 and not any(alive[x] and x not in(a,b,c) and in_tri(P[a],P[b],P[c],P[x]) for x in range(n))
  if good:heapq.heappush(heap,(b,version[b]))
 for b in range(n):evaluate(b)
 left=n
 while left>3:
  while heap:
   b,v=heapq.heappop(heap)
   if alive[b] and version[b]==v:break
  else:raise AssertionError('no valid ear')
  a,c=prev[b],nxt[b];assert visible(P,P[a],P[c]);T.append((a,b,c));history.append({'ear':(a,b,c),'remaining':[i for i in range(n) if alive[i]]})
  alive[b]=False;nxt[a]=c;prev[c]=a;left-=1;evaluate(a);evaluate(c)
 a=next(i for i in range(n) if alive[i]);T.append((a,nxt[a],nxt[nxt[a]]))
 assert tests==n+2*(n-3)
 assert sum(orient(P[a],P[b],P[c]) for a,b,c in T)==area2(P)
 return T,history,tests

def sleeve(P,T,s,t):
 si=next(i for i,(a,b,c) in enumerate(T) if in_tri(P[a],P[b],P[c],s))
 ti=next(i for i,(a,b,c) in enumerate(T) if in_tri(P[a],P[b],P[c],t))
 edge=defaultdict(list)
 for i,tr in enumerate(T):
  for a,b in zip(tr,tr[1:]+tr[:1]):edge[tuple(sorted((a,b)))].append(i)
 adj=defaultdict(list)
 for e,ts in edge.items():
  if len(ts)==2:
   x,y=ts;adj[x].append((y,e));adj[y].append((x,e))
 assert sum(len(v) for v in adj.values())//2==len(T)-1
 prev={si:None};q=deque([si])
 while q:
  i=q.popleft()
  for j,e in adj[i]:
   if j not in prev:prev[j]=(i,e);q.append(j)
 chain=[];i=ti
 while i!=si:
  a,e=prev[i];chain.append((a,i,e));i=a
 chain.reverse();portals=[]
 for a,b,(i,j) in chain:
  k=next(x for x in T[a] if x not in(i,j))
  if orient(P[i],P[j],P[k])>0:i,j=j,i
  portals.append((P[i],P[j]))
 return portals,[si]+[b for a,b,e in chain]

def funnel(portals,s,t):
 """Two persistent deques; no restart and no old-portal replay."""
 L=deque([s]);R=deque([s]);tail=[s];pushes=0;pops=0;apex_moves=0;history=[]
 def record(label):history.append({'event':label,'tail':tail.copy(),'L':list(L),'R':list(R)})
 def left(v):
  nonlocal pushes,pops,apex_moves
  if v==L[-1]:return
  while len(L)>1 and orient(L[-2],L[-1],v)<=0:L.pop();pops+=1
  if len(L)==1:
   while len(R)>1 and orient(R[0],R[1],v)<0:
    R.popleft();pops+=1;apex_moves+=1
    L[0]=R[0]
    if tail[-1]!=R[0]:tail.append(R[0])
  if v!=L[-1]:L.append(v);pushes+=1
 def right(v):
  nonlocal pushes,pops,apex_moves
  if v==R[-1]:return
  while len(R)>1 and orient(R[-2],R[-1],v)>=0:R.pop();pops+=1
  if len(R)==1:
   while len(L)>1 and orient(L[0],L[1],v)>0:
    L.popleft();pops+=1;apex_moves+=1
    R[0]=L[0]
    if tail[-1]!=L[0]:tail.append(L[0])
  if v!=R[-1]:R.append(v);pushes+=1
 if portals:
  l,r=portals[0];L.append(l);R.append(r);pushes=2;record('initial')
  for i,(l,r) in enumerate(portals[1:],1):
   if l==L[-1]:right(r);record('right '+str(i))
   elif r==R[-1]:left(l);record('left '+str(i))
   else:raise AssertionError('successive portals must share endpoint')
  left(t);record('target')
  out=tail[:-1]+list(L)
 else:out=[s,t]
 assert pops<=pushes, (pops,pushes)
 return out,history,{'pushes':pushes,'pops':pops,'apex_moves':apex_moves,'portals':len(portals),'replayed_portals':0}

def distance(a,b):return math.hypot(float(a[0]-b[0]),float(a[1]-b[1]))
def length(path):return sum(distance(a,b) for a,b in zip(path,path[1:]))
def visibility_shortest(P,s,t):
 Q=[s,t]+P;n=len(Q);D=[math.inf]*n;D[0]=0;prev=[None]*n;seen=set();edge=[]
 for i,j in combinations(range(n),2):
  if visible(P,Q[i],Q[j]):edge.append((i,j,distance(Q[i],Q[j])))
 adj=defaultdict(list)
 for i,j,w in edge:adj[i].append((j,w));adj[j].append((i,w))
 for _ in range(n):
  i=min((k for k in range(n) if k not in seen),key=lambda k:D[k]);seen.add(i)
  for j,w in adj[i]:
   if D[i]+w<D[j]:D[j]=D[i]+w;prev[j]=i
 out=[];i=1
 while i is not None:out.append(Q[i]);i=prev[i]
 return out[::-1],D[1],edge

def obstacle_visible(obstacles,a,b):
 if a==b:return True
 v=sub(b,a);tt={F(0),F(1)}
 for O in obstacles:
  for c,d in zip(O,O[1:]+O[:1]):
   w=sub(d,c);den=cross(v,w)
   if den:
    t=cross(sub(c,a),w)/den;u=cross(sub(c,a),v)/den
    if 0<=t<=1 and 0<=u<=1:tt.add(t)
   elif cross(sub(c,a),v)==0:
    k=0 if v[0] else 1
    for z in(c,d):
     t=(z[k]-a[k])/v[k]
     if 0<=t<=1:tt.add(t)
 tt=sorted(tt)
 for x,y in zip(tt,tt[1:]):
  p=(a[0]+(x+y)/2*v[0],a[1]+(x+y)/2*v[1])
  for O in obstacles:
   if inside(O,p) and not any(onseg(c,d,p) for c,d in zip(O,O[1:]+O[:1])):return False
 return True

def obstacle_path(obstacles,s,t):
 Q=[s,t]+[p for O in obstacles for p in O];N=len(Q);D=[math.inf]*N;D[0]=0;seen=set();prev=[None]*N;E=[]
 for i,j in combinations(range(N),2):
  if obstacle_visible(obstacles,Q[i],Q[j]):E.append((i,j,distance(Q[i],Q[j])))
 for _ in range(N):
  u=min((i for i in range(N) if i not in seen),key=lambda i:D[i]);seen.add(u)
  for i,j,w in E:
   if j==u:i,j=j,i
   if i==u and D[i]+w<D[j]:D[j]=D[i]+w;prev[j]=i
 out=[];u=1
 while u is not None:out.append(Q[u]);u=prev[u]
 return out[::-1],D[1],E

def clip(P,a,b):
 out=[]
 if not P:return out
 for p,q in zip(P,P[1:]+P[:1]):
  u=orient(a,b,p);v=orient(a,b,q)
  if (u>=0)!=(v>=0):
   t=u/(u-v);out.append((p[0]+t*(q[0]-p[0]),p[1]+t*(q[1]-p[1])))
  if v>=0:out.append(q)
 clean=[]
 for p in out:
  if not clean or p!=clean[-1]:clean.append(p)
 if len(clean)>1 and clean[0]==clean[-1]:clean.pop()
 return clean

def intersect(l,m):
 a,b=l;c,d=m;v=sub(b,a);w=sub(d,c);den=cross(v,w);assert den!=0
 t=cross(sub(c,a),w)/den
 return(a[0]+t*v[0],a[1]+t*v[1])
def halfplane_deque(lines):
 # Main theorem contract: no pair of boundary lines is parallel; bounded full-dimensional intersection.
 from functools import cmp_to_key
 def half(v):return 0 if v[1]>0 or (v[1]==0 and v[0]>=0) else 1
 def cmp(l,m):
  v=sub(l[1],l[0]);w=sub(m[1],m[0]);h=half(v)-half(w)
  return h if h else (-1 if cross(v,w)>0 else 1)
 assert all(cross(sub(b,a),sub(d,c)) for (a,b),(c,d) in combinations(lines,2))
 Q=deque();history=[];pops=0
 for h in sorted(lines,key=cmp_to_key(cmp)):
  while len(Q)>1 and orient(*h,intersect(Q[-2],Q[-1]))<0:Q.pop();pops+=1
  while len(Q)>1 and orient(*h,intersect(Q[0],Q[1]))<0:Q.popleft();pops+=1
  Q.append(h);history.append(list(Q))
 while True:
  oldlen=len(Q)
  while len(Q)>2 and orient(*Q[0],intersect(Q[-2],Q[-1]))<0:Q.pop();pops+=1
  while len(Q)>2 and orient(*Q[-1],intersect(Q[0],Q[1]))<0:Q.popleft();pops+=1
  if len(Q)==oldlen:break
 assert len(Q)>=3
 raw=[intersect(Q[i-1],Q[i]) for i in range(len(Q))]
 out=[]
 for p in raw:
  if not out or p!=out[-1]:out.append(p)
 if len(out)>1 and out[0]==out[-1]:out.pop()
 assert all(orient(*h,p)>=0 for h in lines for p in out)
 return out,history,pops

def hull(P):
 P=sorted(set(P))
 if len(P)<3:return P
 def chain(P):
  H=[]
  for p in P:
   while len(H)>1 and orient(H[-2],H[-1],p)<=0:H.pop()
   H.append(p)
  return H
 return chain(P)[:-1]+chain(P[::-1])[:-1]
def minkowski(P,Q):
 def rotate(P):
  j=min(range(len(P)),key=lambda i:(P[i][1],P[i][0]));return P[j:]+P[:j]
 P=rotate(P);Q=rotate(Q);i=j=0;out=[];history=[]
 while i<len(P) or j<len(Q):
  p=P[i%len(P)];q=Q[j%len(Q)];out.append((p[0]+q[0],p[1]+q[1]));history.append((i,j,out[-1]))
  if i==len(P):j+=1;continue
  if j==len(Q):i+=1;continue
  v=sub(P[(i+1)%len(P)],p);w=sub(Q[(j+1)%len(Q)],q);c=cross(v,w)
  half=lambda z:0 if z[1]>0 or (z[1]==0 and z[0]>=0) else 1
  if half(v)<half(w):i+=1
  elif half(v)>half(w):j+=1
  elif c>0:i+=1
  elif c<0:j+=1
  else:i+=1;j+=1
 assert all(orient(out[i-1],out[i],out[(i+1)%len(out)])>0 for i in range(len(out)))
 k=min(range(len(out)),key=lambda i:out[i]);out=out[k:]+out[:k]
 return out,history

def diameter(P):
 n=len(P);j=1;best=F(-1);pairs=[];history=[];advances=0
 def score(i,j):return sum(z*z for z in sub(P[i%n],P[j%n]))
 for i in range(n):
  def ht(k):return orient(P[i],P[(i+1)%n],P[k%n])
  while ht(j+1)>ht(j):j=(j+1)%n;advances+=1
  js=[j]
  if ht(j+1)==ht(j):js.append((j+1)%n)
  for k in js:
   for q in(i,(i+1)%n):
    val=score(q,k)
    if val>best:best=val;pairs=[(q,k)]
    elif val==best:pairs.append((q,k))
  history.append((i,j,js))
 assert advances<=2*n
 return best,pairs,history,advances

# Sweep-line monotone partition. This auditable reference uses a sorted active list
# rather than a balanced tree: same helper transitions, O(n^2) reference time.
def monotone_partition(P):
 n=len(P);S=[(x,y+x/F(1000)) for x,y in P];order=sorted(range(n),key=lambda i:(S[i][1],S[i][0]),reverse=True)
 assert len({p[1] for p in S})==n
 below=lambda a,b:S[a][1]<S[b][1]
 typ={}
 for i in range(n):
  prev=(i-1)%n; nxt=(i+1)%n;c=orient(P[prev],P[i],P[nxt])>0
  if below(prev,i) and below(nxt,i):typ[i]='start' if c else 'split'
  elif below(i,prev) and below(i,nxt):typ[i]='end' if c else 'merge'
  else:typ[i]='regular'
 active={};diags=[];hist=[]
 def connect(i,j):
  if i!=j and (i-j)%n not in(1,n-1):
   assert visible(P,P[i],P[j]);diags.append(tuple(sorted((i,j))))
 def handle_helper(i,e):
  if typ[active[e]]=='merge':connect(i,active[e])
 def leftedge(i):
  y=S[i][1];x=S[i][0];candidates=[]
  for e in active:
   a,b=S[e],S[(e+1)%n]; xx=a[0]+(y-a[1])*(b[0]-a[0])/(b[1]-a[1])
   if xx<x:candidates.append((xx,e))
  assert candidates,(i,active)
  return max(candidates)[1]
 for i in order:
  prev=(i-1)%n;t=typ[i]
  if t=='start':active[i]=i
  elif t=='end':handle_helper(i,prev);del active[prev]
  elif t=='split':e=leftedge(i);connect(i,active[e]);active[e]=i;active[i]=i
  elif t=='merge':handle_helper(i,prev);del active[prev];e=leftedge(i);handle_helper(i,e);active[e]=i
  elif below((i+1)%n,i):handle_helper(i,prev);del active[prev];active[i]=i
  else:e=leftedge(i);handle_helper(i,e);active[e]=i
  hist.append({'vertex':i,'type':t,'helpers':dict(active),'diagonals':diags.copy()})
 faces=[list(range(n))]
 for a,b in diags:
  face=next(f for f in faces if a in f and b in f);faces.remove(face);ia=face.index(a);ib=face.index(b)
  if ia>ib:ia,ib=ib,ia
  faces += [face[ia:ib+1],face[ib:]+face[:ia+1]]
 for face in faces:
  lo=min(range(len(face)),key=lambda j:S[face[j]][1]);hi=max(range(len(face)),key=lambda j:S[face[j]][1])
  for step in(1,-1):
   j=hi
   while j!=lo:
    k=(j+step)%len(face);assert S[face[j]][1]>S[face[k]][1];j=k
 return diags,faces,hist,S

def monotone_triangulate(P,face,S):
 """Linear after merging the two sorted chains; reference constructs order by sorting."""
 n=len(face);top=max(face,key=lambda i:S[i][1]);bottom=min(face,key=lambda i:S[i][1]);pos=face.index(top);left=[]
 while face[pos]!=bottom:left.append(face[pos]);pos=(pos+1)%n
 left.append(bottom);L=set(left);order=sorted(face,key=lambda i:S[i][1],reverse=True);stack=order[:2];T=[];hist=[];push=2;pop=0
 def emit(a,b,c):
  if orient(P[a],P[b],P[c])<0:b,c=c,b
  assert orient(P[a],P[b],P[c])>0;T.append((a,b,c))
 for k in range(2,len(order)-1):
  v=order[k];old=stack.copy()
  if (v in L)!=(stack[-1] in L):
   while len(stack)>1:b=stack.pop();pop+=1;emit(v,b,stack[-1])
   stack=[order[k-1],v];pop+=1;push+=2
  else:
   b=stack.pop();pop+=1
   while stack:
    c=stack[-1];z=orient(P[v],P[b],P[c]);valid=(z<0 if v in L else z>0)
    if not valid:break
    emit(v,b,c);b=stack.pop();pop+=1
   stack += [b,v];push+=2
  hist.append({'vertex':v,'before':old,'after':stack.copy(),'triangles':T.copy()})
 v=order[-1]
 while len(stack)>1:b=stack.pop();pop+=1;emit(v,b,stack[-1])
 assert len(T)==n-2
 assert sum(orient(P[a],P[b],P[c]) for a,b,c in T)==area2([P[i] for i in face])
 return T,hist,{'pushes':push,'pops':pop}

# A square split by a diagonal: DCEL next/prev/twin/face checks before and after deletion.
def dcel_check():
 cycles={'f1':['AB','BC','CA'],'f2':['AC','CD','DA'],'out':['BA','AD','DC','CB']};nexts={};prev={};face={}
 for f,C in cycles.items():
  for i,e in enumerate(C):nexts[e]=C[(i+1)%len(C)];prev[e]=C[i-1];face[e]=f
 before={e:{'twin':e[::-1],'next':nexts[e],'prev':prev[e],'face':face[e]} for e in nexts}
 for e in nexts:assert nexts[prev[e]]==e and prev[nexts[e]]==e
 # Remove CA and AC: splice BC->CD and DA->AB, with reciprocal prev pointers.
 nexts['BC']='CD';prev['CD']='BC';nexts['DA']='AB';prev['AB']='DA'
 for e in ('CA','AC'):del nexts[e];del prev[e];del face[e]
 for e in ('AB','BC','CD','DA'):face[e]='inside'
 for e in nexts:assert nexts[prev[e]]==e and prev[nexts[e]]==e
 return before,{e:{'twin':e[::-1],'next':nexts[e],'prev':prev[e],'face':face[e]} for e in nexts}

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

def run():
 report={};P=list(map(pt,[(0,0),(8,0),(8,6),(6,6),(6,2),(4,2),(4,5),(2,5),(2,6),(0,6)]));s=pt((1,5));t=pt((7,5))
 T,eh,tests=ears(P);portals,ts=sleeve(P,T,s,t);path,fh,counts=funnel(portals,s,t);op,d,edges=visibility_shortest(P,s,t)
 assert abs(length(path)-d)<1e-9,(path,op,length(path),d)
 assert all(visible(P,a,b) for a,b in zip(path,path[1:]))
 diags,faces,mh,S=monotone_partition(P);mt=[];stackhist=[]
 for face in faces:
  tt,hh,cc=monotone_triangulate(P,face,S);mt+=tt;stackhist+=hh
 assert len(mt)==len(P)-2
 report['capstone']={'polygon':P,'start':s,'target':t,'area':area2(P)/2,'ears':T,'ear_history':eh,'ear_tests':tests,'portals':portals,'triangle_path':ts,'funnel_path':path,'funnel_history':fh,'funnel_counts':counts,'visibility_path':op,'length':d,'visibility_edges':len(edges),'monotone_diagonals':diags,'monotone_faces':faces,'partition_history':mh,'monotone_triangles':mt,'stack_history':stackhist}
 C=list(map(pt,[(-1,1),(2,-1),(5,2),(2,4)]));states=[C]
 for a,b in [(pt((0,2)),pt((0,0))),(pt((3,0)),pt((3,2))),(pt((0,0)),pt((3,0))),(pt((3,2)),pt((0,2)))]:C=clip(C,a,b);states.append(C)
 assert area2(C)==F(71,6);report['clipping']={'states':states,'final':C,'area':area2(C)/2}
 H=list(map(pt,[(0,0),(4,0),(5,2),(2,5),(-1,3)]));lines=list(zip(H,H[1:]+H[:1]));lines+=[(pt((-3,-2)),pt((1,-1)))]
 hp,hh,pc=halfplane_deque(lines);assert set(hp)==set(H);report['halfplanes']={'lines':lines,'output':hp,'pops':pc,'history':hh}
 # Pairwise nonparallel lines can share a polygon vertex: remove cyclic duplicates.
 triangle=list(map(pt,[(0,0),(4,0),(0,4)]));tangent_lines=list(zip(triangle,triangle[1:]+triangle[:1]))+[(pt((0,0)),pt((1,-2)))];tangent,_,_=halfplane_deque(tangent_lines)
 assert len(tangent)==3 and set(tangent)==set(triangle) and area2(tangent)==16
 report['halfplanes']['concurrent_boundary_regression']={'lines':tangent_lines,'output':tangent,'distinct_vertices':3}
 C=hull(list(map(pt,[(0,0),(4,0),(6,2),(4,5),(1,6),(-1,3)])));best,pairs,ch,adv=diameter(C);brute=max(sum(z*z for z in sub(a,b)) for a,b in combinations(C,2));assert best==brute
 report['calipers']={'polygon':C,'squared_diameter':best,'pairs':pairs,'history':ch,'advances':adv}
 O=list(map(pt,[(4,2),(7,2),(7,4),(4,4)]));B=list(map(pt,[(0,0),(2,0),(0,1)]));minus=hull([(-x,-y) for x,y in B]);M,mh=minkowski(O,minus);oracle=hull([(a[0]+b[0],a[1]+b[1]) for a in O for b in minus]);assert set(M)==set(oracle)
 report['minkowski']={'obstacle':O,'robot':B,'negative_robot':minus,'sum':M,'merge_history':mh}
 O1=list(map(pt,[(2,-3),(3,-3),(3,1),(2,1)]));O2=list(map(pt,[(5,-1),(6,-1),(6,3),(5,3)]));op,od,oe=obstacle_path([O1,O2],pt((0,0)),pt((8,0)))
 assert abs(od-(2*math.sqrt(5)+2+2*math.sqrt(2)))<1e-10
 report['obstacle_visibility']={'obstacles':[O1,O2],'path':op,'length':od,'edge_count':len(oe)}
 # Exact clipping boundary cases: singleton, segment, and empty output.
 def rectangle_clip(poly,lo,hi,bottom,top):
  for aa,bb in [(pt((lo,top)),pt((lo,bottom))),(pt((hi,bottom)),pt((hi,top))),(pt((lo,bottom)),pt((hi,bottom))),(pt((hi,top)),pt((lo,top)))]:poly=clip(poly,aa,bb)
  return poly
 raw=list(map(pt,[(-1,1),(2,-1),(5,2),(2,4)]))
 singleton=rectangle_clip(raw,5,6,1,3);empty=rectangle_clip(raw,6,7,1,3)
 segment=rectangle_clip(list(map(pt,[(0,0),(1,0),(1,1),(0,1)])),1,2,0,1)
 assert singleton==[pt((5,2))] and empty==[] and set(segment)=={pt((1,0)),pt((1,1))}
 report['clipping']['degenerate_cases']={'singleton':singleton,'segment':segment,'empty':empty}
 before,after=dcel_check();report['dcel']={'before':before,'after':after}
 rng=random.Random(1931);ok=0;maxcount=0;skipped=[]
 for it in range(100):
  n=rng.randrange(6,19);angles=sorted(rng.random()*2*math.pi for _ in range(n));Q=[pt((round((r:=rng.randrange(40,100))*math.cos(a)),round(r*math.sin(a)))) for a in angles]
  if len(set(Q))!=n or area2(Q)<=0:skipped.append({'case':it,'reason':'duplicate coordinates or nonpositive orientation'});continue
  bad=simplicity_failure(Q)
  if bad is not None:skipped.append({'case':it,'reason':'non-simple generated polygon','intersecting_edges':bad});continue
  tt,_,_=ears(Q)
  si=rng.randrange(len(tt));ti=rng.randrange(len(tt));ss=tuple(sum(Q[j][k] for j in tt[si])/3 for k in(0,1));tx=tuple(sum(Q[j][k] for j in tt[ti])/3 for k in(0,1))
  portals,_=sleeve(Q,tt,ss,tx);p,_,cnt=funnel(portals,ss,tx);_,dd,_=visibility_shortest(Q,ss,tx)
  assert abs(length(p)-dd)<1e-8,('funnel',it,Q,ss,tx,p,dd,cnt)
  assert all(visible(Q,a,b) for a,b in zip(p,p[1:]));assert cnt['pops']<=cnt['pushes'];maxcount=max(maxcount,cnt['pushes']+cnt['pops']);ok+=1
  ds,fs,hs,sp=monotone_partition(Q)
  for f in fs:monotone_triangulate(Q,f,sp)
  hp=hull(Q);v,_,_,_=diameter(hp);assert v==max(sum(z*z for z in sub(a,b)) for a,b in combinations(hp,2))
 # Edge-merge sum independently compared with all vertex sums.
 sum_cases=0
 for it in range(40):
  U=hull([pt((rng.randrange(-20,21),rng.randrange(-20,21))) for _ in range(12)])
  V=hull([pt((rng.randrange(-20,21),rng.randrange(-20,21))) for _ in range(12)])
  got,_=minkowski(U,V);want=hull([(a[0]+b[0],a[1]+b[1]) for a in U for b in V]);assert set(got)==set(want);sum_cases+=1
 report['random_checks']={'seed':1931,'polygon_cases':ok,'generated_cases':100,'skipped':skipped,'funnel_vs_visibility':'passed','partition_and_stack_area_counts':'passed','calipers_vs_all_pairs':'passed','max_funnel_push_plus_pop':maxcount,'no_portal_replay':True,'minkowski_vs_all_vertex_sums':sum_cases}
 OUTPUT.parent.mkdir(parents=True,exist_ok=True)
 OUTPUT.write_text(json.dumps(encode(report),ensure_ascii=False,indent=2)+'\n');print(json.dumps(encode({'capstone_path':path,'length':d,'funnel_counts':counts,'partition':diags,'tests':report['random_checks']}),ensure_ascii=False,indent=2))
if __name__=='__main__':run()