# This file constructs the kernels, following Section 4.4 from the main paper.
# It also proves inequalities (3a) and (3b) from Section 6.2.

# In several places, we'll use variables divided by I, to avoid having
# to deal with complex numbers. (There's no reason we couldn't in principle,
# but they aren't in the rings defined by setup.sage.)

taudI = Kpx/Kx # tau divided by i, (2.21)
Utau  = 4/piv^2*Kx^2 # (2.20)
Vtau  = 4/piv^2*x*Kx^2 # (2.20)
Wtau  = 4/piv^2*(1-x)*Kx^2 # (2.20)
E2tau = 4/piv^2*Kx*(3*Ex - (2-x)*Kx) # (2.22)
E2Stau = E2tau - 6/(piv*taudI) # (2.8)
E4tau = (Utau^2 + Vtau^2 + Wtau^2)/2 # (2.12)
E6tau = (Utau+Vtau)*(Utau+Wtau)*(Wtau-Vtau)/2 # (2.12)
Deltatau = (Utau*Vtau*Wtau)^2/256 # (2.12)
sqrtDeltatau = Utau*Vtau*Wtau/16 # right after (2.12)
E8tau = E4tau^2
E10tau = E4tau*E6tau
E14tau = E4tau^2*E6tau

# same as above, but for z rather than tau
zdI = Kpy/Ky # z divided by i
Uz  = 4/piv^2*Ky^2
Vz  = 4/piv^2*y*Ky^2
Wz  = 4/piv^2*(1-y)*Ky^2
E2z = 4/piv^2*Ky*(3*Ey - (2-y)*Ky)
E2Sz = E2z - 6/(piv*zdI)
E4z = (Uz^2 + Vz^2 + Wz^2)/2
E6z = (Uz+Vz)*(Uz+Wz)*(Wz-Vz)/2
Deltaz = (Uz*Vz*Wz)^2/256
sqrtDeltaz = Uz*Vz*Wz/16
E8z = E4z^2
E10z = E4z*E6z
E14z = E4z^2*E6z

# From Section 6.2
Ehat00 = -6912*l2*Deltatau - 36*E2tau*E4tau*E6tau + 16*E4tau^3 + 20*E6tau^2 + 108*E4tau*sqrtDeltatau*(Vtau*Lx + Wtau*Lpx)
Ehat01dI = -piv*(6*E2tau^2*E4tau*E6tau - 5*E2tau*E4tau^3 - 7*E2tau*E6tau^2 + 6*E4tau^2*E6tau)
Ehat10dI = 12*piv*(-E2tau*E4tau*E6tau + E6tau^2 + 720*Deltatau)
Ehat11 = 2*piv^2*(E2tau^2*E4tau*E6tau - 2*E2tau*E6tau^2 - 1728*E2tau*Deltatau + E4tau^2*E6tau)
Ehat0 = -taudI*Ehat01dI + Ehat00
Ehat1dI = taudI*Ehat11 + Ehat10dI
p = 101/100
# Inequalities (3a) and (3b):
I3a = Ehat0 - p*Ehat1dI
I3b = -Ehat1dI
# Removing positive factors:
I3a = I3a/(2304*Kx^11/(25*piv^13))
I3b = I3b/(9216*Kx^11/piv^12)
# Switch sign of I3a so that both are supposed to be positive:
I3a = -I3a
# Convert from rational functions to polynomials:
assert denominator(I3a) == 1
I3a = numerator(I3a)
assert denominator(I3b) == 1
I3b = numerator(I3b)

# Basis from Proposition 4.4
phim2taudI = taudI # m2 means -2, note also divided by i
phi0taudI = taudI*E2tau - 3/piv
phi2taudI = taudI*E2tau^2 - 6/piv*E2tau

# Basis from Proposition 4.7 (phit = phi tilde)
phitm2z = -zdI^2 # m2 means -2
phit0z = -zdI^2*E2Sz
phit2z = -zdI^2*E2Sz^2

# Basis from Proposition 4.5
psi0tau = 1
xi2tau = Utau+Wtau
xi2Stau = -Utau-Vtau # from (2.11)
xi4tau = Utau^2+Wtau^2-2*Vtau^2
xi4Stau = Utau^2+Vtau^2-2*Wtau^2 # from (2.11)
psi2tau = xi2tau*Lx + xi2Stau*Lpx 
psi4tau = xi4tau*Lx + xi4Stau*Lpx

