"""Exact halfspace bounds from the existing 23/53 quadratic certificate."""
from fractions import Fraction as F
from pathlib import Path
import runpy

ns = runpy.run_path(str(Path(__file__).resolve().with_name('previous_certificate.py')))
A = ns['add']((F(4,5),ns['h0']),(F(1,5),ns['hy']),(-F(7,53),ns['outer'](ns['ell'],ns['ell'])),(-F(1,5),ns['n']))
v = [F(i,53) for i in (26,4,16,7,3,10)]
W = [[F(0)]*5 for _ in range(6)]
for j in range(3):
    W[j][j]=1
    W[3][j]=-1
W[4][3]=1
W[5][4]=1
K=[[sum(W[i][j]*A[i][k]*W[k][h] for i in range(6) for k in range(6)) for h in range(5)] for j in range(5)]

def solve(mat,rhs):
    a=[list(row)+[x] for row,x in zip(mat,rhs)]
    for j in range(len(a)):
        p=next(i for i in range(j,len(a)) if a[i][j])
        a[j],a[p]=a[p],a[j]
        piv=a[j][j]
        a[j]=[x/piv for x in a[j]]
        for i in range(len(a)):
            if i!=j:
                m=a[i][j]
                a[i]=[x-m*y for x,y in zip(a[i],a[j])]
    return [row[-1] for row in a]

def bound(name,coeff,constant=0):
    coeff=list(map(F,coeff));constant=F(constant)
    center=constant+sum(x*y for x,y in zip(coeff,v))
    assert center>0,(name,center)
    a=[sum(coeff[i]*W[i][j] for i in range(6)) for j in range(5)]
    sol=solve(K,a)
    denom=sum(x*y for x,y in zip(a,sol))
    gap=center**2/denom
    print(name, 'count gap >=',gap, 'gamma gap >=',float(gap/2))
    return gap

if __name__=='__main__':
    checks=[
      ('r>=1/10',[0,0,0,0,0,1],-F(1,10)),
      ('w>=1/10',[0,0,0,1,0,0],-F(1,10)),
      ('z>=1/5',[0,0,1,0,0,0],-F(1,5)),
      ('y<=w',[0,-1,0,1,0,0],0),
      ('2y<=z',[0,-2,1,0,0,0],0),
      ('r<=2w',[0,0,0,2,0,-1],0),
      ('r+z<=6w',[0,0,-1,6,0,-1],0),
      ('2b>=r/3',[0,0,0,0,2,-F(1,3)],0),
      ('2b<=r',[0,0,0,0,-2,1],0),
      ('2b>=1/20',[0,0,0,0,2,0],-F(1,20)),
      ('r>=3/25',[0,0,0,0,0,1],-F(3,25)),
    ]
    for args in checks:bound(*args)