Q2 = Q8sw;
midcol = "brown"

# In these calculations x and y will be shifted by 1/2, so that
# the center of the square is (0,0). In these coordinates,
# u = x-y and v=-x-y.

Q2x = partial(Q2,x)
Q2y = partial(Q2,y)
Q2u = (Q2x-Q2y)/2
Q2v = -(Q2x+Q2y)/2
print("Computed first partials.")
timesofar()

Q2uv = -(partial(Q2u,x)+partial(Q2u,y))/2

print("Computed second partials.")
timesofar()

#sanity check
# ellderiv(ellderiv(Q2,x),y) - ellderiv(ellderiv(Q2,y),x)

Q2uvx = partial(Q2uv,x)
Q2uvy = partial(Q2uv,y)
Q2uuv = (Q2uvx-Q2uvy)/2
Q2uvv = -(Q2uvx+Q2uvy)/2

print("Computed third partials.")
timesofar()

Q2uuuv = (partial(Q2uuv,x)-partial(Q2uuv,y))/2
Q2uvvx = partial(Q2uvv,x)
Q2uvvy = partial(Q2uvv,y)
Q2uuvv = (Q2uvvx-Q2uvvy)/2
Q2uvvv = -(Q2uvvx+Q2uvvy)/2

print("Computed fourth partials.")
timesofar()

T.<xR,yR,uR,vR,LxR,LyR,LpxR,LpyR,A1xR,A2xR,A3xR,A4xR,A1yR,A2yR,A3yR,A4yR> = PolynomialRing(R,16)

# first: Q2uuuv

midcutoff = 25
Erx = R(-1,1)*10^(-17)*sum(xR^j for j in range(midcutoff+1) )
A1xSh = trunc(A1xS.subs(x=x+1/2).subs(l2=L2),x,midcutoff).subs(x=xR,y=yR) + Erx   # shifted
A2xSh = trunc(A2xS.subs(x=x+1/2).subs(l2=L2),x,midcutoff).subs(x=xR,y=yR) + Erx
A3xSh = trunc(A3xS.subs(x=x+1/2).subs(l2=L2),x,midcutoff).subs(x=xR,y=yR) + Erx
A4xSh = trunc(A4xS.subs(x=x+1/2).subs(l2=L2),x,midcutoff).subs(x=xR,y=yR) + Erx
LpxSh = -L2-sum(2^j*x^j/j for j in range(1,midcutoff+1)).subs(x=xR,y=yR)
LxSh = -L2 - sum(2^j*(-1)^j*x^j/j for j in range(1,midcutoff+1)).subs(x=xR,y=yR)

Ery = R(-1,1)*10^(-17)*sum(yR^j for j in range(midcutoff+1) )
A1ySh = trunc(A1yS.subs(y=y+1/2).subs(l2=L2),y,midcutoff).subs(x=xR,y=yR) + Ery   # shifted
A2ySh = trunc(A2yS.subs(y=y+1/2).subs(l2=L2),y,midcutoff).subs(x=xR,y=yR) + Ery
A3ySh = trunc(A3yS.subs(y=y+1/2).subs(l2=L2),y,midcutoff).subs(x=xR,y=yR) + Ery
A4ySh = trunc(A4yS.subs(y=y+1/2).subs(l2=L2),y,midcutoff).subs(x=xR,y=yR) + Ery
LpySh = -L2-sum(2^j*y^j/j for j in range(1,midcutoff+1)).subs(x=xR,y=yR)
LySh = -L2 - sum(2^j*(-1)^j*y^j/j for j in range(1,midcutoff+1)).subs(x=xR,y=yR)

