from fractions import Fraction as F
from pathlib import Path
import json
N=8;ZERO=F(0)
def av(*v):return[sum(t,ZERO)for t in zip(*v)]
def sc(c,v):return[c*x for x in v]
def sub(a,b):return av(a,sc(-1,b))
def sy(a,b):return[[(a[i]*b[j]+a[j]*b[i])/2 for j in range(N)]for i in range(N)]
def ms(*v):return[[sum(A[i][j]for A in v)for j in range(N)]for i in range(N)]
def mm(c,A):return[[c*x for x in row]for row in A]
x,y,z,w,bal,r,u,t=[[F(i==j)for j in range(N)]for i in range(N)]
ell=av(x,y,z,w);d=sub(sc(F(1,2),x),bal);q=av(sc(F(1,2),x),bal);h=sc(2,bal)
base=[[F([F(1,4),1,F(1,2),1,1,0,0,0][i])if i==j else ZERO for j in range(N)]for i in range(N)]
K=mm(F(1,2),ms(sy(d,d),mm(-1,sy(r,r))));J=mm(F(1,2),sy(q,sub(d,r)))
fR=ms(sy(r,av(x,z,sc(-3,y))),mm(-3,sy(r,r)))
fD=ms(sy(d,av(x,z,sc(-3,y))),mm(-3,sy(d,d)),mm(-1,sy(d,h)),mm(F(5,8),sy(h,h)))
def build(box):
 vl,vh,rl,rh=box
 P=F(10,9)*(F(53,10)+vh);A=P+vh/(1+vh-F(9,40)/rl);C=P;alpha=1/(1+vl);beta=1/(9*(1+vl))
 forms=[sub(h,sc(vl,d)),sub(sc(vh,d),h),sub(r,sc(rl,d)),sub(sc(rh,d),r),sub(sc(F(3,8),z),y),sub(sc(F(13,20),d),y),sub(sc(3,d),z),sub(d,sc(2,w)),av(sc(6,w),sc(-1,r),sc(-1,z)),sub(w,y),sub(y,sc(F(9,50),d)),sub(r,u),sub(z,t),u,t,x,y,z,w,d,h,sub(w,sc(F(1,10),ell)),sub(sc(F(1,2),ell),d)]
 terms=[ms(mm(-1,fR),sy(sub(z,y),u)),ms(mm(-1,fD),sy(sub(z,y),u),mm(-2,sy(h,t))),K,J]
 terms +=[mm(-1,sy(ell,a))for a in forms]+[mm(-1,sy(u,a))for a in forms]+[mm(-1,sy(t,a))for a in forms]
 terms +=[mm(-1,sy(forms[0],forms[1])),mm(-1,sy(forms[2],forms[3]))]+[mm(-1,sy(d,a))for a in forms[:4]]
 terms +=[ms(sy(z,u),mm(-1,sy(r,t))),mm(-1,sy(ell,ell))]
 n=len(terms);B=[]
 for c,vals in[(1,{0:-1,1:-1-A,2:-alpha,3:-1}),(1,{0:-1,1:-1-C,2:-1,3:-beta}),(1,{0:-1,1:-2}),(1,{0:-F(5,2),1:-F(5,2)})]:
  a=[ZERO]*n
  for i,v in vals.items():a[i]=F(v)
  B.append((F(c),a))
 return terms,B

def pd(A):
 A=[row[:]for row in A]
 for i in range(N):
  pivot=A[i][i];assert pivot>0,(i,pivot)
  for j in range(i+1,N):
   for k in range(j,N):A[j][k]-=A[j][i]*A[k][i]/pivot;A[k][j]=A[j][k]

def verify(data):
 G=F(data['gamma']);root=list(map(F,data['root']));assert root==list(map(F,['3/10','5/4','4/5','1']))
 boxes=[];area=ZERO
 for row in data['cases']:
  box=list(map(F,row['box']));a=list(map(F,row['parameters']));terms,B=build(box)
  assert len(a)==len(terms)and min(a[:-2])>=0 and a[-1]==1-2*G
  assert all(c+sum(x*y for x,y in zip(b,a))>=0 for c,b in B)
  pd(ms(base,*[mm(v,A)for v,A in zip(a,terms)]))
  assert all(root[2*i]<=box[2*i]<box[2*i+1]<=root[2*i+1]for i in range(2))
  boxes.append(box);area+=(box[1]-box[0])*(box[3]-box[2])
 for i,A in enumerate(boxes):
  for B in boxes[i+1:]:assert any(A[2*k+1]<=B[2*k]or B[2*k+1]<=A[2*k]for k in range(2))
 assert area==(root[1]-root[0])*(root[3]-root[2])
 print('PASS early CD top-capacity boundary:',len(boxes),'exact8D certificates; gamma <=',G,float(G))
if __name__=='__main__':verify(json.loads(Path(__file__).with_name('cd-boundary-certificate.json').read_text()))