# Basis from Proposition 4.8 (psit = psi tilde)
psit0z = Ly
psit2z = Wz
psit4z = Uz^2-Vz^2

# Notation from Section 4.4 (proof of Theorem 4.3)
N = 1728
jtau = 256*(1-x+x^2)^3/(x^2*(1-x)^2) # see Section 2.1.3
jz = 256*(1-y+y^2)^3/(y^2*(1-y)^2) # see Section 2.1.3
jtauz = jtau - jz
fm2tau = E10tau/Deltatau # m2 means -2
f0tau = 1
f2tau = E14tau/Deltatau
f4tau = E4tau
f6tau = E6tau
f8tau = E8tau
f10tau = E10tau
f12tau = Deltatau
f14tau = E14tau
fm2z = E10z/Deltaz # m2 means -2
f0z = 1
f2z = E14z/Deltaz
f4z = E4z
f6z = E6z
f8z = E8z
f10z = E10z
f12z = Deltaz
f14z = E14z

# Upsilon 8+ divided by i (minus sign in front from dividing by i)
Upsilon8pdI = -matrix([[f6tau,0,0],[0,f4tau,0],[0,0,f2tau]])*matrix([[jtauz,0,1],[-2*jtauz,-2,0],[1,0,0]])*matrix([[f12z,0,0],[0,f10z,0],[0,0,f8z]])/(36/piv^2*Deltaz*jtauz)

# Upsilon 8-
Upsilon8m = matrix([[f4tau,0,0],[0,f2tau,0],[0,0,f0tau]])*matrix([[-2*N,0,0],[0,1,-1],[0,N-jtau,jtau]])*matrix([[f10z,0,0],[0,f8z,0],[0,0,f6z]])/(2*N*piv*Deltaz*jtauz)

# K 8+ (minus sign in front from dividing by i)
K8p = -(matrix([phim2taudI,phi0taudI,phi2taudI])*Upsilon8pdI*matrix([[phitm2z],[phit0z],[phit2z]]))[0,0]

# K 8-
K8m = (matrix([psi0tau,psi2tau,psi4tau])*Upsilon8m*matrix([[psit0z],[psit2z],[psit4z]]))[0,0]

K8 = (K8p+K8m)/2
Khat8 = (K8p-K8m)/2

# Upsilon 24+ divided by i (minus sign in front from dividing by i)
Upsilon24pdI = -matrix([[f14tau,0,0],[0,f12tau,0],[0,0,f10tau]])*matrix([[6,0,N/jtauz-6],[-12*jtau+5*N,-2*N/jtauz,12*jtau-7*N],[N/jtauz+6,0,-6]])*matrix([[f4z,0,0],[0,f2z,0],[0,0,f0z]])/(36*N/piv^2*Deltaz)

# Upsilon 24-
Upsilon24m = matrix([[f12tau,0,0],[0,f10tau,0],[0,0,f8tau]])*matrix([[-2*N,-2*N*jtauz,0],[0,jtau+2*jtauz,-1],[0,N-2*jtauz-jtau,1]])*matrix([[f2z,0,0],[0,f0z,0],[0,0,fm2z]])/(2*N*piv*Deltaz*jtauz)

# K 24+ (minus sign in front from dividing by i)
K24p = -(matrix([phim2taudI,phi0taudI,phi2taudI])*Upsilon24pdI*matrix([[phitm2z],[phit0z],[phit2z]]))[0,0]

# K 24-
K24m = (matrix([psi0tau,psi2tau,psi4tau])*Upsilon24m*matrix([[psit0z],[psit2z],[psit4z]]))[0,0]

K24 = (K24p+K24m)/2
Khat24 = (K24p-K24m)/2

# Check that the 8-dimensional kernel agrees with kernels.sage:
Q8 = Khat8
Q8 = Q8*Ky^2/Kx^4*Kpx*piv^4
Q8 = Q8*(x-y)*(1-x-y)*(1-x*y)*(1-(1-x)*y)*(1-x*(1-y))*(1-(1-x)*(1-y))
# Swap x and y with 1-x and 1-y, respectively:
Q8 = Q8.subs(x=1-x,y=1-y,Lx=Lpx,Lpx=Lx,Ly=Lpy,Lpy=Ly,Kx=Kpx,Kpx=Kx,Ex=Epx,Epx=Ex,Ky=Kpy,Kpy=Ky,Ey=Epy,Epy=Ey)
# Apply Legendre identity:
Q8 = Q8.subs(Epx=(piv/2+Kx*Kpx-Ex*Kpx)/Kx,Epy=(piv/2+Ky*Kpy-Ey*Kpy)/Ky)
assert Q8 == Q8mathnum

