Proof of Lemma 1

import time

# function that generates list of possible d's for quadratic integers 
# |a + b*sqrt(-d)| <= m

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*sqrt(-d) with |x| <= m
# and either Re(x)>0 or x in Z and x>0

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

def bmax(m,d):
    if mod(-d,4)==1: 
        return floor(2*m/sqrt(d))
    else: return floor(m/sqrt(d))

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))

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))

# test listing bounded imaginary quadratic integers

m = 3
dlist = makedlist(m)

qilist = [1..floor(m)] # rational integers

# append non-rational integers
for d in dlist:
    K.<s> = QQ.extension(x^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):
            el = a + b*omega
            qilist.append(el[0] + el[1]*sqrt(-d))

print(qilist)
print("number of elements:", len(qilist))

# irreducibility

m = 4

tlist = []

# rational integers:
for a in range(1, m + 1):
    norma = a^2
    if 16/norma in ZZ:
        c = -4/a # is always in Z
        t = - (a+c)
        tlist.append(t)

# non rational integers:
dlist = makedlist(m)

for d in dlist:
    K.<s> = QQ.extension(x^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):
            a_inproof = a + omega*b
            norma = a_inproof.norm()
            if 16/norma in ZZ:
                c_inproof = -4/a_inproof
                if c_inproof.is_integral():
                    t = - (a_inproof + c_inproof)
                    tlist.append(t[0] + t[1]*sqrt(-d))

print(set(tlist))
print("len:", len(set(tlist)))