Proof of Lemma 7

##### parameters

j = 3 # type 0 or 3

printd = False # print details
printybound = True # print lower bound for y

kstart = 2

if j == 0: 
    print("type:", j)
    cstart = 1.12
    kmax = 11
    constrouche = 2.71 * 10^16
    exporouche = 31
elif j == 3: 
    print("type:", j)
    cstart = 2.27
    kmax = 11
    constrouche = 9.84 * 10^15
    exporouche = 30
else:
    print("Attention! Type should be 0 or 1")

##### set up    

def roundup(x):
    return n(ceil(100*x)/100,  digits = 4)

S.<s> = PowerSeriesRing(QQ,'s',31)
R.<x> = S[]
T.<t> = LaurentSeriesRing(QQ,'t',31)

f = x^4 - (1/s) * x^3 - 6 * x^2 + (1/s) * x + 1
fprime = x^3 + 3 * (1/s) * x^2 - 12 * x + (1/s)

# get approximation to alpha with Newton's method
x0 = 0
k = 15

for i in range(1,k+1):
    xnew = x0 - f(x0)/fprime(x0)
    x0 = xnew

if j == 0: alpha = x0
elif j == 3: alpha = - (x0+1)/(x0-1)

bs = alpha.truncate(30).power_series()
bt = bs.subs(s=t^-1)
print("B(1/t) =", bt)


##### Steps: 

clist = []
lblist = []

c = cstart

for k in range(kstart, kmax + 1):

    print("\n*** Step", k, "***")
    c0 = c

    pade = bs.pade(k-1,k-1)
    ut = numerator(pade).subs(s=t^-1)
    vt = denominator(pade).subs(s=t^-1)
    denominators = [denominator(x) for x in ut.coefficients()] + [denominator(x) for x in vt.coefficients()]
    ut = ut * lcm(denominators)
    vt = vt * lcm(denominators)
    if printd: print("U(1/t)=", ut), print("V(1/t)=", vt)
    
    # bound U(1/t) - V(1/t) * B(1/t)
    uminusvb = ut - vt * bt
    if printd: print("U(1/t) - V(1/t) * B(1/t) =", uminusvb)
    expmin = 2*k - 1
    if expmin == - max(uminusvb.exponents()):
        if printd: print("2*k - 1 =", expmin)
    else: 
        print("Attention: Something went wrong with exponents in U - VB!")
    expos = uminusvb.exponents()
    coeffs = uminusvb.coefficients()
    c1 = add([abs(coeffs[i]) * 100^(expos[i] + expmin) for i in range(0,len(coeffs))])
    if printd: print('c_1 =', c1.n())
    
    # bound x - alpha*y
    c2 = 8.86 * c0^4
    if printd: print("c_2 =", c2)

    # bound V(1/t)
    expos = vt.exponents()
    coeffs = vt.coefficients()
    c3 = add([abs(coeffs[i]) * 100^(expos[i]) for i in range(0,len(coeffs))])
    if printd: print("c_3 =", c3.n())
    
    # compute c
    c = roundup(c1 + c2*c3*100^(-2*k + 2) + constrouche * c3 * 100^-(exporouche + 1 - 2*k))
    print("c =", c)
    clist.append(c)
    lby = 100^k/c
    lblist.append(lby)
    if printybound:
        print("new lower bound for y =", lby)

    ### Check expression does not vanish
    if printd: print("Check |F_(A(1/t)y, y)|>1:")
    
    denom = (t^(k-1)*vt)^4
    Pt = t^(4*(k-1)) * (ut^4 - t * ut^3 * vt - 6 * ut^2 * vt^2 + t * ut * vt^3 + vt^4)
    if printd: print("P(t) =", Pt)
    if printd: print("denominator = ", denom)
    if max(Pt.exponents()) == 2*k - 2:
        if printd: print("Shape of f(A(1/t)) in paper checks out.")
    else: 
        print("Attention: Shape of f(A(1/t)) is weird!")

    # Constant for lower bound for |P(t)|
    coeffs = Pt.coefficients()
    expos = Pt.exponents()
    clbPt = abs(coeffs[-1]) - add([abs(coeffs[j])*100^(-(2*k - 2) + expos[j]) for j in range(0,len(coeffs)-1)])
    if printd: print("|P(t)| >", clbPt.n(), "* |t| ^ ", 2*k-2)

    # Constant for upper bound for |denominator|
    coeffs = denom.coefficients()
    expos = denom.exponents()
    cubdenom = sum([abs(coeffs[j])*100^(-4*(k-1) + expos[j]) for j in range(0,len(coeffs))])
    if printd: print("|denominator| <", cubdenom.n(), "* |t| ^", 4*(k-1))

    # Overall lower bound
    lb = clbPt * cubdenom^(-1) * c0^(-4) * 100^(2*k - 2)

    if lb > 1:
        if printd: print("|F_t(x,y)| >", lb, "> 1")
    else: 
        print("Attention: Lower bound for |F(x,y)| was too small!")

# produce table
if j == 0:
    kstart = 4
    liststart = 2
if j == 3:
    kstart = 2
    liststart = 0

print("$k$", end=" & $" )
for k in range(kstart, kmax):
    print(k, end="$ & $")
print(kmax, end="$ \\\\ \\hline \n")

print("$c$", end=" & $" )
for k in range(liststart, len(clist)-1):
    print(latex(clist[k].n(digits=4)), end="$ & $")
print(latex(clist[-1].n(digits=4)), end="$ \\\\ \n")


print("$|y| \\geq \\ldots$", end=" & $" )
for k in range(liststart, len(lblist)-1):
    print(latex(lblist[k].n(digits=4)), end="$ & $")
print(latex(lblist[-1].n(digits=4)), "$")