def is_powerof(a,b): # checks if a is a power of b
    x = a
    while mod(x,b) == 0:
        x = x/b
    if x == 1:
        return True
    else:
        return False
    
def remove3(a): # computes a/3^(v_3(a))
    x = a
    while mod(x,3) == 0:
        x = x/3
    return x

# LLL reduction and finishing proof for plus sign case, even x

xmax = 297
M = 8.11 * 10^27

C0 = 10^int(4*log(M)/log(10)) # approx. M^4 but with full precision
Cmax = C0 * 10^50

Uxmax = 1
C = C0 #always keeps increasing, never goes back down! This just works better

for x in range(4, xmax + 1, 2):
    lowerBound = 1
    done = false
    while not done:
        C = C*10
        A = Matrix([[1, 0, 0, round(C*log(3))],
                    [0, 1, 0, round(C*log(2^(x+1) + 1))],
                    [0, 0, 1, round(C*log(2))],
                    [0, 0, 0, round(C*log(2^x + 1))]])

        B = A.LLL()
        Bstar, mu = B.gram_schmidt()
        c = min([norm(N(b)) for b in Bstar])
        if c^2 > 7*M^2:
            lowerBound = 1/C * (sqrt(c^2 - 3*M^2) - 2*M)
            done = True
        elif C == Cmax:
            print('did not work for x =', x)
            break
    Ux = - log(lowerBound)
    #print('x =', x, 'Ux = ', Ux)
    Uxmax = max(Uxmax, Ux)
    rhsfactor = remove3(2^(x+1) + 1)
    for z in range(1, floor((Ux+1)/(x-0.2)) + 1):
        for y in range (1, floor(Ux + 1 - z*(x-0.2)) + 1):
            qxyz = remove3(2^y * (2^x + 1)^z -1)
            if is_powerof(qxyz, rhsfactor):
                if not ([y,z] == [1,1]):
                    print('[x,y,z] =', [x,y,z])
print('Uxmax =', Uxmax)

# LLL reduction and finishing proof for minus sign case, odd x

xmax = 297
M = 8.11 * 10^27

C0 = 10^int(4*log(M)/log(10)) # approx. M^4 but with full precision
Cmax = C0 * 10^50

Uxmax = 1
C = C0 #always keeps increasing, never goes back down! This just works better

for x in range(3, xmax + 1, 2):
    done = false
    while not done:
        C = C*10
        A = Matrix([[1, 0, 0, round(C*log(3))],
                    [0, 1, 0, round(C*log(2^(x+1) - 1))],
                    [0, 0, 1, round(C*log(2))],
                    [0, 0, 0, round(C*log(2^x - 1))]])

        B = A.LLL()
        Bstar, mu = B.gram_schmidt()
        c = min([norm(N(b)) for b in Bstar])
        if c^2 > 7*M^2:
            lowerBound = 1/C * (sqrt(c^2 - 3*M^2) - 2*M)
            done = True
        elif C == Cmax:
            print('did not work for x =', x)
            break
    Ux = - log(lowerBound)
    #print('x =', x, 'Ux = ', Ux)
    Uxmax = max(Uxmax, Ux)
    rhsfactor = remove3(2^(x+1) - 1)
    for z in range(1, floor((Ux+1)/(x-0.2)) + 1):
        for y in range (1, floor(Ux + 1 - z*(x-0.2)) + 1):
            qxyz = remove3(2^y * (2^x - 1)^z + 1)
            if is_powerof(qxyz, rhsfactor):
                if not ([y,z] == [1,1] or [y,z] == [x+2,1]):
                    print('[x,y,z] =', [x,y,z])
print('Uxmax =', Uxmax)