"""Exact rational audit of the maximum-joint-core packing/swapped saturated-AP certificates.
Combinatorial inputs are documented separately; no floating point is used here.
"""
from fractions import Fraction as F
from pathlib import Path
import json
N=7;ZERO=F(0)
def addv(*v):return[sum(t,ZERO)for t in zip(*v)]
def scale(c,v):return[c*x for x in v]
def sub(a,b):return addv(a,scale(-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=[[F(i==j)for j in range(N)]for i in range(N)]
ell=addv(x,y,z,w);q=sub(scale(F(1,2),x),bal);d=addv(scale(F(1,2),x),bal);h=scale(2,bal);H=sub(d,r)
base=[[F([F(1,4),1,F(1,2),1,1,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(q,q),mm(-1,sy(r,r))))
fR=ms(sy(r,addv(x,z,scale(-3,y))),mm(-3,sy(r,r)))
fT=ms(sy(r,addv(scale(F(7,2),q),scale(-1,z),scale(2,y))),mm(F(-23,8),sy(r,r)),mm(F(-1,2),sy(q,q)))
fA=sy(H,addv(scale(2,z),scale(-1,d),scale(2,q),scale(-4,r)))
fQ=ms(mm(F(1,8),sy(q,q)),mm(-1,sy(q,z)),mm(2,sy(q,y)))

def build(box):
 vl,vh,rl,rh,zl,zh,sl,sh,yl,yh=box
 assert 0<=sl<sh<=1
 small=sl*rl<1-rl
 L=4*(1-rl)/(5*rl-4);beta=1/(1+vl)
 cp_alpha=F(179,1000)
 kc=4/(1-cp_alpha)/(5-4/rl)
 A=3*cp_alpha;Cp=cp_alpha;hc=cp_alpha-F(5,8);ec=F(3,8)
 fCP=ms(fR,mm(cp_alpha,sy(r,addv(z,q,scale(-3,r)))),mm(-cp_alpha**2/2,sy(r,r)))
 At=2+(2-rl)/(1+vl)+L*beta;Ct=1+L;alpha=1/(1+vl)
 epsilon=F(5,64);eta_rank=F(1,10)
 J=ZERO if small else max(ZERO,(F(7,4)-F(23,4)*epsilon+zh-2*yl)/rl-F(3,2)*eta_rank)/(2*sl)
 Ax=1+1/rl+2*(1-rl)/(1+vl)+J+(2+J)*(rh/(1+vl)+L*beta)
 Cx=1+(2+J)*L
 distance=max(sl-F(1,2),F(1,2)-sh,ZERO);gap=distance**2/2
 ft=ms(fT,mm(F(-1,2),sy(q,sub(q,r))),mm(-L,K))
 forms=[sub(h,scale(vl,q)),sub(scale(vh,q),h),sub(r,scale(rl,q)),sub(scale(rh,q),r),sub(z,scale(zl,q)),sub(scale(zh,q),z),
 sub(scale(F(3,8),z),y),sub(scale(F(13,20),q),y),sub(scale(2,w),q),addv(scale(6,w),scale(-1,r),scale(-1,z)),
 addv(z,scale(-zl/rh,u),scale(F(-3,2),r)),sub(scale(F(23,20),sub(z,scale(zl/rh,u))),d),x,y,z,w,q,h,H,sub(w,y),sub(y,scale(yl,q)),sub(scale(yh,q),y)]
 uforms=[sub(r,u),ell]+forms
 joint=ms(mm(F(-1,2),sy(z,u)),mm(F(-1,2)*(zh/rl)**2,sy(u,u)))
 terms=[ms(mm(-1,fR),sy(sub(z,y),u)),ms(mm(-1,ft),joint),
  ms(mm(-1,fCP),mm(-hc,sy(z,u))),
  ms(mm(-1,ft),mm(gap,sy(r,r)),mm(-(2-sl),sy(z,u))),
  ms(mm(-1,fQ),mm(-eta_rank*(sh-2+eta_rank/2),sy(r,r)),mm(-F(3,2)*eta_rank,sy(q,r)),mm(-eta_rank,sy(z,u)),mm(-F(23,8)*epsilon**2,sy(q,q)),mm((2+J)*L,K),mm(-1/(2*rl),sy(z,u)),mm(F(-1,2)*(zh/rl)**2,sy(u,u))),K]
 terms +=[mm(-1,sy(ell,a))for a in forms]+[mm(-1,sy(u,a))for a in uforms]
 terms +=[mm(-1,sy(forms[i],forms[i+1]))for i in[0,2,4]]+[mm(-1,sy(q,a))for a in forms[:6]+forms[-2:]]+[mm(-1,sy(forms[-2],forms[-1])),mm(-1,sy(ell,ell))]
 # Each affine budget is constant + coefficients dot parameters.
 n=len(terms);B=[]
 def row(c,vals):
  a=[ZERO]*n
  for i,v in vals.items():a[i]=F(v)
  B.append((F(c),a))
 # Parameters b,c,e,dd,xx,mu; p=1-b-c-e-dd-xx.
 row(1,{0:-1,1:-1-At,2:-1-A,3:-1-At,4:-1-Ax,5:-alpha})
 row(1,{0:-1,1:-1-Ct,2:-1-Cp,3:-1-Ct,4:-1-Cx,5:-1})
 row(1,{0:-1,1:-1,2:-1,3:-1,4:-2-J})
 row(1,{0:-1,1:-F(15,2),2:-1,3:-5,4:-F(39,4)-F(7,2)*eta_rank})
 row(0,{3:1,4:-J})
 row(1,{0:-1,1:-1,2:-1-ec,3:-1,4:-1})
 if small:row(0,{4:-1})
 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','7/8','1','21/16','3','0','1','187/1000','13/20']))
 boxes=[];vol=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[:-1])>=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(5))
  boxes.append(box);v=F(1)
  for i in range(5):v*=box[2*i+1]-box[2*i]
  vol+=v
 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(5))
 expected=F(1)
 for i in range(5):expected*=root[2*i+1]-root[2*i]
 assert vol==expected
 print('PASS two-orientation deficient CP large:',len(boxes),'exact7D certificates; gamma <=',G,float(G))
if __name__=='__main__':verify(json.loads(Path(__file__).with_name('cp-large-certificate.json').read_text()))