#!/usr/bin/env python3
"""Exact arithmetic verification of the improved count comparison, gamma4<=23/53."""
from fractions import Fraction as F
N=6

def zero(): return [[F(0) for _ in range(N)] for _ in range(N)]
def mat(diag, off):
 a=zero()
 for i,v in enumerate(diag): a[i][i]=F(v)
 for i,j,v in off: a[i][j]=a[j][i]=F(v)
 return a

def outer(a,b): return [[x*y for y in b] for x in a]
def add(*terms): return [[sum(c*a[i][j] for c,a in terms) for j in range(N)] for i in range(N)]
def quad(a,v): return sum(v[i]*a[i][j]*v[j] for i in range(N) for j in range(N))

def ldl_psd(a):
 l=zero(); d=[F(0)]*N
 for j in range(N):
  l[j][j]=1
  d[j]=a[j][j]-sum(l[j][k]**2*d[k] for k in range(j))
  assert d[j]>=0, (j,d[j])
  for i in range(j+1,N):
   r=a[i][j]-sum(l[i][k]*l[j][k]*d[k] for k in range(j))
   if d[j]: l[i][j]=r/d[j]
   else: assert r==0, (i,j,r)
 assert all(a[i][j]==sum(l[i][k]*d[k]*l[j][k] for k in range(N)) for i in range(N) for j in range(N))
 return d

h0=mat([F(3,8),1,F(1,2),1,F(3,2),F(-1,2)],[(0,4,F(-1,4))])
hy=mat([F(1,4),1,F(1,2),1,1,3],[(0,5,F(-1,2)),(1,5,F(3,2)),(2,5,F(-1,2))])
hz=mat([F(1,4),1,F(1,2),1,1,3],[(0,5,F(-1,2)),(1,5,1)])
ell=[F(1),F(1),F(1),F(1),F(0),F(0)]
c=[F(1,2),F(0),F(0),F(0),F(-1),F(-1)]
s=[F(-1,2),F(0),F(0),F(0),F(1),F(2)]
n=add((F(1,2),outer(c,s)),(F(1,2),outer(s,c)))
for name,h,k in [('y',hy,F(1,5)),('z',hz,F(0))]:
 a=add((F(4,5),h0),(F(1,5),h),(F(-7,53),outer(ell,ell)),(-k,n))
 d=ldl_psd(a)
 print(name,'exact LDL pivots:',[str(x) for x in d])

v=[F(i,53) for i in (26,4,16,7,3,10)]
assert sum(v[:4])==1 and v[1]<v[2]
assert v[5]==v[0]/2-v[4]
assert 1-quad(h0,v)==1-quad(hy,v)==F(46,53)
assert quad(n,v)==0
assert F(110,127)<F(46,53)
# Failure of y<4w or x<16w incurs larger old-count quadratic loss.
assert F(17,127)>F(7,53)
assert F(65,484)>F(7,53)
print('PASS: all rational certificates verify gamma(4)<=23/53.')
print('The optimized relaxation witness has count46/53 and coordinates(26,4,16,7,3,10)/53.')