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.