def checkMiddle(f):
    fnum = f.numerator()
    fden = f.denominator()
    C = fden.content()
    fnum = fnum/C
    fden = fden/C
    tb = termbound(fnum)
    dLx = fnum.degree(Lx)
    dLy = fnum.degree(Ly)
    err1 = sum(tb[i][j]*abs(log(R(4/10)))^(i+j)*8388608*10^(-24)/(4/10*6/10)^8 for j in range(dLy+1) for i in range(dLx+1) )
    # print("err1 = ",err1)
    g = fnum.subs(x=xR+R(1/2),y=yR+R(1/2),Lx=LxR,Ly=LyR,Lpx=LpxR,Lpy=LpyR,A1x=A1xR,A2x=A2xR,A3x=A3xR,A4x=A4xR,A1y=A1yR,A2y=A2yR,A3y=A3yR,A4y=A4yR,piv=Pi,l2=L2)
    g = clean(g,xR,yR,midcutoff,midcutoff)
    g = subsH(g,LpxR,LpxSh)
    g = subsH(g,LpyR,LpySh)
    g = subsH(g,A1xR,A1xSh)
    g = subsH(g,A1yR,A1ySh)
    g = subsH(g,A2xR,A2xSh)
    g = clean(g,xR,yR,midcutoff,midcutoff)
    g = subsH(g,A2yR,A2ySh)
    g = subsH(g,A3xR,A3xSh)
    g = subsH(g,A3yR,A3ySh)
    g = subsH(g,A4xR,A4xSh)
    g = subsH(g,A4yR,A4ySh)
    g = clean(g,xR,yR,midcutoff,midcutoff)
    ab = absbound(g,-1,R(1/10))
    # print("ab = ",ab)
    err2 = sum(ab[i][j]*(i+j)*(2/10)^26/26*abs(log(R(4/10)))^(i+j-1)/(4/10*6/10)^8 for j in range(dLy+1) for i in range(dLx+1) )
    # print("err2 = ",err2)
    g = subsH(g,LxR,LxSh)
    g = subsH(g,LyR,LySh)
    ab2 = absbound(g,midcutoff,R(1/10))   # now there are no log terms so this will just be [[ number ]]
    # print("ab2 = ",ab2)
    err3 = ab2[0][0]/(4/10*6/10)^8;
    # print("err3 = ",err3)
    totalerr = err1*R((-1,1)) + err2*R((-1,1)) + err3*R((-1,1))
    g = clean(g,xR,yR,midcutoff,midcutoff)
    WW = fden.subs(x=xR+R(1/2),y=yR+R(1/2))
    M = 0
    sN = sqhalfside*100
    for j in range(-sN,sN,2):
        # print(j, M)
        for k in range(-sN,sN,2):
            x00 = R((j/100,(j+1)/100))
            y00 = R((k/100,(k+1)/100))
            gg = R(subsH2var(g,xR,yR,x00,y00)/subsH2var(WW,xR,yR,x00,y00)) + totalerr
            M = max(M,R(abs(gg).upper()))
    # print("M = ",M)
    return M

print("Computing bound on first coefficient.")
M1 = checkMiddle(Q2uuuv)
timesofar()
print("Computing bound on second coefficient.")
M2 = checkMiddle(Q2uuvv)
timesofar()
print("Computing bound on third coefficient.")
M3 = checkMiddle(Q2uvvv)
timesofar()

print("Checking positivity in the middle.")

Q2uv00 = evalatcenter(Q2uv)
Q2uuv00 = evalatcenter(Q2uuv)
Q2uvv00 = evalatcenter(Q2uvv)

#p = Q2uv00 + Q2uuv00*u/2 + Q2uvv00*v/2 - M1*u^2/6 - M2*|u*v|/4 - M3*v^2/6
p = Q2uv00 + Q2uuv00*(x-y)/2 + Q2uvv00*(-x-y)/2 - M1*(x-y)^2/6 - M3*(-x-y)^2/6
# subtract absolute value of - M2*(x-y)*(-x-y)/4 in loop

m = minctr

for ml in midlist:
    xl = ml[0]; xu = ml[1]; yl = ml[2]; yu = ml[3];
    for i in range(100*xl,100*xu):
        for j in range(100*yl,100*yu):
            x00 = R((i/100,(i+1)/100))
            y00 = R((j/100,(j+1)/100))
            assert (p.subs(x=x00,y=y00)- M2*abs((x00-y00)*(-x00-y00))/4)  > m
    CPlist.append([[1/2+xl,1/2+xu],[1/2+yl,1/2+yu],midcol])

timesofar()