Proof of Lemma 3

import time

# function that generates list of possible d's for quadratic integers 
# |a + b*omgega| <= m
# where omega = sqrt(-d) if -d = 2 or 3 mod 4
# and omega = (1 + sqrt(-d))/2 if -d = 1 mod 4

def makedlist(m):
    dlist=[]
    for d in range(1, floor(4 * m^2 - 1) +1):
        if mod(-d,4)==2 or mod(-d,4)==3:
            if sqrt(d) <= m:
                if is_squarefree(d):
                    dlist.append(d)
        elif mod(-d,4)==1:
            if is_squarefree(d):
                dlist.append(d)
    return dlist

#preparations to list all x = a + b*omega with |x| < m
# and either Re(x)>0; or x in Z and x>0

# define omega for give d
def defomega(d,s): # s = sqrt(-d) inputted
    if mod(-d,4)==1: 
        return (1 + s)/2
    else: return s

# compute upper bound for b depending on m and d
def bmax(m,d):
    if mod(-d,4)==1: 
        return floor(2*m/sqrt(d))
    else: return floor(m/sqrt(d))

# compute lower bound for a depending on m, d, b
def amin(m,d,b):
    if mod(-d,4)==1: 
        return - floor(sqrt(m^2 - d*b^2/4) + b/2)
    else: return - floor(sqrt(m^2 - d*b^2))

# compute upper bound for a depending on m, d, b
def amax(m,d,b):
    if mod(-d,4)==1: 
        return floor(sqrt(m^2 - d*b^2/4) - b/2)
    else: return floor(sqrt(m^2 - d*b^2))

# function that checks if (x,y) can be a solution in field Q(sqrt(-d))
# d = 0 means that x,y are both rational integers
def checksol(d, x, y, s):
    if d == 0:
        tempvar = polygen(QQ)
        K.<s3> = QQ.extension(tempvar^2 + 3)
        mulist = [1, (1 + s3) / 2, (-1 + s3) / 2, -1, (1 - s3) / 2, (-1 - s3) / 2, - I, I]
    elif d == 1:
        mulist = [-1, 1, -s , s]
    elif d == 3:
        mulist = [1, (1 + s) / 2, (-1 + s) / 2, -1, (1 - s) / 2,(-1 - s) / 2]
    else:
        mulist = [-1,1]
    for mu in mulist:
        denom = x^3*y - x*y^3
        if denom != 0:
            t = (x^4 - 6* x^2 * y^2 + y^4 - mu) / denom
            if t.is_integral():
                print("    d=", d, ", t =", t, ", [x, y] =", [x, y])

#computations
mx = 3
                
start=time.time()

# rational x
print("rational x:")
for x in range(1, mx - 1 + 1): # |x| < mx
    print("x =", x)
    m = max(1 + x^4, 6.86)
    print("checking rational y's ...")
    for y in range(x, floor(m) + 1):
        checksol(0, x, y, 1) 
    dlist = makedlist(m)
    print("now non-rational y's:")
    print("dlist =", dlist)
    for d in dlist:
        tempvar = polygen(QQ)
        K.<s> = QQ.extension(tempvar^2 + d)
        omega = defomega(d,s)
        for b in range(1, bmax(m,d) + 1):
            for a in range(amin(m,d,b), amax(m,d,b) + 1):
                y = a + omega*b
                if x^2 <= y.norm():
                    checksol(d,x,y,s)
                    
print("-----------------------------")

# non rational x
print("non-rational x:")
dlist = makedlist(mx)
print("dlist =", dlist) 
for d in dlist:
    print("d =", d)
    tempvar = polygen(QQ)
    K.<s> = QQ.extension(tempvar^2 + d)
    omega = defomega(d,s)
    for bx in range(1, bmax(mx,d) + 1):
        for ax in range(amin(mx,d,bx), amax(mx,d,bx) + 1):
            x = ax + omega*bx
            my = max(1 + x.norm()^2, 6.86)
            for y in range(ceil(sqrt(x.norm())),floor(my) + 1): # rational y's
                checksol(d,x,y,s)
            for by in range(1, bmax(my,d) + 1): # non-rational y's
                for ay in range(amin(my,d,by), amax(my,d,by) + 1):
                    y = ay + omega*by
                    if x.norm() <= y.norm():
                        checksol(d,x,y,s)
                    
print("time:", time.time()-start)