Revision history [back]

Here is my take to switch your code to work over cyclotomic field as I explained in the comments:

r = 2

FF.<w> = CyclotomicField(4*r)           # we have q = w^4
KK.<L1,L2> = LaurentPolynomialRing(FF)          # L1 stands for w^L, L2 stands for w^(L^2)

base_field = KK

recip_diff = w^2 - w^(-2)     #  = q^(1/2) - q^(-1/2)

def quantum(num,Lcoef=0):
    return (w^(2*num)*L1^Lcoef - w^(-2*num)*L1^(-Lcoef)) // recip_diff

def quantum_factorial(num):
    if num <= 1:
        return 1
    return quantum(num) * quantum_factorial(num-1)

def qexp(x,r):
    ans = matrix(r^2,r^2,0)
    for n in range(r):
        ans += (w^(n*(n-1)) / quantum_factorial(n)) * x^n
    return ans

E = matrix(base_field, r,r,0)
F = matrix(base_field, r,r,0)
K = matrix(base_field, r,r,0)

for i in range(r):
    if i < r-1:
        E[i,i+1] = quantum(r-(i+1),1)
    if i > 0:
        F[i,i-1] = quantum(i)
    K[i,i] = w^((r-2*i-1)*2) * L1^2

H = matrix(base_field, r^2, r^2, 0)
P = matrix(base_field, r^2, r^2, 0)
for i in range(r):
    for j in range(r):
        k = r*i + j
        H[k,k] = w^((r-2*i-1)*(r-2*j-1)) * L1^(r-2*i-1 + r-2*j-1) * L2                          # = q^((r-2*i-1+L)*(r-2*j-1+L)/4)
        P[k,r*j + i] = 1


R = P * H * qexp(recip_diff * E.tensor_product(F),r)

crossE = E.tensor_product(K) + matrix.identity(r).tensor_product(E)
crossF = F.tensor_product(matrix.identity(r)) + K.inverse().change_ring(KK).tensor_product(F)
crossK = K.tensor_product(K)

# linear algebra does not work over our ring KK, so we temporarily switch to SR
kernel = crossE.change_ring(SR).right_kernel()
eigen_vectors = crossK.change_ring(SR).eigenvectors_right()
new_basis = []
for i in range(r):
    v = kernel.intersection(span(eigen_vectors[i][1], kernel.base_ring())).basis()[0]
    new_basis.append(v)
    for j in range(r - 1):
        new_basis.append(crossF * new_basis[-1])

new_basis = matrix(KK,new_basis).transpose()     # return back to KK
inv_basis = new_basis.inverse()

newR = inv_basis * R * new_basis