#A vector space basis of the quantum symplectic sphere
#Sophie Emma Mikkelsen
#Department of Mathematics and Computer Science,
#University of Southern Denmark, 
#Campusvej 55, 5230 Odense M, Denmark
#mikkelsen@imada.sdu.dk

#We import the library sympy to be able to work with symbols 
from sympy import *
q=symbols('q')

import time
start_time = time.time()

#Below the generators of O(S_q^{4n-1}) are collected into a set denoted X and afterwards made into a list.
#The generators are denoted xi,yi for i=1,...n
#The adjoint of xi and yi are denoted axi and ayi respectively.

#To run the program for O(S_q^{4n-1}) one have to change the number 3 below to n+1.
X=var('x(1:3) ax(1:3) y(1:3) ay(1:3)' , commutative=False)
list=list(X)

#Change n below to the one according to the algebra O(S_q^{4n-1}) you wish to investigate.
n=2
#Now the program is ready to be executed for the indicated n. 


#We divide "list" into smaller pieces which containes the generators xi, axi, yi, ayi respectively. 
xlist=list[:n]
axlist=list[n:2*n]
ylist=list[2*n:3*n]
aylist=list[3*n:4*n]

#An empty list is made which will be filled with all the commutation relations below.
relations=[]

#A commutation relation is added to the list "relations" as a list  itself, which contains two elements.
#The first is the one on the left hand side of the relation and the second the one on the right hand side. 
#The relations are computed in the order corresponding to the numbers N_1,N_2,N_(3,i), N_4,N_5 and N_6. 

########################################

#Inverted pairs containing an xi or yi and an axi or ayi. This corresponds to the number N_1.
for i in range(1,n+1):
    summ=0
    for k in range(0, i-1):
        summ=summ+axlist[k]*xlist[k]
    relations.append([xlist[i-1]*axlist[i-1],axlist[i-1]*xlist[i-1]+(1-q**(2))*summ])

summ1=0
for k in range(0, n):
    summ1=summ1+axlist[k]*xlist[k]
for i in range(1,n+1):
    summ2=0
    for k in range(i, n):
        summ2=summ2+aylist[k]*ylist[k]
    relations.append([ylist[i-1]*aylist[i-1],aylist[i-1]*ylist[i-1]+(1-q**(2))*(q**(2*(n+1-i))*axlist[i-1]*xlist[i-1]+summ1+summ2)])


for i in range(0,n):
    relations.append([xlist[i]*aylist[i],q**(2)*aylist[i]*xlist[i]])
    relations.append([ylist[i]*axlist[i],q**(2)*axlist[i]*ylist[i]])

for i in range(0,n):
    for j in range(0,n):
        if i is not j:
            relations.append([xlist[i]*axlist[j],q*axlist[j]*xlist[i]])


for i in range(0,n):
    for j in range(0,n):
        if i is not j:
            relations.append([ylist[i]*aylist[j],q*aylist[j]*ylist[i]-(q**(2)-1)*q**(2*n+2-(i+1)-(j+1))*axlist[i]*xlist[j]])
            

for i in range(0,n):
    for j in range (0,n):
        if i<j:
            relations.append([xlist[i]*aylist[j],q*aylist[j]*xlist[i]])
for i in range(0,n):
    for j in range (0,n):
        if i<j:
            relations.append([ylist[j]*axlist[i],q*axlist[i]*ylist[j]])


for i in range (0,n):
    for j in range (0,n):
        if i>j:
            relations.append([xlist[i]*aylist[j],q*aylist[j]*xlist[i]+(q**(2)-1)*q**(i-j)*aylist[i]*xlist[j]])

for i in range (0,n):
    for j in range (0,n):
        if i>j:
            relations.append([ylist[j]*axlist[i],q*axlist[i]*ylist[j]+(q**(2)-1)*q**(i-j)*axlist[j]*ylist[i]])
            

#The sphere relation. This corresponds to the number N_2.
summ3=0
for k in range(0, n-1):
    summ3=summ3+axlist[k]*xlist[k]
summ4=0
for k in range(0, n):
    summ4=summ4+aylist[k]*ylist[k]
relations.append([axlist[n-1]*xlist[n-1],1-summ3-summ4])


#Number of inverted pairs xi,yi and the number of inverted pairs axi,ayi. This corresponds to the number N_(3,i). 
for i in range(1,n+1):
    summ=0
    for k in range(0, n-(i-1)-1):
        summ=summ+(q**(n-(i-1)-1-k)*xlist[k]*ylist[k])
    relations.append([ylist[n-(i-1)-1]*xlist[n-(i-1)-1],q**(2)*xlist[n-(i-1)-1]*ylist[n-(i-1)-1]+(q**(2)-1)*summ])



