from fractions import Fraction as F
from pathlib import Path
import json
N=7;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=[[F(i==j)for j in range(N)]for i in range(N)]
ell=av(x,y,z,w);q=sub(sc(F(1,2),x),bal);d=av(sc(F(1,2),x),bal);h=sc(2,bal)
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,av(x,z,sc(-3,y))),mm(-3,sy(r,r)))
fT=ms(sy(r,av(sc(F(7,2),q),sc(-1,z),sc(2,y))),mm(F(-23,8),sy(r,r)),mm(F(-1,2),sy(q,q)))
def build(box):
 vl,vh,rl,rh,zl,zh=box
 L=4*(1-rl)/(5*rl-4);beta=2*rh/(1+vl)+F(9,25)/(1+vl)**2
 At=2+2/(1+vl)+L*beta;Ct=1+L;alpha=1/(1+vl)
 ft=ms(fT,mm(F(-1,2),sy(q,sub(q,r))),mm(-L,K))
 forms=[sub(h,sc(vl,q)),sub(sc(vh,q),h),sub(r,sc(rl,q)),sub(sc(rh,q),r),sub(z,sc(zl,q)),sub(sc(zh,q),z),sub(sc(F(3,8),z),y),sub(sc(F(13,20),q),y),sub(sc(2,w),q),av(sc(6,w),sc(-1,r),sc(-1,z)),av(z,sc(-zl/rh,u),sc(F(-3,2),r)),x,y,z,w,q,h,sub(d,r),sub(w,y),sub(y,sc(F(187,1000),q))]
 uforms=[sub(r,u),ell]+forms
 exclusion=ms(sy(z,u),mm(-1,sy(r,z)),mm(F(20,23),sy(r,d)))
 terms=[ms(mm(-1,fR),sy(sub(z,y),u)),ms(mm(-1,ft),mm(F(-1,2),sy(z,u)),mm(F(-1,2)*(zh/rl)**2,sy(u,u))),K,mm(-1,exclusion)]
 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]]+[mm(-1,sy(ell,ell))]
 n=len(terms);B=[]
 for c,vals in[(1,{0:-1,1:-1-At,2:-alpha}),(1,{0:-1,1:-1-Ct,2:-1}),(1,{0:-1,1:-F(15,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','7/8','1','21/16','3']))
 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(3))
  boxes.append(box);v=F(1)
  for i in range(3):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(3))
 expected=F(1)
 for i in range(3):expected*=root[2*i+1]-root[2*i]
 assert vol==expected
 print('PASS improved d>=23Z/20 exclusion:',len(boxes),'exact7D certificates; gamma <=',G,float(G))
if __name__=='__main__':verify(json.loads(Path(__file__).with_name('good-domain-new-certificate.json').read_text()))