# Check that the 24-dimensional kernel agrees with kernels.sage:
Q24 = Khat24
Q24 = Q24*6*piv^4*Ky^10/Kx^12*Kpx
Q24 = Q24*(1-y)^2*y^2*(x-y)*(1-x-y)*(1-x*y)*(1-(1-x)*y)*(1-x*(1-y))*(1-(1-x)*(1-y))
# Swap x and y with 1-x and 1-y, respectively:
Q24 = Q24.subs(x=1-x,y=1-y,Lx=Lpx,Lpx=Lx,Ly=Lpy,Lpy=Ly,Kx=Kpx,Kpx=Kx,Ex=Epx,Epx=Ex,Ky=Kpy,Kpy=Ky,Ey=Epy,Epy=Ey)
# Apply Legendre identity:
Q24 = Q24.subs(Epx=(piv/2+Kx*Kpx-Ex*Kpx)/Kx,Epy=(piv/2+Ky*Kpy-Ey*Kpy)/Ky)
assert Q24 == Q24mathnum

# Now for the truncated 24-dimensional kernel.

psitildetrunc = (Ehat0 - zdI*Ehat1dI)/(3456*piv)
psitrunc = (256/y^2 - 256/y + 24 + 4*10^9/970299*y^2)*psitildetrunc;
Khat24trunc = Khat24 - psitrunc

Q24trunc = Khat24trunc
Q24trunc = Q24trunc*piv^14/Kx^11*Ky^12*54
Q24trunc = Q24trunc*(1-y)^2*y^2*(x-y)*(1-x-y)*(1-x*y)*(1-(1-x)*y)*(1-x*(1-y))*(1-(1-x)*(1-y))
# Don't swap x,y with 1-x,1-y in the truncated case.
# Apply Legendre identity:
Q24trunc = Q24trunc.subs(Epx=(piv/2+Kx*Kpx-Ex*Kpx)/Kx,Epy=(piv/2+Ky*Kpy-Ey*Kpy)/Ky)
assert Q24trunc == Q24truncmathnum

Q24truncpsi = -psitrunc
Q24truncpsi = Q24truncpsi*piv^14/Kx^11*Ky^12*54
Q24truncpsi = Q24truncpsi*(1-y)^2*y^2*(x-y)*(1-x-y)*(1-x*y)*(1-(1-x)*y)*(1-x*(1-y))*(1-(1-x)*(1-y))
# Don't	swap x,y with 1-x,1-y in the truncated case.
# Apply Legendre identity:
Q24truncpsi = Q24truncpsi.subs(Epx=(piv/2+Kx*Kpx-Ex*Kpx)/Kx,Epy=(piv/2+Ky*Kpy-Ey*Kpy)/Ky)
assert Q24truncpsi == Q24truncmathnumpsi

Q24truncmain = Khat24
Q24truncmain = Q24truncmain*piv^14/Kx^11*Ky^12*54
Q24truncmain = Q24truncmain*(1-y)^2*y^2*(x-y)*(1-x-y)*(1-x*y)*(1-(1-x)*y)*(1-x*(1-y))*(1-(1-x)*(1-y))
# Don't	swap x,y with 1-x,1-y in the truncated case.
# Apply Legendre identity:
Q24truncmain = Q24truncmain.subs(Epx=(piv/2+Kx*Kpx-Ex*Kpx)/Kx,Epy=(piv/2+Ky*Kpy-Ey*Kpy)/Ky)
assert Q24truncmain == Q24truncmathnummain

# Finally, we prove inequalities (3a) and (3b) from Section 6.2.
# This amounts to showing that I3a and I3b are positive for 0<x<1.
# (Note that we changed the sign of 3a above.)
# This calculation is much like our two-variable computations,
# but easier. We'll use the two-variable functions from setup.sage
# and set y=0. In each case, we use Lemma 4.10 from the code documentation.
# The interval (0,1/2) can be done in a single step, and then
# we switch x with 1-x and break up (0,1/2) into several smaller intervals.