for i in range(1,n+1):
    summ=0
    for k in range(0, n-(i-1)-1):
        summ=summ+(q**(n-(i-1)-1-k)*aylist[k]*axlist[k])
    relations.append([axlist[n-(i-1)-1]*aylist[n-(i-1)-1],q**(2)*aylist[n-(i-1)-1]*axlist[n-(i-1)-1]+(q**(2)-1)*summ])


#Inverted pairs involving xi,yi or axi, ayj, i not j. This corresponds to the number N_4. 
for i in range(0,n):
    for j in range(0,n):
        if i is not j:
            relations.append([ylist[j]*xlist[i],q*xlist[i]*ylist[j]])
            relations.append([axlist[j]*aylist[i],q*aylist[i]*axlist[j]])



#Inverted pairs involving yi, ayj and inverted pars involving xi, axj where i and j are different. This corresponds to the number N_5 and N_6. 

for j in range(0,n):
    for i in range(0,n):
        if 0<=i<j:
            relations.append([xlist[i]*xlist[j],q**(-1)*xlist[j]*xlist[i]])
            relations.append([ylist[i]*ylist[j],q*ylist[j]*ylist[i]])
            
        
for i in range(0,n):
    for j in range(0,n):
        if 0<=j<i:
            relations.append([axlist[i]*axlist[j],q**(-1)*axlist[j]*axlist[i]])
            relations.append([aylist[i]*aylist[j],q*aylist[j]*aylist[i]])

########################################

# reduce_algebra(expr) takes as an input a monomial. For every element [a,b] in the list "relation" it checks whether it contains an a and then replace it by b.
#This is done using "expr.subs()". 
def reduce_algebra(expr):
    old_expr = None
    while old_expr is None or not expr.equals(old_expr):
        old_expr = expr
        for x in relations:
            expr = expr.subs(x[0], x[1]).expand()
    return expr


#the_same(k,l,m) takes tree elements k,l,m in the set {xi,axi,yi,ayi| i=1,...n}. Then it reduces the monomial kl and lm.
#In the end it checks whether the two reduced expressions are the same. 
def the_same(k,l,m):
    dexpr1=k*l
    for x in relations:
        dexpr1 = dexpr1.subs(x[0], x[1]).expand()

    dexpr2=l*m
    for x in relations:
        dexpr2= dexpr2.subs(x[0], x[1]).expand()

    expr1 = expand(dexpr1*m)
    expr2 = expand(k*dexpr2)

    reduce_expr1=reduce_algebra(expr1)
    reduce_expr2=reduce_algebra(expr2)

    return reduce_expr1==reduce_expr2

#We now list the generators in the correct order in according to the basis in the Conjecture
#ie. ay1,ay2,...,ayn,ax1,ax2,...,axn,xn,...,x2,x1,yn,...,y2,y1. The list is denoted "g". 
g=[*aylist, *axlist]
for i in range(0,n):
    g.append(xlist[n-1-i])
for i in range(0,n):
    g.append(ylist[n-1-i])

#We print the list g 
print(g)

#At once we determine the ambiguites and use the function the_same() to determine if the ambiguites are resolvable.
#For every ambiguity the result True or False is included in the list "s".
#The list m includes the numbers from 0 to the number of elements in g minus 1.
#This makes it possible to subtract the elements from the list "g" depending on where they are places in "g".

m=[]
for i in range(0,len(g)):
    m.append(i)
s=[]
for i in m:
    for j in m:
        for k in m:
            if 0<j<i:
                if 0<=k<j:
                    s.append(the_same(g[i],g[j],g[k]))
                    #The following prints the ambiguity. 
                    #If the program runs for a large n this could be omitted. 
                    print (g[i],g[j],g[k])
                    

#Below all the ambiguities involving an axn*xn is determined and we check whether they are resolvable. 
z=int(len(g)/2)                   
for i in range(0,z):
    s.append(the_same(axlist[n-1],xlist[n-1],g[i]))
    #The following prints the ambiguity. 
    #If the program runs for a large n this could be omitted. 
    print(axlist[n-1],xlist[n-1],g[i])

for i in range(z, len(g)):
    s.append(the_same(g[i],axlist[n-1],xlist[n-1]))
    #The following prints the ambiguity. 
    #If the program runs for a large n this could be omitted. 
    print(g[i],axlist[n-1],xlist[n-1])

#We print the number of ambiguities. 
print(len(s))

#We now check if any of the ambiguities are not resolvable. 
if False in s:
    print(False)
else:
    print(True)

#If the program respond True in the end, then there is no False in the list "s", hence all the ambiguities are resolvable. 

#Running time: 
print("--- %s seconds ---" % (time.time() - start_time))