print("Starting nonsingular computation.")

quadcols = { "SW" : "yellow", "SE" : "orange", "NE" : "purple" , "NW" : "cyan" }

def procQuadrant(qu, nscutoff, alst):
    col = quadcols[qu]
    print("Starting quadrant "+qu+".")
    Q2 = Qlst[qu];
    minxexp = minxexps[qu]
    G1 = G1lst[qu]
    G1low = trunc(G1.subs(A1y=trunc(A1yS,y,1),A2y=trunc(A2yS,y,1),A3y=trunc(A3yS,y,1),A4y=trunc(A4yS,y,1),Lpy=trunc(LpyS,y,1)),y,1)
    G1main = G1 - G1low
    tbG1main = termbound(G1main)
    erbG1main = varerrorbound(1/2,1/2,minxexp,minxexp,nscutoff,10)
    G1maint = plugin(G1main,nscutoff,nscutoff)
    G1maind = G1maint/(x^minxexp*y^minyexp)
    assert G1maind.denominator() == 1
    erG1main = condense(tbG1main, erbG1main)
    # G1low = G1dlow*x^minxexp where G1dlow has the same termbound as G1low
    tbG1low = termbound(G1low)
    G1lowt = plugin(G1low,nscutoff,1)    # only need degree 1 in y
    erbG1low = R((-1,1))*(1/2)^(nscutoff-minxexp-2-1)   # need to multiply by x^2 
    G1lowd = G1lowt/(x^minxexp)
    assert G1lowd.denominator() == 1
    erG1low = condense(tbG1low, erbG1low)
    # so G1low/x^minxexp = G1lowd + erG1low*x^2
    # so now we have (G0low + G02*y^2 + erG0*y^4)*(G1low + G1maind*y^2*x^minxexp + erG1main*y^2*x^minxexp) - G0low*G1low
    # = y^2*x^minxexp*(G0low + G02*y^2 + erG0*y^4)*(G1maind + erG1main) + G1low*y^2*(G02 + erG0*y^2)
    # = y^2*x^minxexp*(G0low + G02*y^2 + erG0*y^4)*(G1maind + erG1main) + y^2*x^minxexp*(G1lowd + erG1low*x^2)*(G02 + erG0*y^2)
    Grem =  (num(G0low + G02*y^2) + erG0*yR^4)*(num(G1maind) + erG1main ) + (num(G02) + erG0*yR^2)*( num(G1lowd) + erG1low*xR^2 )
    assert Grem.degree(LyR) < 2
    Q2modify = Q2 + (G0low*G1low)
    Q3 = plugin(Q2modify,nscutoff,nscutoff)
    timesofar()
    tb2 = termbound(Q2modify)
    tb2main, tb2extra = split(tb2)
    for xexp in range(minxexp):
        assert Q3.coefficient({x:xexp}) == 0
    for yexp in range(minyexp):
        assert Q3.coefficient({y:yexp}) == 0
    Q3 = Q3/(x^minxexp*y^minyexp)
    assert Q3 == numerator(Q3)
    Q3 = numerator(Q3)
    QR3 = Q3.subs(x=xR,y=yR,Lx=LxR,Ly=LyR,piv=Pi,l2=L2) + Grem
    Q3extra = Q3.coefficient({Ly:2})/y^2
    assert denominator(Q3extra) == 1
    Q3extra = numerator(Q3extra)
    Q3main = Q3 - Q3extra*Ly^2*y^2
    QR3main = Q3main.subs(x=xR,y=yR,Lx=LxR,Ly=LyR,piv=Pi,l2=L2) + Grem
    QR3extra = Q3extra.subs(x=xR,y=yR,Lx=LxR,Ly=LyR,piv=Pi,l2=L2)
    QR3 = QR3main + QR3extra*LyR^2*yR^2
    def procregion(xmin,xmax,ymin,ymax,xctr,yctr,nscutoff2x,nscutoff2y,sg,interv):
        er1 = varerrorbound(R(xmax),R(ymax),minxexp,minyexp,nscutoff,10)
        er1extra = varerrorbound(R(xmax),R(ymax),minxexp,max(minyexp-2,0),nscutoff-2,10-2)
        QR4 = subsHxy(QR3+condense(tb2,er1),xctr+uR,yctr+vR)
        err2 = varerrorbounduv(QR4,R(max(xmax-xctr,xctr-xmin)),R(max(ymax-yctr,yctr-ymin)),nscutoff2x,nscutoff2y,10)
        if (interv == True):
            QR4main = subsHxy(QR3main+condense(tb2main,er1),xctr+uR,yctr+vR)
            QR4extra = subsHxy(QR3extra+condense(tb2extra,er1extra),xctr+uR,yctr+vR)
            err2main = varerrorbounduv(QR4main,R(max(xmax-xctr,xctr-xmin)),R(max(ymax-yctr,yctr-ymin)),nscutoff2x,nscutoff2y,10)
            err2extra = varerrorbounduv(QR4extra,R(max(xmax-xctr,xctr-xmin)),R(max(ymax-yctr,yctr-ymin)),nscutoff2x,nscutoff2y-2,10-2)
            QR5main = clean(sg*QR4main,uR,vR,nscutoff2x,nscutoff2y) + err2main
            QR5extra = clean(sg*QR4extra,uR,vR,nscutoff2x,nscutoff2y-2) + err2extra
            def Q6(x,y):
                return ( evalRdiruv(QR5main,x,y,xctr,yctr,0,(R(x).lower()<R(y).lower())) 
                         + evalRdiruv(QR5extra,x,y,xctr,yctr,0,(R(x).lower()<R(y).lower()))*y*4/eR^2 )
        else:
            QR5 = clean(sg*QR4,uR,vR,nscutoff2x,nscutoff2y)
            def Q6(x,y):
                return evalRdiruv(QR5,x,y,xctr,yctr,0,(R(x).lower()<R(y).lower()))
        assert CP(Q6,xmin,xmax,ymin,ymax,qu,col)
    for a in alst:
        procregion(*a)  # standard python trick to make a tuple the argument
    print("Finished quadrant "+qu+".")


procQuadrant("SW",45,[(0,41/100,42/100,49/100,41/200,46/100,10,3,-1,False),
                      (42/100,1/2,0,41/100,46/100,41/200,3,10,+1,True),
                      (50/101,1/2,41/100,49/100,99/200,45/200,2,10,+1,False)])

procQuadrant("SE",45,[(0,41/100,42/100,49/100,41/200,46/100,10,3,-1,False),
                      (42/100,1/2,0,41/100,46/100,41/200,3,10,+1,True),
                      (50/101,1/2,41/100,49/100,99/200,45/200,2,10,+1,False)])

print("Finished nonsingular computation.")