# First we do (3a).
# Convert to A1,A2,A3,A4:
Q1 = I3a.subs(Kpx=A1x+A2x*Lx, Kx=-piv*A2x, Epx=A3x+A4x*Lx, Ex=piv*(A4x-A2x))
# Plug in power series for A1,A2,A3,A4,Lpx:
Q = plugin(Q1,100,0)
# Remove factor of x^4:
assert Q/x^4 == numerator(Q/x^4)
Q = numerator(Q/x^4)
# Compute the associated error bound and add it in.
# (This bound is symmetric about 0, so QR+consense(tb,er)
# and QR-condense(tb,er) are equivalent.)
tb = termbound(Q1)
er = varerrorbound(1/2,0,4,0,100,1)
QR = Q.subs(x=xR,Lx=LxR,piv=Pi,l2=L2)
QR = QR+condense(tb,er)
# We can now evaluate at x=z as follows:
def Q(z):
    return evalRdir(QR,z,0)
# The desired inequality holds on (0,1/2).
assert Q(R(0,1/2)) > 0
# Note that removing the factor of x^4 means we can show
# strict positivity.

# Now we switch x with 1-x and take the same approach:
Q1 = I3a.subs(x=1-x,Lx=Lpx,Lpx=Lx,Kx=Kpx,Kpx=Kx,Ex=Epx,Epx=Ex)
Q1 = Q1.subs(Kpx=A1x+A2x*Lx, Kx=-piv*A2x, Epx=A3x+A4x*Lx, Ex=piv*(A4x-A2x))
Q = plugin(Q1,100,0)
assert Q/x^2 == numerator(Q/x^2)
Q = numerator(Q/x^2)
tb = termbound(Q1)
er = varerrorbound(1/2,0,2,0,100,1)
QR = Q.subs(x=xR,Lx=LxR,piv=Pi,l2=L2)
QR = QR+condense(tb,er)
def Q(z):
    return evalRdir(QR,z,0)
assert Q(R(0,2/7)) > 0
assert Q(R(2/7,3/7)) > 0
assert Q(R(3/7,1/2)) > 0
# The only difference is divisibility by x^2 rather than x^4,
# and checking the inequality on subintervals because
# interval arithmetic on (0,1/2) is insufficient to prove
# positivity.

# Inequality (3b) starts the same way:
Q1 = I3b.subs(Kpx=A1x+A2x*Lx, Kx=-piv*A2x, Epx=A3x+A4x*Lx, Ex=piv*(A4x-A2x))
Q = plugin(Q1,100,0)
assert Q/x^4 == numerator(Q/x^4)
Q = numerator(Q/x^4)
tb = termbound(Q1)
er = varerrorbound(1/2,0,4,0,100,1)
QR = Q.subs(x=xR,Lx=LxR,piv=Pi,l2=L2)
QR = QR+condense(tb,er)
def Q(z):
    return evalRdir(QR,z,0)
assert Q(R(0,1/2)) > 0

# The only difference is that the interval arithmetic
# has trouble with the log(x)^2 terms in the second half
# of the interval when x is near 0. To address this issue,
# we use the same trick as in the two-variable calculations.
# We group the log(x)^2 with a factor of x and use the fact
# that x log(x)^2 < 4/e^2 for 0<x<1.

Q1 = I3b.subs(x=1-x,Lx=Lpx,Lpx=Lx,Kx=Kpx,Kpx=Kx,Ex=Epx,Epx=Ex)
Q1 = Q1.subs(Kpx=A1x+A2x*Lx, Kx=-piv*A2x, Epx=A3x+A4x*Lx, Ex=piv*(A4x-A2x))
Q = plugin(Q1,100,0)
assert Q/x^2 == numerator(Q/x^2)
Q = numerator(Q/x^2)
tb = termbound(Q1)
er = varerrorbound(1/2,0,2,0,100,2)
QR = Q.subs(x=xR,Lx=LxR,piv=Pi,l2=L2)
QR = QR+condense(tb,er)
assert QR.degree(LxR) == 2
assert QR.degree(xR) == 98
assert QR.coefficient({LxR:2,xR:0}) == 0
assert QR.coefficient({LxR:2,xR:1}) == 0
# 0 <= x^2 log(x) <= 4/exp(2) for 0<x<1
QR2 = QR.coefficient({LxR:0}) + QR.coefficient({LxR:1})*LxR + sum([R(0,4/eR^2)*QR.coefficient({LxR:2,xR:i})*xR^(i-1) for i in range(2,100)])
def Q(z):
    return evalRdir(QR2,z,0)
assert Q(R(0,2/7)) > 0
def Q(z):
    return evalRdir(QR,z,0)
assert Q(R(2/7,3/7)) > 0
assert Q(R(3/7,1/2)) > 0