Determinants over cyclotomic fields are broken?

savecancel

I launch the code below which generates two unimodular matrices L and U until it thinks that norm(det(L*U)) is not in [-1,1] which is probably a bug:

def seek_bug(K,n, rng=[-3,3], density=1.0):
    print("Searching")
    z=K.gen()
    OK = K.fraction_field().ring_of_integers()

    units = [z^((1-(i))/2) * (1-z^i)/(1-z) for i in list( range(1,K.degree(),2) ) ] if K.degree()>2 else [1,z]   
    assert( all([OK(u).is_unit() for u in units]) ), "Wrong units!"

    L = Matrix.identity(K,n)
    U = Matrix.identity(K,n)
    for i in range(n):
        for j in range(n):
            if i<j:
                if uniform(0,1)<density:
                    L[i,j]=OK.random_element(rng[0],rng[1])
            elif j<i:
                if uniform(0,1)<density:
                    U[n-j-1,n-i-1]=OK.random_element(rng[0],rng[1])
            else:
                U[i,i]*=prod([ units[randrange(len(units))] for i in range(6) ])
                L[i,j]*=prod([ units[randrange(len(units))] for i in range(6) ])
    B = L*U
    if not( norm(det(B)) in [-1,1] ):
        return L, U, B
    return None, None, None

n=2
K.<z> = CyclotomicField(16)
L, U, B = None, None, None

while L is None:
    L, U, B = seek_bug(K,n)
print(norm(det(L)), norm(det(U)))
assert det(L)*det(U) == det(L*U), "This should work... But it doesn't"

The last assertion fails which implies det(L)det(U) != det(LU) (Note: even without taking the norm!), but this should not be the thing. It fails in power-of-two cyclotomics. Did I make a mistake, or are the determinants in sage are broken? Checked this on 3 machines with sage 9.4, 9.5 and 9.6: it fails everywhere.