I'm doing some calculations involving matrices whose entries look something like q^L, where q is a root of unity and L is meant to be a generic complex number that I would prefer to leave as a symbol. I would actually prefer to have the final result be displayed in terms of q as well, but that isn't as big of a sticking point.
For small dimensions, if I work over the symbolic ring (SR) and leave q as a symbol too, everything behaves as it should. The problem is when moving to higher dimensions, it not only gets prohibitively slow, but there are a few entries that I know (by other means) should simplify to zero but aren't, even if I replace the symbol q with it's concrete complex value.
I was recommended Cyclotomic Fields as a way to help Sage deal with fixed roots of unity more efficiently, but I can't figure out the right way (if there is one) to work over a Cyclotomic Field where the exponents can include symbols without just putting everything back into SR which defeats the purpose.
Code leaving everything in SR since that worked okay.okay. The choice to use simplify_rational() followed by simplify_expand() was just because simplify_full() was a lot slower without better results.
r = 2
var('q L')
base_field = SR
recip_diff = q^(1/2) - q^(-1/2)
def quantum(num):
return (q^(num/2) - q^(-num/2)) / 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 += (q^(n*(n-1)/4) / 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)+L)
if i > 0:
F[i,i-1] = quantum(i)
K[i,i] = q^((r-2*i-1+L)/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] = 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)
R = R.simplify_rational().expand()
# need to calculate alternate basis for V tensor V
crossE = E.tensor_product(K) + matrix.identity(r).tensor_product(E)
crossF = F.tensor_product(matrix.identity(r)) + K.inverse().tensor_product(F)
crossK = K.tensor_product(K)
crossE = crossE.simplify_rational().expand()
crossF = crossF.simplify_rational().expand()
crossK = crossK.simplify_rational().expand()
kernel = crossE.right_kernel()
eigen_vectors = crossK.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(new_basis).transpose().simplify_rational().expand()
inv_basis = new_basis.inverse().simplify_rational().expand()
newR = inv_basis * R * new_basis
