from sympy import *

####################################################
# Report Bugs to Sam.Nelson@cmc.edu
####################################################

#######################
# Quandles and Racks
#######################


def racktest(M):   
    ### Test whether M is a rack matrix ###
    q = True
    for i in range(1,len(M)+1):
        c = []
        for j in range(1,len(M)+1):
            for k in range(1,len(M)+1):
                if M[M[i-1][j-1]-1][k-1]!=M[M[i-1][k-1]-1][M[j-1][k-1]-1]:
                    q = False
                    break
            c = c + [M[j-1][i-1]]
        if not permtest(c):
            q = False
            break
    return q


def quandletest(M):   
    ### Test whether M is quandle matrix ###
    q = True
    for i in range(1,len(M)+1):
        if M[i-1][i-1] != i:
            q = False
            break
        c = []
        for j in range(1,len(M)+1):
            for k in range(1,len(M)+1):
                if M[M[i-1][j-1]-1][k-1]!=M[M[i-1][k-1]-1][M[j-1][k-1]-1]:
                    q = False
                    break
            c = c + [M[j-1][i-1]]
        if not permtest(c):
            q = False
            break
    return q


def subrack(L,M):
    ### find subrack of M generated by L ###
    out = []
    for x in L:
        if not x in out:
            out.append(x)
    c = True
    while c:
        c = False
        for x in out:
            for y in out:
                z = M[x-1][y-1]
                if not z in out:
                    out.append(z)
                    c = True
    return out 



def isotest(M,N):
    ### Test whether M is isomorphic to N ###
    z = []
    for i in range(1,len(M)+1):
        z = z + [0]
    L = [z]
    out = []
    while len(L) != 0 and len(out) == 0:
        w = L[0]
        L[0:1] = []  
        if w:
            i = hfindzero(w)
            if not i:
                if permtest(w): out.append(w)
            else:
                for j in pavail(w):
                    phi = list(w)
                    phi[i-1] = j
                    v = homfill(M,N,phi)
                    if v: L.append(tuple(v))
    if len(out) == 0:
        return False
    else:
        return True    
     

def reducelist(L):
    ### Remove isomorphic copies from L ###
    out = [tm(L[0])]
    W = L
    while len(W)>0:
        x = W[0]
        W[0:1] = []
        newrack = True
        for y in out:
            if isotest(x,y): newrack = False
        if newrack: out.append(tm(x))
    return out



def homlist(M,N):   
    ###Lists homomorphisms f:M->N###
    z = []
    for i in range(1,len(M)+1):
        z = z + [0]
    L = [z]
    out = []
    while len(L) != 0:
        w = L[0]
        L[0:1] = []  
        if w:
            i = hfindzero(w)
            if not i:
                out.append(w)
            else:
                for j in range(1,len(N)+1):
                    phi = list(w)
                    phi[i-1] = j
                    v = homfill(M,N,phi)
                    if v: L.append(tuple(v))
    return out     



def homfill(M,N,phi):   
    ### Fills in entries in a homomorphism ###
    f = phi
    c = True
    out = True
    while c == True:
        c = False
        for i in range(1,len(M)+1):
            for j in range(1,len(M)+1):
                if f[i-1] !=0 and f[j-1] != 0:
                    if f[M[i-1][j-1]-1] == 0 and M[i-1][j-1] != 0:
                       f[M[i-1][j-1]-1] = N[f[i-1]-1][f[j-1]-1]
                       c = True
                    elif f[M[i-1][j-1]-1] != 0 and f[M[i-1][j-1]-1] != N[f[i-1]-1][f[j-1]-1]:
                       out = False
                       c = False
    if out == True: 
       return f
    else:
       return out


def hfindzero(f):   
    ### find zero in homomorphism template ###
    j = -1
    for i in range(0,len(f)):
        if f[i] == 0: 
            j = i+1
            break
    if j < 0: out = False
    else: out = j
    return out


def rackrank(M):
    ### find rack rank ###
    rr = []
    for i in range(1,len(M)+1):
        c = 1
        j = M[i-1][i-1]
        while j != i:
            c = c+1
            j = M[j-1][j-1]
        rr.append(c)
    out = 1
    for x in rr:
        out = lcm(out,x)
    return out


########################
# Generating racks
########################


def reptest(p):   
    ### Test whether p has repeated non-zero entries ###
    q = True
    L = []
    for i in range(1,len(p)+1):
        if p[i-1] != 0:
            if p[i-1] in L:
                q = False
            else:
                L.append(p[i-1])
    return q



def sdfill(N):
    ### fill with self-distributiviity ###
    c = True
    con = False
    M = lm(N)
    while c:
        c = False
        for i in range(1,len(M)+1):
            for j in range(1,len(M)+1):
                for k in range(1,len(M)+1):
                    if M[i-1][j-1] != 0 and M[M[i-1][j-1]-1][k-1] != 0:
                        if M[i-1][k-1] != 0 and M[j-1][k-1] != 0:
                            if M[M[i-1][j-1]-1][k-1] != M[M[i-1][k-1]-1][M[j-1][k-1]-1]:
                                if M[M[i-1][j-1]-1][k-1] == 0:
                                    M[M[i-1][j-1]-1][k-1] = M[M[i-1][k-1]-1][M[j-1][k-1]-1]
                                    c = True
                                elif M[M[i-1][k-1]-1][M[j-1][k-1]-1] == 0:
                                    M[M[i-1][k-1]-1][M[j-1][k-1]-1] = M[M[i-1][j-1]-1][k-1]
                                    c = True
                                else:
                                    con = True
    if con: M = False
    return M
     


def findzero(M):
    ### Find position of first zero entry in a matrix ###
    out = False
    for i in range(1,len(M)+1):
        if not out:
            for j in range(1,len(M[0])+1):
                if M[i-1][j-1] == 0:
                    out = (i,j)
                    break
    return out



def pavail(v):
    ### List available entries ###
    L = []
    for i in range(1,len(v)+1):
        if not i in v: L.append(i)
    return tuple(L)



def racklist(N):
    ### find all racks matching pattern ###
    L = []
    L.append(tm(N))
    out = []
    while len(L)>0:
       M = L[0]
       L[0:1] = []
       q = findzero(M)
       if q:
           for i in pavail(getcolumn(M,q[1])):
               M2 = lm(M)
               M2[q[0]-1][q[1]-1] = i
               M3 = sdfill(M2)
               if M3:
                    col = True 
                    for j in range(1,len(M)+1):
                        if not reptest(getcolumn(M3,j)):
                            col = False
                    if col:  L.append(tm(M3))
       else:
           out.append(tuple(M))
    return out


def rackaut(M):
    out = []
    for x in homlist(M,M):
        if permtest(x): out.append(x)
    return out


################################
# Mod n linear algebra stuff
###############################

def getcolumn(M,j):   # get column n from matrix M
    ### Gets column n from matrix M### 
    c =[]
    for i in range(1,len(M)+1):
        c.append(M[i-1][j-1])
    return c



def tm(M):
    ### tuple matrix### 
    out  = []
    for x in M: out.append(tuple(x))
    return tuple(out)



def lm(M):
    ### list matrix### 
    out = []
    for x in M: out.append(list(x))
    return list(out)



def mtranspose(M):
    ### transpose matrix### 
    out=[]
    for i in range(1,len(M[0])+1):
        r = []
        for j in range(1,len(M)+1):
            r.append(M[j-1][i-1])
        out.append(r)
    return out



def matrixdisplay(M):
    ### readable matrix display### 
    m = 1
    for i in range(0,len(M)):
        for j in range(0,len(M[i])):
            if len(str(M[i][j]))>m: m = len(str(M[i][j]))
    print '\n'
    for i in range(0,len(M)):
        for j in range(0,len(M[i])):
            print repr(M[i][j]).center(m+1),
        print '\n'



def mm(M,N,n):
    ### multiply matrices in Zn (n==0 for Z or Q)### 
    if len(M[0]) != len(N): 
        return False
    else:
       out = []
       for i in range(0,len(M)):
           r = []
           for j in range(0,len(N[0])):
               t = 0
               for k in range(0,len(M[0])):
                   if n != 0 and n != 1: 
                       t = (t + M[i][k]*N[k][j]) % n
                   else:
                       t = (t + M[i][k]*N[k][j]) 
               r.append(t)
           out.append(r)
    return out



def vp(v,w,n):
    ### add two vectors mod n### 
    if len(v) != len(w): 
        return False
    else:
        out = []
        for j in range(0,len(v)):
            if n != 0:
                out.append((v[j]+w[j]) % n)        
            else:
                out.append(v[j]+w[j])
    return out


def sm(a,v,n):
    ### scalar multiply a vector mod n### 
    out = []
    for j in range(0,len(v)):
        if n == 0:
            out.append(a*v[j])        
        else:
            out.append((a*v[j]) % n)
    return out



def invn(k,n):
    ### k inverse mod n### 
    out = False
    if n == 0:
        return Rational(1)/k
    for j in range(1,n):
        if j*k %n == 1: out = j
    return out



def lexinc(v,n):
    ### increment in dictionary order mod n### 
    p = len(v)-1   
    while p != -1 and v[p] == n-1:
        p = p-1
    if p == -1: 
        return False
    else:
        w = list(v)
        w[p] = w[p]+1 % n
        for j in range(p+1,len(w)):
            w[j] = 0
    return w



def znm(n,m):
    ### generate (Zn)^m### 
    if m==1:
        out=[]
        for i in range(0,n):
            out.append([i])
        return out
    f = []
    for j in range(0,m):
        f.append(0)
    L = []
    while f:
       L.append(tuple(f))
       f = lexinc(f,n)
    return L
  


def matnm(n,m):
    ### generate M_m(n)### 
    out = []
    L=znm(n,m*m)
    for x in L:
        t=[]
        for j in range(0,m):
            t.append(x[j*m:j*m+m])
        out.append(tm(t))
    return out
  


def znmopmat(n,m):
    ###  operation matrix for (Zn)^m### 
    M = []
    L = znm(n,m)
    for x in range(1,len(L)+1):
        r = []
        for y in range(1,len(L)+1):
            z = vp(L[x-1],L[y-1],n)
            for k in range(1,len(L)+1):
                if tuple(z) == L[k-1]: r.append(k)
        M.append(r)
    return M



def auttest(M,f):
    ###  test whether f is an automorphism of M### 
    for i in range(1,len(M)+1):
        for j in range(1,len(M)+1):
            if f[M[i-1][j-1]-1] != M[f[i-1]-1][f[j-1]-1]:
                return False
    return True



def autconjtest(M,f,g):
    ###  test whether automorphisms f,g of M are conjugate### 
    A=rackaut(M)
    for a in A:
       if permcomp(f,a) == permcomp(g,a):
           return True
    return False



def conjclasslist(M):
    Au = rackaut(M)
    out = [Au[0]]
    for b in Au:
        z = True
        for x in out:
            for y in Au:
                if permcomp(b,y) == permcomp(y,x):
                    z = False
        if z: out.append(b)
    return out



def submod(L,n):
    ### find submodule of (Zn)^m generated by L### 
    out = []
    for x in L:
        if not x in out:
            out.append(x)
    c = True
    while c:
        c = False
        for x in out:
            for a in range(0,n-1):
                z = tuple(sm(a,x,n))
                if not z in out:
                    out.append(z)
                    c = True
            for y in out:
                z = tuple(vp(x,y,n))
                if not z in out:
                    out.append(z)
                    c = True
    return out



def normalize(M,n):
    ### Divide each row by leading entry ### 
    out=lm(M)
    for i in range(0,len(M)):
        k=lead(out[i])
        if k!= len(M):
            x=invn(out[i][k],n)
            if x: 
                out[i]=sm(x,out[i],n)
    return out



def lead(x):
    ### Position of leading entry in x### 
    i=0
    while i<len(x) and x[i] == 0:
        i = i+1
    return i



def leadsort(M):
    ### Order by leading entry position### 
    N = M
    for j in range(0,len(M)):
        for i in range(0,j):
            if lead(N[i])>lead(N[j]):
                N=rowswitch(N,i,j)
    return N




def gaussjord(M,n):
    ### Gauss-Jordan elimination ### 
    out=leadsort(normalize(lm(M),n))
    for i in range(0,len(M)):
        if lead(out[i])!=len(M):
            for j in range(lead(out[i])+1,len(out)):
                out[j]= vp(out[j],sm(-1*out[j][lead(out[i])],out[i],n),n)
            out=leadsort(normalize(out,n))
    for i in range(0,len(M)):
        if lead(out[i])!=len(M):
            for j in range(0,lead(out[i])):
                out[j]= vp(out[j],sm(-1*out[j][lead(out[i])],out[i],n),n)
            out=leadsort(normalize(out,n))
    return out




def rowswitch(M,i,j):
    ### switch rows i and j### 
    out = []
    for k in range(0,len(M)):
        if k==i: out.append(M[j])
        if k==j: out.append(M[i])
        if k!= i and k!= j:out.append(M[k])
    return out


def rowreplace(M,i,v):
    ### replace row i with vector v### 
    out = []
    for k in range(0,i-1):
        out.append(M[k])
    out.append(v)
    for k in range(i,len(M)):
        out.append(M[k])
    return out



def getkernel(M1,n):
    ###  get a basis for the kernel of M1 mod n### 
    M = gaussjord(M1,n)
    L = leadingentries(M)
    L1 = []
    out = []
    for x in L:
        if x[0] != len(M[0]):
            L1.append(x[0])
    for i in range(0,len(M[0])):
        if not i in L1:
            v =[]
            for j in range(0,i):
                if j in L1:
                    for k in range(0,len(L1)):
                        if j==L[k][0]:
                            v.append(M[L[k][1]][i])
                else:
                    v.append(0)
            v.append(n-1)
            for j in range(i+1,len(M[0])):
                v.append(0)
            out.append(tuple(v))
    return out



def zeromatrix(n):
    ### return n by n zero matrix### 
    out = []
    for i in range(0,n): 
        r = []
        for j in range(0,n): r.append(0)
        out.append(r)
    return out



def nmzeromatrix(n,m):
    ### return n by m zero matrix### 
    out = []
    for i in range(0,n): 
        r = []
        for j in range(0,m): r.append(0)
        out.append(r)
    return out



def minor(M,i,j):
    ###  i,j minor of M ### 
    out=[]
    for k in range(0,len(M)):
        if k != i:
            r=[]
            for l in range(0,len(M[0])):
                if l != j: r.append(M[k][l])
            out.append(r)
    return out



def detn(M,n):
    ### determinant of M mod n ### 
    m,m1=len(M),len(M[0])
    if m!=m1: return False
    if m == 1: return ((M[0][0]) % n)
    out = 0
    for k in range(0,m):
        out = (out+M[0][k]*(-1)**(k % 2)*detn(minor(M,0,k),n)) %n
    return out



def cofacn(M,n):
    ### Cofactor matrix mod n ### 
    out=[]
    for i in range(0,len(M)):
        r = []
        for j in range(0,len(M[0])):
            r.append((-1)**(i+j)*detn(minor(M,i,j),n) % n)
        out.append(r)
    return out



def adjn(M,n):
    ### Adjoint ### 
    return mtranspose(cofacn(M,n))



def rowint(X,Y,i):
    ###  ### 
    x,y=X,Y
    if x[i]<0: x=sm(-1,x,0)
    if y[i]<0: y=sm(-1,y,0)
    while (x[i]!=0 and y[i]!= 0):
        if x[i] <= y[i]:
            y=vp(y,sm(-1,x,0),0)
        else:
            x=vp(x,sm(-1,y,0),0)
    if x[i]!=0:
        return (x,y)
    else:
        return (y,x)



###########################
# Rack homology
###########################


def reducedrackcocycles(R,n):
    ###Find a basis for the mod n 2-cocycles of rack R###
    M=[]
    m=len(R)
    r=[]
    N=rackrank(R)
    for i in range(0,m**2):
        r.append(0)
    for i in range(0,m**3+m):
        M.append(tuple(r))
    M=lm(M)
    j = 0
    for a in range(1,m+1):
        y = a
        for i in range(0,N):
            M[j][m*(y-1)+y-1]=(M[j][m*(y-1)+y-1]+1) %n
            y = R[y-1][y-1]
        j=j+1
        for b in range(1,m+1):
            for c in range(1,m+1):
                M[j][m*(a-1)+b-1]=(M[j][m*(a-1)+b-1]+1)%n
                M[j][m*(R[a-1][b-1]-1)+c-1]=(M[j][m*(R[a-1][b-1]-1)+c-1]+1)%n
                M[j][m*(a-1)+c-1]=(M[j][m*(a-1)+c-1]-1) %n
                M[j][m*(R[a-1][c-1]-1)+R[b-1][c-1]-1]=(M[j][m*(R[a-1][c-1]-1)+R[b-1][c-1]-1]-1) %n
                j=j+1
    out = getkernel(M,n)
    return out



def rcocycles(R,n):
    ### reduced rack cocycles over Z_n ###
    m,N=len(R),rackrank(R)
    bl=[]
    for k in range(0,m*m):
        bl.append(Z)
    working,out=[bl],[]
    while len(working)>0:
        f=working[0]
        working[0:1]=[]
        k=hfindz(f)
        if k == -1:
            if reducedcocycletest(R,f,n):
                out.append(tuple(f))
        else:
            for j in range(0,n):
                f1=list(f)
                f1[k]=j
                f3=cocyclefill(f1,R,n)
                if f3:
                    working.append(tuple(f3))
    return out





from random import randrange

def randcocycle(R,n):
    ### reduced rack cocycles over Z_n ###
    m,N=len(R),rackrank(R)
    bl=[]
    for k in range(0,m*m):
        bl.append(Z)
    working=[bl]
    while len(working)>0:
        f=working[0]
        working[0:1]=[]
        k=hfindz(f)
        if k == -1:
            if reducedcocycletest(R,f,n):
                return f
        else:
            j = randrange(0,n)
            f1=list(f)
            f1[k]=j
            f3=cocyclefill(f1,R,n)
            if f3:
                working.append(tuple(f3))
            else: 
                working.append(bl)


def cocyclefill(f,R,n):
    ### fill rack cocycle ###
    v,m=list(f),len(R)
    c=True
    while c:
        c=False
        for x in range(1,m+1):
            for y in range(1,m+1):
                for z in range(1,m+1):
                    if v[m*(x-1)+y-1] ==Z and v[m*(R[x-1][y-1]-1)+z-1]!=Z and v[m*(x-1)+z-1]!=Z and v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1]!=Z:
                        v[m*(x-1)+y-1] = (v[m*(x-1)+z-1] + v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1] - v[m*(R[x-1][y-1]-1)+z-1]) %n
                        c=True
                    if v[m*(x-1)+y-1] !=Z and v[m*(R[x-1][y-1]-1)+z-1]==Z and v[m*(x-1)+z-1]!=Z and v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1]!=Z:
                        v[m*(R[x-1][y-1]-1)+z-1] = (v[m*(x-1)+z-1] + v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1] - v[m*(x-1)+y-1]) %n
                        c=True
                    if v[m*(x-1)+y-1] !=Z and v[m*(R[x-1][y-1]-1)+z-1]!=Z and v[m*(x-1)+z-1]==Z and v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1]!=Z:
                        v[m*(x-1)+z-1] = ( v[m*(x-1)+y-1] + v[m*(R[x-1][y-1]-1)+z-1] - v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1] ) %n
                        c=True
                    if v[m*(x-1)+y-1] !=Z and v[m*(R[x-1][y-1]-1)+z-1]!=Z and v[m*(x-1)+z-1]!=Z and v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1]==Z:
                        v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1] = ( v[m*(x-1)+y-1] + v[m*(R[x-1][y-1]-1)+z-1] - v[m*(x-1)+z-1] ) %n
                        c=True
                    if v[m*(x-1)+y-1] !=Z and v[m*(R[x-1][y-1]-1)+z-1]!=Z and v[m*(x-1)+z-1]!=Z and v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1]!=Z:
                        if (v[m*(x-1)+y-1] + v[m*(R[x-1][y-1]-1)+z-1])%n != (v[m*(x-1)+z-1]+ v[m*(R[x-1][z-1]-1)+R[y-1][z-1]-1])%n:
                            return False
    return tuple(v)




def reducedcocycletest(R,v,n):
    ### Test whether v is a reduced 2-cocyle mod n ###
    ### use n==0 for integer case ###
    N=rackrank(R)
    m=len(R)
    out = True
    for x in range(1,len(R)+1):
        for y in range(1,len(R)+1):
            for z in range(1,len(R)+1): 
   



def grackhominv1(G,T,v,n):
    ###Gauss code rack cocycle invariant###
    N = rackrank(T)
    m = len(T)
    u = Symbol('u')
    ff = 0
    w = G
    for i in range(0,N):
        D=gauss2pd(w)
        H = pdhomlist(D,T)
        for h in H:
            bw=0
            for x in D:
                if x[0]==1:
                    bw = (bw + v[m*(h[x[1]-1]-1)+h[x[4]-1]-1]) %n
                if x[0] == -1:
                     bw = (bw - v[m*(h[x[3]-1]-1)+h[x[4]-1]-1]) %n
            ff = ff + z**(bw % n)
        w = gaussinc(w,1)
    return ff



def grackhominv2(G,T,v,n):
    ###Gauss code rack cocycle invariant###
    N = rackrank(T)
    m = len(T)
    u = Symbol('u')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            D=gauss2pd(w)
            H = pdhomlist(D,T)
            for h in H:
                bw=0
                for x in D:
                    if x[0]==1:
                        bw = (bw + v[m*(h[x[1]-1]-1)+h[x[4]-1]-1]) %n
                    if x[0] == -1:
                        bw = (bw - v[m*(h[x[3]-1]-1)+h[x[4]-1]-1]) %n
                ff = ff + z**(bw % n)
            w = gaussinc(w,1)
        w = gaussinc(w,2)
    return ff


def grackhominv3(G,T,v,n):
    ###Gauss code rack cocycle invariant###
    N = rackrank(T)
    m = len(T)
    u = Symbol('u')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for l in range(0,N):
                D=gauss2pd(w)
                H = pdhomlist(D,T)
                for h in H:
                    bw=0
                    for x in D:
                        if x[0]==1:
                            bw = (bw + v[m*(h[x[1]-1]-1)+h[x[4]-1]-1]) %n
                        if x[0] == -1:
                            bw = (bw - v[m*(h[x[3]-1]-1)+h[x[4]-1]-1]) %n
                    ff = ff + z**(bw % n)
                w = gaussinc(w,1)
            w = gaussinc(w,2)
        w = gaussinc(w,3)
    return ff


def grackhominv(G,T,v,n):
    ###Gauss code rack cocycle invariant###
    if len(G) == 1:
        return grackhominv1(G,T,v,n)
    if len(G) == 2:
        return grackhominv2(G,T,v,n)
    if len(G) == 3:
        return grackhominv3(G,T,v,n)


def grackhominvset(G,T,v,n):
    ###Gauss code rack cocycle invariant###
    N = rackrank(T)
    m = len(T)
    z = Symbol('z')
    w = G
    out = []
    for i in range(0,N):
        out.append(0)
    for i in range(0,N):
        D = gauss2pd(w)
        H = pdhomlist(D,T)
        for h in H:
            bw=0
            for x in D:
                if x[0] ==  1:
                    bw = (bw + v[ m*(h[x[1]-1]-1) +h[x[2]-1] -1 ]) %n
                if x[0] == -1:
                    bw = (bw - v[ m*(h[x[3]-1]-1) +h[x[2]-1] -1 ]) %n
            out[gausswrithe(w[0]) %N] = out[gausswrithe(w[0]) %N] + z**bw
        w = gaussinc(w,1)
    return out



##########################################
# Permuation group/Number theory stuff
##########################################

def permtest(p):   # test whether p is a permutation
    ### Test whether p is a permutation### 
    q = True
    for i in range(1,len(p)+1):
        if not i in p:
            q = False
    return q


def permcomp(p,q):
    ### Compose permutations### 
    out = []
    for i in range(1,len(p)+1):
        out.append(q[p[i-1]-1])
    return tuple(out)


def subgroup(L):
    ### complete subgroup of permutation group### 
    out = []
    for x in L:
        if not tuple(x) in out: out.append(tuple(x))
    c = True
    while c:
        c = False
        for x in out:
            for y in out:
                z = tuple(permcomp(x,y))
                if not z in out:
                    out.append(tuple(z))
                    c = True
    return out




def permorder(p):
    ### order of a permutation p### 
    i = 1
    x = list(p)
    id = []
    for j in range(1,len(p)+1):
        id.append(j)
    while x != id:
        x = list(permcomp(x,p))
        i = i + 1
    return i 


def gcd(x,y):
    ### Greatest common divisor### 
    if y == 0:
        return x
    return gcd(y, x % y)


def lcm(x,y):
    ### Least common multiple### 
    return(x*y/gcd(x,y))


def groupexponent(L):
    ### Compute exponent of group### 
    S = subgroup(L)
    out = 1
    for x in S:
        out = lcm(out,permorder(x))
    return out


def conjquandle(L):
   ### Conjugation quandle of group generated by L### 
   S=subgroup(L)
   M=zeromatrix(len(S))
   for i in range(1,len(S)+1):
       for j in range(1,len(S)+1):
           x=permcomp(permcomp(perminv(S[j-1]),S[i-1]),S[j-1])
           for k in range(1,len(S)+1):
               if S[k-1] == x: M[i-1][j-1] = k
           
   return M

def nconjquandle(L,n):
   ### n-fold Conjugation quandle of n-conj class generated by L### 
   S=[]
   for x in L: 
       if tuple(x) not in S: 
           S.append(tuple(x))
   c=True
   while c:
       c=False
       for x in S:
           for y in S:
               #z = permcomp(permcomp(perminv(x),y),x)
               z=nconj(x,y,n)
               if not tuple(z) in S: 
                   S.append(tuple(z))
                   c=True
   M=zeromatrix(len(S))
   for i in range(1,len(S)+1):
       for j in range(1,len(S)+1):
           #x=permcomp(permcomp(perminv(S[j-1]),S[i-1]),S[j-1])
           x=nconj(S[i-1],S[j-1],n)
           for k in range(1,len(S)+1):
               if S[k-1] == x: M[i-1][j-1] = k
           
   return M

def nconj(x,y,n):
    ###  y^-n x y^n ### 
    z=x
    for i in range(0,n):
        z = permcomp(permcomp(perminv(y),z),y)
    return z



###########################
# Planar Diagram stuff
###########################


def pdhomlist(PD,N):   # list quandle homomorphisms
    ###Lists homomorphisms from knot rack/quandle of planar diagram###
    z = []
    for i in range(1,2*len(PD)+1):
        z = z + [0]
    L = [z]
    out = []
    while len(L) != 0:
        w = L[0]
        L[0:1] = []  
        if w:
            i = hfindzero(w)
            if not i:
                out.append(w)
            else:
                for j in range(1,len(N)+1):
                    phi = list(w)
                    phi[i-1] = j
                    v = pdhomfill(PD,N,phi)
                    if v: L.append(tuple(v))
    return out     


def pdhomfill(PD,N,phi):   # fill in homomorphism
    ###Fills in entries in a homomorphism###
    f = phi
    c = True
    out = True
    while c == True:
        c = False
        for X in PD:
            if f[X[2]-1] == 0 and f[X[4]-1] != 0:
                f[X[2]-1] = f[X[4]-1]
                c = True
            if f[X[4]-1] == 0 and f[X[2]-1] != 0:
                f[X[4]-1] = f[X[2]-1]
                c = True
            if f[X[4]-1] != f[X[2]-1] and (f[X[2]-1] != 0 and f[X[4]-1] != 0):
                out = False
            if X[0] == 1:
                if f[X[1]-1] != 0 and f[X[4]-1] != 0:
                    if f[X[3]-1] == 0: 
                        f[X[3]-1] = N[f[X[1]-1]-1][f[X[4]-1]-1]
                        c = True
                    elif f[X[3]-1] != N[f[X[1]-1]-1][f[X[4]-1]-1]:
                        out = False
            if X[0] == -1:
                if f[X[3]-1] != 0 and f[X[4]-1] != 0:
                    if f[X[1]-1] == 0: 
                        f[X[1]-1] = N[f[X[3]-1]-1][f[X[4]-1]-1]
                        c = True
                    elif f[X[1]-1] != N[f[X[3]-1]-1][f[X[4]-1]-1]:
                        out = False
    if out == True: 
       return f
    else:
       return out


def quandlecount(PD,M):
    ###quandle counting invariant###
    return len(pdhomlist(PD,M))


def rackcount(PD,M):
    ###rack counting invariant list version###
    N = rackrank(M)
    out = range(0,N)
    w = PD
    for n in range(0,N):
        out[writhe(w) % N] = quandlecount(w,M)
        w = pdwritheinc(w)
    return out


def rackcountpoly(PD,M):
    ###rack counting invariant polynomial version###
    N = rackrank(M)
    q = Symbol('q')
    x = Symbol('x')
    x = 0
    w = PD
    for n in range(0,N):
        x = x + quandlecount(w,M)*q**(writhe(w) % N)
        w = pdwritheinc(w)
    return x


def grackcountpoly(g,M):
    ###rack counting invariant polynomial version###
    N = rackrank(M)
    q = Symbol('q')
    x = Symbol('x')
    x = 0
    w = g
    for n in range(0,N):
        x = x + quandlecount(gauss2pd(w),M)*q**(gausswrithe(w[0])% N)
        w = gaussinc(w,1)
    return x



def writhe(PD):
    ###writhe of a planar diagram###
    w = 0
    for X in PD:
        if X[0] == 1: w = w+1
        if X[0] == -1: w = w-1
    return w


def pdwritheinc(PD):
    ###increment writhe of planar diagram ###
    out = []
    inced = False
    for X in PD:
       if X[2] == 2*len(PD) and X[4] == 1 and inced == False:
           out.append((X[0],X[1],X[2]+2,X[3],X[4]))
           out.append((1,X[2]+1,X[2]+1,X[2]+2,X[2]))
           inced = True
       elif X[1] == 2*len(PD) and X[3] ==1 and inced == False:
           out.append((X[0],X[1]+2,X[2],X[3],X[4]))
           out.append((1,X[1]+1,X[1]+1,X[1]+2,X[1])) 
           inced = True
       elif X[4] == 2*len(PD) and X[2]==1 and inced == False:
           out.append((X[0],X[1],X[2],X[3],X[4]+2))
           out.append((1,X[4]+1,X[4]+1,X[4]+2,X[4]))
           inced = True
       else: out.append(X)
    return out



################################
# Knot table -- planar diagrams
################################

def knot(n):
    ### planar diagam ###
    return gauss2pd(gknot[n])


##############################
# Gauss code table
##############################

gknot = {
  (0):[[-1,1]],
(3,1):[[-1,2,-3,1,-2,3]],
(4,1):[[-1.5,2.5,-3,4,-2.5,1.5,-4,3]],
(5,1):[[-1,2,-3,4,-5,1,-2,3,-4,5]],
(5,2):[[-1,2,-3,4,-5,3,-2,1,-4,5]],
(6,1):[[-1,2,-3.5,4.5,-5.5,6.5,-2,1,-6.5,5.5,-4.5,3.5]],
(6,2):[[-1,2.5,-3.5,4.5,-5.5,1,-6,3.5,-4.5,5.5,-2.5,6]],
(6,3):[[-1,2,-3.5,4.5,-5.5,1,-6.5,3.5,-4.5,6.5,-2,5.5]],
(6,4):[[-1,2,-3,1,-2,3,-4,5,-6,4,-5,6]],
(6,5):[[-1,2,-3,1,-2,3,4.5,-5.5,6.5,-4.5,5.5,-6.5]],
(7,1):[[-1,2,-3,4,-5,6,-7,1,-2,3,-4,5,-6,7]],
(7,2):[[-1.5,6.5,-7.5,1.5,-2.5,3.5,-4.5,5.5,-6.5,7.5,-5.5,4.5,-3.5,2.5]],
(7,3):[[-1,2,-3,4,-7,6,-5,1,-2,3,-4,5,-6,7]],
(7,4):[[-1,2,-3,4,-5,6,-7,1,-4,3,-2,7,-6,5]],
(7,5):[[-1.5,2.5,-3.5,4.5,-5.5,1.5,-2.5,3.5,-6.5,7.5,-4.5,5.5,-7.5,6.5]],
(7,6):[[-1.5,2.5,-3,4,-5.5,6.5,-4,3,-7.5,1.5,-6.5,5.5,-2.5,7.5]],
(7,7):[[-1.5,2.5,-3,4,-2.5,5.5,-6,7,-5.5,1.5,-4,3,-7,6]],
(8,1):[[-1.5,2.5,-3,4,-2.5,1.5,-5.5,6.5,-7.5,8.5,-4,3,-8.5,7.5,-6.5,5.5]],
(8,2):[[1,-2,3.5,-4.5,5.5,-6.5,7.5,-8.5,2,-1,8.5,-3.5,4.5,-5.5,6.5,-7.5]],
(8,3):[[-1,2,-3,4,-5.5,6.5,-7.5,8.5,-4,3,-2,1,-8.5,7.5,-6.5,5.5]],
(8,4):[[-1.5,2,-3,4,-5,1.5,-6.5,7.5,-8.5,5,-4,3,-2,6.5,-7.5,8.5]],
(8,5):[[-1.5,2.5,-3,4,-5,6,-7,8,-2.5,1.5,-6,7,-8,3,-4,5]],
(8,6):[[1,-2,3.5,-4.5,5.5,-6.5,7.5,-8.5,2,-1,8.5,-7.5,6.5,-3.5,4.5,-5.5]],
(8,7):[[-1.5,2.5,-3.5,1.5,-4,5,-6,7,-8,4,-2.5,3.5,-5,6,-7,8]],
(8,8):[[1,-2,3,-4.5,5.5,-6.5,4.5,-7,8,-5.5,6.5,-3,2,-1,7,-8]],
(8,9):[[-1.5,2.5,-3.5,4.5,-5,6,-7,8,-2.5,3.5,-4.5,1.5,-8,5,-6,7]],
(8,10):[[-1,2,-3,4,-5,6,-7,1,-2,8.5,-6,7,-8.5,3,-4,5]],
(8,11):[[1,-2,3.5,-4.5,5.5,-6.5,7.5,-8.5,2,-1,8.5,-5.5,4.5,-3.5,6.5,-7.5]],
(8,12):[[-1,2,-3,4,-5.5,6.5,-4,3,-7.5,8.5,-2,1,-8.5,7.5,-6.5,5.5]],
(8,13):[[1.5,-2.5,3,-4,5,-6,7,-3,8.5,-1.5,2.5,-8.5,4,-7,6,-5]],
(8,14):[[-1,2,-3.5,4.5,-5.5,6.5,-7.5,5.5,-4.5,8.5,-2,1,-6.5,7.5,-8.5,3.5]],
(8,15):[[1.5,-2.5,3.5,-4.5,5.5,-3.5,6.5,-7.5,8.5,-6.5,2.5,-1.5,7.5,-8.5,4.5,-5.5]],
(8,16):[[1.5,-2,3,-4.5,5.5,-6,2,-7.5,4.5,-5.5,8.5,-1.5,7.5,-3,6,-8.5]],
(8,17):[[1.5,-2.5,3,-4,5.5,-1.5,6.5,-3,4,-7.5,8.5,-5.5,2.5,-6.5,7.5,-8.5]],
(8,18):[[-1,2.5,-3.5,4,-5,6.5,-2.5,7,-4,8.5,-6.5,1,-7,3.5,-8.5,5]],
(8,19):[[-1,2,-3,-4,5,1,-2,-6,4,7,-8,-5,6,3,-7,8]],
(8,20):[[-1,2.5,3,-4,-5.5,1,6.5,-3,4,-7.5,8.5,5.5,-2.5,-6.5,7.5,-8.5]],
(8,21):[[-1,2.5,-3.5,-4.5,5.5,-6,-7.5,1,4.5,-5.5,8.5,7.5,-2.5,3.5,6,-8.5]],
(9,2):[[-1.5,2.5,-3.5,4.5,-5.5,6.5,-7.5,8.5,-9.5,1.5,-2.5,9.5,-8.5,7.5,-6.5,5.5,-4.5,3.5]],
(9,24):[[-1.5,2.5,-3,4.5,-5.5,6,-7,8,-2.5,1.5,-6,7,-8,9.5,-4.5,5.5,-9.5,3]]
}

glink = {
(0,2):[[-1,1],[-2,2]],
(0,3):[[-1,1],[-2,2],[-3,3]],
(2,0,1):[[-1,2],[-2,1]],
(4,0,1):[[-1,2,-3,4],[-4,3,-2,1]],
(5,0,1):[[-1,2.5,-3.5,4],[1,-4,5.5,-2.5,3.5,-5.5]],
(6,0,1):[[-1,2,-3,4],[1,-5.5,6.5,-4,3,-6.5,5.5,-2]],
(6,0,2):[[-1,2,-5.5,4,-3.5,6],[1,-2,3.5,-4,5.5,-6]],
(6,0,3):[[-1,2,-3,4,-5,6],[1,-2,3,-4,5,-6]],
(6,0,4):[[-1,2,-3,4],[-5,1,-6,3],[-2,6,-4,5]],
(6,0,5):[[-1,2.5,-3.5,4],[1,-5,6,-2.5],[3.5,-6,5,-4]],
(6,1,1):[[-1,-2,3,4],[1,-5,-3,6],[2,5,-4,-6]],
(7,0,1):[[-1.5,2.5,-3,4,-5.5,1.5,-6,3,-7.5,5.5],[-2.5,7.5,-4,6]],
(7,0,2):[[-1,2,-3,4,-5,1,-2,5,-6,7],[3,-7,6,-4]],
(7,0,3):[[-1.5,2.5,-3.5,4.5,-5.5,1.5,-2.5,3.5,-6,7],[-4.5,5.5,-7,6]],
(7,0,4):[[-1.5,2.5,-3.5,4.5,-5.5,3.5,-2.5,1.5,-6,7],[-4.5,5.5,-7,6]],
(7,0,5):[[-1,2.5,-3.5,1,-4,5,-6,7],[-2.5,4,-7,6,-5,3.5]],
(7,0,6):[[-1,2.5,-3.5,4.5,-2.5,5,-6,7],[1,-7,6,-5,3.5,-4.5]],
(7,0,7):[[-1,2,-3,4],[1,-5.5,6.5,-2],[3,-7.5,5.5,-6.5,7.5,-4]],
(7,1,1):[[-1,2.5,3,-4.5,-5,1,6.5,-3,-7.5,5],[-2.5,-6.5,4.5,7.5]],
(7,1,2):[[-1,2.5,3,-4.5,-5,1,-6,-3,7,5],[-2.5,6,4.5,-7]]
}



#############################
# Coxter rack stuff
#############################


def cr(v,w,a,M,n):
    ###coxeter rack operation###
    if (len(M) == len(v)) and (len(v) == len(w)):
        return sm((-a) %n,vp(sm(2*mm(mm([list(w)],M,n),mtranspose([list(v)]),n)[0][0]*invn(mm(mm([list(w)],M,n),mtranspose([list(w)]),n)[0][0],n),list(w),n),sm(-1,list(v),n),n),n)
    else:
        return False


def coxrack(n,m,a,M):
    ###Coxeter rack on (Zn)^m with scalar a and bilinear form <x,y> = xMy^t###
    L = znm(n,m)
    C = []
    for x in L:
        if gcd(mm(mm([list(x)],M,n),mtranspose([list(x)]),n)[0][0],n) == 1:
            C.append(tuple(x))
    out = []
    for i in C:
        r = []
        for j in C:
            z = cr(list(i),list(j),a,M,n)
            for k in range(1,len(C)+1):
                if tuple(z) == C[k-1]:
                    r.append(k)
        out.append(r)
    return out


def coxrack2(n,m,a,M):
    ###Coxeter rack on (Zn)^m with scalar a and bilinear form <x,y> = xMy^t###
    L = znm(n,m)
    C = []
    for x in L:
        if gcd(mm(mm([list(x)],M,n),mtranspose([list(x)]),n)[0][0],n) == 1:
            C.append(tuple(x))
    out = []
    for i in C:
        r = []
        for j in C:
            z = cr(list(i),list(j),a,M,n)
            for k in range(1,len(C)+1):
                if tuple(z) == C[k-1]:
                    r.append(k)
        out.append(r)
    return (out,C)



def coxinv(PD,n,m,a,M):
    ###compute Coxeter rack invariant###
    L = znm(n,m)
    C = []
    for x in L:
        if gcd(mm(mm([list(x)],M,n),mtranspose([list(x)]),n)[0][0],n) == 1:
            C.append(tuple(x))
    T = coxrack(n,m,a,M)
    N = rackrank(T)
    q = Symbol('q')
    t = Symbol('t')
    s = Symbol('s')
    ff = Symbol('ff')
    ff = 0
    w = PD
    for zz in range(0,N):
        H = pdhomlist(w,T)
        for h in H:
            h2 = subrack(h,T)
            SM = []
            for j in h2:        
                SM.append(C[j-1])
            ff = ff + q**(writhe(w) % N)*t**(len(h2))*s**(len(submod(SM,n)))
        w = pdwritheinc(w)
    return ff    



###########################
# symplectic quandles
###########################



def symplecticquandle(n,m,M):
    ###Symplectic quandle on (Zn)^m with bilinear form <x,y> = xMy^t###
    L = znm(n,m)
    out = []
    for i in range(1,n**m+1):
        r = []
        for j in range(1,n**m+1):
            a = mm([L[i-1]],mm(M,mtranspose([L[j-1]]),n),n)[0][0]
            z = vp(L[i-1],sm(a,L[j-1],n),n)
            for k in range(1,n**m+1):
                if tuple(z) == L[k-1]:
                    r.append(k)
        out.append(r)
    return out


def sympinv(PD,n,m,M):
    ###compute symnplectic quandle invariant###
    L = znm(n,m)
    T = symplecticquandle(n,m,M)
    t = Symbol('t')
    s = Symbol('s')
    ff = 0
    H = pdhomlist(PD,T)
    for h in H:
        h2 = subrack(h,T)
        SM = []
        for j in h2:        
            for i in L:
                if not i in SM and L[j-1]==i:
                    SM.append(i)
        ff = ff + t**(len(h2))*s**(len(submod(SM,n)))
    return ff    





############################
# Rack polynomials
############################



def rackpoly(M):
    s = Symbol('s')
    t = Symbol('t')
    f = 0
    for x in range(1,len(M)+1):
        c,r = 0,0 
        for i in range(1,len(M)+1):
            if M[x-1][i-1] == x: r = r+1 
            if M[i-1][x-1] == i: c = c+1
        f = f+s**c*t**r
    return f


def subrackpoly(S,M):
    s = Symbol('s')
    t = Symbol('t')
    f = 0
    for x in subrack(S,M):
        c,r = 0,0 
        for i in range(1,len(M)+1):
            if M[x-1][i-1] == x: r = r+1 
            if M[i-1][x-1] == i: c = c+1
        f = f+s**c*t**r
    return f

def rackopn(x,y,M,n):
    z = x
    for i in range(0,n):
        z=M[z-1][y-1]
    return z


def mnsubrackpoly(S,M,m,n):
    s = Symbol('s')
    t = Symbol('t')
    f = 0
    for x in subrack(S,M):
        c,r = 0,0 
        for i in range(1,len(M)+1):
            if rackopn(x,i,M,n) == x: r = r+1 
            if rackopn(i,x,M,m) == i: c = c+1
        f = f+s**c*t**r
    return f



def grackpolyinv1(G,T):
    ###Gauss code  rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q = Symbol('q')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        H = pdhomlist(gauss2pd(w),T)
        for h in H:
            ff = ff + q**(gausswrithe(w[0]) % N)*z**subrackpoly(h,T)
        w = gaussinc(w,1)
    return ff


def gmnrackpolyinv1(G,T,m,n):
    ###Gauss code rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q = Symbol('q')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        H = pdhomlist(gauss2pd(w),T)
        for h in H:
            ff = ff + q**(gausswrithe(w[0]) % N)*z**mnsubrackpoly(h,T,m,n)
        w = gaussinc(w,1)
    return ff



def grackpolyinv2(G,T):
    ### Gauss code rack counting invariant ###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            H = pdhomlist(gauss2pd(w),T)
            for h in H:
                ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*z**subrackpoly(h,T)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff


def gmnrackpolyinv2(G,T,m,n):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            H = pdhomlist(gauss2pd(w),T)
            for h in H:
                ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*z**mnsubrackpoly(h,T,m,n)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff



def grackpolyinv3(G,T):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                H = pdhomlist(gauss2pd(w),T)
                for h in H:
                    ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*z**subrackpoly(h,T)
                w = gaussinc(w,3) 
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff


def gmnrackpolyinv3(G,T,m,n):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                H = pdhomlist(gauss2pd(w),T)
                for h in H:
                    ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*z**mnsubrackpoly(h,T,m,n)
                w = gaussinc(w,3) 
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff




def grackpolyinv4(G,T):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    q4 = Symbol('q4')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for l in range(0,N):
                    H = pdhomlist(gauss2pd(w),T)
                    for h in H:
                        ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*q4**(gausswrithe(w[3]) % N)*z**subrackpoly(h,T)
                    w=gaussinc(w,4)   
                w = gaussinc(w,3) 
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff


def gmnrackpolyinv4(G,T,m,n):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    q4 = Symbol('q4')
    t = Symbol('t')
    s = Symbol('s')
    z = Symbol('z')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for l in range(0,N):
                    H = pdhomlist(gauss2pd(w),T)
                    for h in H:
                        ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*q4**(gausswrithe(w[3]) % N)*z**mnsubrackpoly(h,T,m,n)
                    w=gaussinc(w,4)   
                w = gaussinc(w,3) 
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff



def grackpolyinv(G,T):
    ###Gauss code rack counting invariant###
    k = len(G)
    if k == 1:
        out = grackpolyinv1(G,T)
    elif k == 2:
        out = grackpolyinv2(G,T)
    elif k == 3:
        out = grackpolyinv3(G,T)
    elif k == 4:
        out = grackpolyinv4(G,T)
    return out


def gmnrackpolyinv(G,T,m,n):
    ###Gauss code rack counting invariant###
    k = len(G)
    if k == 1:
        out = gmnrackpolyinv1(G,T,m,n)
    elif k == 2:
        out = gmnrackpolyinv2(G,T,m,n)
    elif k == 3:
        out = gmnrackpolyinv3(G,T,m,n)
    elif k == 4:
        out = gmnrackpolyinv4(G,T,m,n)
    return out





############################
# Gauss codes
############################





def gauss2pd(G):
    ###convert Gauss code to planar diagram###
    out = []
    semiarccount = 0
    complen = []
    for i in range(0,len(G)):
        semiarccount = semiarccount + len(G[i])
        complen.append(len(G[i]))
    nx, cr = [], [] 
    for k in range(0,semiarccount):
        nx.append(0)
        cr.append(0)
    current = 1
    i = 0
    for x in range(1,len(G)+1):
        for y in range(1,len(G[x-1])+1):
            i = i+1
            if y == complen[x-1]:
                nx[i-1] = current
            else:
                nx[i-1] = i+1
            j = 0
            for z in range(1,len(G)+1):
                for w in range(1,len(G[z-1])+1):
                    j = j +1
                    if G[z-1][w-1] + G[x-1][y-1] == 0:
                        cr[i-1] = j 
        current = current + complen[x-1]
    i = 0
    for x in range(1,len(G)+1):
        for y in range(1,len(G[x-1])+1):
            i = i+1
            if G[x-1][y-1] < 0:
                if G[x-1][y-1] % 1 == 0:
                    out.append([ 1,i,nx[cr[i-1]-1],nx[i-1],cr[i-1] ])
                else:
                    out.append([-1,i,cr[i-1],nx[i-1],nx[cr[i-1]-1] ])
    return out


def gaussinc(G,c):
    ###increment writhe in component c###
    M = 0
    for x in G:
        M = M + len(x)
    j = M/2 + 1
    out = []
    count = 0
    for x in G:
        if count == c-1:
            out.append(x+[-j,j])
        else:
            out.append(x)
        count = count + 1
    return out

def gausswrithe(g):
    ###self-crossing sum of component###
    out = 0
    for i in range(0,len(g)):
       if g[i]<0:
           for j in range(0,len(g)):
               if g[i]+g[j] == 0:
                   if type(g[i]) == int:
                       out = out + 1
                   if type(g[i]) == float:
                       out = out - 1
    return out


def gaussrackcount2(G,T):
    ###Gauss code rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    x = Symbol('x')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            x = x + quandlecount(gauss2pd(w),T)*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x
    

def gcoxinv2(G,n,m,a,M):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    L = znm(n,m)
    C = []
    for x in L:
        if gcd(mm(mm([list(x)],M,n),mtranspose([list(x)]),n)[0][0],n) == 1:
            C.append(tuple(x))
    T = coxrack(n,m,a,M)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    t = Symbol('t')
    s = Symbol('s')
    ff = Symbol('ff')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            H = pdhomlist(gauss2pd(w),T)
            for h in H:
                h2 = subrack(h,T)
                SM = []
                for j in h2:        
                    SM.append(C[j-1])
                ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*t**(len(h2))*s**(len(submod(SM,n)))
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff


def gaussrackcount3(G,T):
    ###Gauss code rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    x = Symbol('x')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                x = x + quandlecount(gauss2pd(w),T)*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x


def gaussrackcount4(G,T):
    ###Gauss code rack counting invariant###
    comp = len(G)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    q4 = Symbol('q4')
    x = Symbol('x')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for l in range(0,N):
                    x = x + quandlecount(gauss2pd(w),T)*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*q4**(gausswrithe(w[3]) % N)
                    w = gaussinc(w,4)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x
    


def gcoxinv3(G,n,m,a,M):
    ###Gauss code Coxeter rack counting invariant###
    comp = len(G)
    L = znm(n,m)
    C = []
    for x in L:
        if gcd(mm(mm([list(x)],M,n),mtranspose([list(x)]),n)[0][0],n) == 1:
            C.append(tuple(x))
    T = coxrack(n,m,a,M)
    N = rackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    t = Symbol('t')
    s = Symbol('s')
    ff = Symbol('ff')
    ff = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                H = pdhomlist(gauss2pd(w),T)
                for h in H:
                    h2 = subrack(h,T)
                    SM = []
                    for j in h2:        
                        SM.append(C[j-1])
                    ff = ff + q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*t**(len(h2))*s**(len(submod(SM,n)))
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return ff


def grackcount(G,T):
    ###Gauss code rack counting invariant###
    n = len(G)
    if n == 1:
        out = rackcountpoly(gauss2pd(G),T)
    elif n == 2:
        out = gaussrackcount2(G,T)
    elif n == 3:
        out = gaussrackcount3(G,T)
    elif n == 4:
        out = gaussrackcount4(G,T)
    return out



def gcoxinv(G,n,m,a,M):
    ###Gauss code rack counting invariant###
    k = len(G)
    if k == 1:
        out = coxinv(gauss2pd(G),n,m,a,M)
    elif k == 2:
        out = gcoxinv2(G,n,m,a,M)
    elif k == 3:
        out = gcoxinv3(G,n,m,a,M)
    return out



########################
# (t,s) racks
########################



def tsrack(n,t,s):
    ###(t,s) rack###
    if (s*(t+s-1)) % n != 0 or gcd(n,t) != 1:
        return False
    out = []
    for i in range(1,n+1):
        temp = []
        for j in range(1,n+1):
            if (t*i+s*j) % n == 0:
                temp.append(n)
            else:
                temp.append((t*i+s*j) % n)
        out.append(tuple(temp))
    return out


def nmtsrack(n,m,t,s):
    ###(t,s) rack###
    if (s*(t+s-1)) % n != 0 or gcd(n,t) != 1: return False
    L=znm(n,m)
    out = []
    for i in L:
        temp = []
        for j in L:
            k=vp(sm(t,i,n),sm(s,j,n),n)
            for l in range(1,len(L)+1):
                if tuple(k)==tuple(L[l-1]):
                    temp.append(l)
        out.append(tuple(temp))
    return out


def nmtsrackinv(G,n,m,t,s):
    ### nmstrack invariant ###
    if (s*(t+s-1)) % n != 0 or gcd(n,t) != 1: return False
    L=znm(n,m)
    T=nmtsrack(n,m,t,s)
    N=rackrank(T)
    q,z=Symbol('q'),Symbol('z')
    out=0
    g=G
    for w in range(0,N):
       p=pdhomlist(gauss2pd(g),T)
       for f in p:
           i = []
           for x in f:
               if not L[x-1] in i: i.append(L[x-1])
           I=submod(i,n)
           out=out+z**(len(I))*q**(gausswrithe(g) % N)
       g=gaussinc(g,1)
    return out




##########################
# Column Groups
##########################


def columngroup(L,M):
    ###compute column group for subrack L of M###
    L2 = subrack(L,M)
    g = []
    for x in L2:
        y = getcolumn(M,x)
        if not y in g:
            g.append(y)
    G = subgroup(g)
    return G


def cgpoly(PD,M):
    ###column group enhanced quandle invariant###
    u,v = Symbol('u'), Symbol('v')
    H = pdhomlist(PD,M)
    out = 0
    for f in H:
        I = subrack(f,M)
        O = columngroup(I,M)
        out = out+u**len(I)*v**(len(O))
    return out



def pdcgrackpoly(PD,M):
    ###column group enhanced rack invariant PD version###
    N = rackrank(M)
    q,u,v = Symbol('q'),Symbol('u'),Symbol('v')
    x = 0
    w = PD
    for n in range(0,N):
        H = pdhomlist(w,M)
        for f in H:
            I = subrack(f,M)
            O = columngroup(I,M)
            x = x + u**len(I)*v**(len(O))*q**(writhe(w) % N)
        w = pdwritheinc(w)
    return x



def gausscgrackpoly2(G,T):
    ###Gauss code rack counting invariant for 2-comp links###
    N = rackrank(T)
    q1,q2,u,v = Symbol('q1'),Symbol('q2'),Symbol('u'),Symbol('v')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            H = pdhomlist(gauss2pd(w),T)
            for f in H:
                I = subrack(f,T)
                O = columngroup(I,T)
                x = x + u**len(I)*v**(len(O))*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x


def gausscgrackpoly3(G,T):
    ###Gauss code rack counting invariant for 2-comp links###
    N = rackrank(T)
    q1,q2,q3,u,v = Symbol('q1'),Symbol('q2'),Symbol('q3'),Symbol('u'),Symbol('v')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                H = pdhomlist(gauss2pd(w),T)
                for f in H:
                    I = subrack(f,T)
                    O = columngroup(I,T)
                    x = x + u**len(I)*v**(len(O))*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x
    


def gcgrackpoly(G,T):
    ###Gauss code cg rack poly for 1,2,3 comp links###
    if len(G) == 1: return pdcgrackpoly(gauss2pd(G),T)
    if len(G) == 2: return gausscgrackpoly2(G,T)
    if len(G) == 3: return gausscgrackpoly3(G,T)



def bcolumngroup(L,B):
    ###compute column group for subbirack L of M###
    L2 = subbiq(L,B)
    g = []
    for x in L2:
        y = getcolumn(B[0],x)
        if not y in g:
            g.append(y)
        y = getcolumn(B[1],x)
        if not y in g:
            g.append(y)
    G = subgroup(g)
    return G



def gbcg1(G,M):
    ### blackboard birack column group polynomial 1 ###
    N=birackrank(M)
    v = Symbol('v')
    g,out = G,0
    for w in range(0,N):
        H = biqpdhomlist(gauss2pd(g),M)
        for f in H:
            C = bcolumngroup(subbiq(f,M),M)
            out = out+v**(len(C))
        G=gaussinc(G,1)
    return out


def gbcg2(G,M):
    ### blackboard birack column group polynomial 2 ###
    N = birackrank(M)
    v = Symbol('v')
    g,out = G,0
    for i in range(0,N):
        for j in range(0,N):
            H = biqpdhomlist(gauss2pd(g),M)
            for f in H:
                C = bcolumngroup(subbiq(f,M),M)
                out = out+v**(len(C))
            g = gaussinc(g,2)
        g = gaussinc(g,1)
    return out



def gbcg3(G,M):
    ### blackboard birack column group polynomial 3 ###
    N = birackrank(M)
    v = Symbol('v')
    g,x = G,0
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                H = biqpdhomlist(gauss2pd(g),M)
                for f in H:
                    C = bcolumngroup(subbiq(f,M),M)
                    out = out+v**(len(C))
                g=gaussinc(g,3)
            g = gaussinc(g,2)
        g = gaussinc(g,1)
    return x


def gbcg4(G,M):
    ### blackboard birack column group polynomial 4 ###
    N = birackrank(M)
    v = Symbol('v')
    g,x = G,0
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for l in range(0,N):
                    H = biqpdhomlist(gauss2pd(g),M)
                    for f in H:
                        C = bcolumngroup(subbiq(f,M),M)
                        out = out+v**(len(C))
                    g=gaussinc(g,4)
                g=gaussinc(g,3)
            g = gaussinc(g,2)
        g = gaussinc(g,1)
    return x
    

def gbcg(G,T):
    ###Gauss code birack counting invariant###
    n = len(G)
    if n == 1:
        out = gbcg1(G,T)
    elif n == 2:
        out = gbcg2(G,T)
    elif n == 3:
        out = gbcg3(G,T)
    elif n == 4:
        out = gbcg4(G,T)
    return out








def biqcolgroup2(L1,N):
    ###compute column group for subbiquandle L of N###
    U,L=N[0],N[1]
    Nb=barops(N)
    Ub,Lb=(Nb[0],Nb[1])
    L2 = subbiq(L1,N)
    g = []
    for x in L2:
        y = getcolumn(U,x)
        if not y in g:
            g.append(y)
        y = getcolumn(L,x)
        if not y in g:
            g.append(y)
        y = getcolumn(Ub,x)
        if not y in g:
            g.append(y)
        y = getcolumn(Lb,x)
        if not y in g:
            g.append(y)
    G = subgroup(g)
    return G





#######################################
# Biquandles and biracks
#######################################

def perminv(v):
    ###inverse of a permutation###
    out = []
    for i in range(1,len(v)+1):
       out.append(0)
    for i in range(1,len(v)+1):
       out[v[i-1]-1] = i
    return out


def pavail(v):
    ###List available entries###
    L = []
    for i in range(1,len(v)+1):
        if not i in v: L.append(i)
    return tuple(L)


def getcolumn(M,j):   # get column n from matrix M
    ###Gets column n from matrix M###
    c =[]
    for i in range(1,len(M)+1):
        c = c + [M[i-1][j-1]]
    return c


def tm(M):
    ###tuple matrix###
    out  = []
    for x in M: out.append(tuple(x))
    return tuple(out)


def lm(M):
    ###list matrix###
    out = []
    for x in M: out.append(list(x))
    return list(out)


def reptest(p):   # test whether p has repeated non-zero entries
    ###Test whether p has repeated non-zero entries###
    q = True
    L = []
    for i in range(1,len(p)+1):
        if p[i-1] != 0:
            if p[i-1] in L:
                q = False
            else:
                L.append(p[i-1])
    return q


def bfindzero(M):
    ###Find position of first zero entry in a birack matrix###
    for i in range(1,len(M[0])+1):
        for j in range(1,len(M[0])+1):
            if M[0][i-1][j-1] == 0:
                return (0,i,j)
            if M[1][i-1][j-1] == 0:
                 return (1,i,j)
    return False


def ybfill(M):
    ### fill with Yang-Baxter axiom ###
    U,L,n=lm(M[0]),lm(M[1]),len(M[0])
    keepgoing = True
    while keepgoing:
        keepgoing = False
        for i in range(1,n+1):
            for j in range(1,n+1):
                for k in range(1,n+1):
                    if U[i-1][j-1] != 0 and L[k-1][j-1] != 0 and U[i-1][L[k-1][j-1]-1] !=0 and U[j-1][k-1] != 0:
                        if U[U[i-1][j-1]-1][k-1] != 0 and U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1] == 0:
                            U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1] = U[U[i-1][j-1]-1][k-1]
                            keepgoing = True
                        if U[U[i-1][j-1]-1][k-1] == 0 and U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1] != 0:
                            U[U[i-1][j-1]-1][k-1] = U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1] 
                            keepgoing = True
                        if U[U[i-1][j-1]-1][k-1] != U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1]:
                            return False
                    if L[j-1][i-1] != 0 and U[i-1][j-1] != 0 and L[k-1][U[i-1][j-1]-1] !=0 and U[j-1][k-1] != 0 and L[k-1][j-1] != 0 and U[i-1][L[k-1][j-1]-1] != 0:
                        if U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1] != 0 and L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1] == 0:
                            L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1] = U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1]
                            keepgoing = True
                        if U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1] == 0 and L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1] != 0:
                            U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1] = L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1]
                            keepgoing = True
                        if U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1] != L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1]:
                            return False
                    if L[k-1][j-1] != 0 and L[j-1][i-1]!= 0 and U[i-1][j-1] != 0 and L[k-1][U[i-1][j-1]-1] != 0:
                        if L[L[k-1][j-1]-1][i-1] != 0 and L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1] == 0:
                            L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1] = L[L[k-1][j-1]-1][i-1]
                            keepgoing = True
                        if L[L[k-1][j-1]-1][i-1] == 0 and L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1] != 0:
                            L[L[k-1][j-1]-1][i-1] = L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1]
                            keepgoing = True
                        if L[L[k-1][j-1]-1][i-1] != L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1]:
                            return False
        for i in range(0,n):
            if not reptest(getcolumn(L,i+1)): return False
            if not reptest(getcolumn(U,i+1)): return False
    return ( tm(U),tm(L) ) 


def barops(M):
    ###get barred ops from birack matrix###
    U,L,n = M[0],M[1],len(M[0])
    Ub,Lb = zeromatrix(n), zeromatrix(n)
    for i in range(1,n+1):
        for j in range(1,n+1):
            Ub[U[i-1][j-1]-1][L[j-1][i-1]-1]=i 
            Lb[L[j-1][i-1]-1][U[i-1][j-1]-1]=j 
    return (Ub,Lb)



def invops(M):
    ###get barred ops from birack matrix###
    B1,B2,n = M[0],mtranspose(M[1]),len(M[0])
    Bi1,Bi2 = zeromatrix(n), zeromatrix(n)
    for i in range(1,n+1):
        for j in range(1,n+1):
            Bi1[B1[i-1][j-1]-1][B2[i-1][j-1]-1]=i 
            Bi2[B1[i-1][j-1]-1][B2[i-1][j-1]-1]=j 
    if bfindzero([Bi1,Bi2]): return False
    return (Bi1,Bi2)



def sideops(M):
    ### get sideways operations from birack matrix ###
    U,L=M[0], M[1]
    n=len(U)
    S1,S2=zeromatrix(n),zeromatrix(n)
    for i in range(1,n+1):
        for j in range(1,n+1):
            S1[i-1][L[j-1][i-1]-1]=j
            S2[i-1][L[j-1][i-1]-1]=U[i-1][j-1]
    if findzero(S1) or findzero(S2):
        return False
    return([S1,S2])


def biqlist(M):
    ###find all biquandles matching pattern###
    L = []
    L.append( ( tm(M[0]),tm(M[1]) ) )
    out = []
    while len(L)>0:
       N = (tm(L[0][0]),tm(L[0][1]))
       L[0:1] = []
       q = bfindzero(N)
       if q:
           for i in pavail(getcolumn(N[q[0]],q[2])):
               M2 = [lm(N[0]),lm(N[1])]
               M2[q[0]][q[1]-1][q[2]-1] = i
               M3 = ybfill(M2)
               if M3: L.append( ( tm(M3[0]),tm(M3[1]) ) )
       else:
           if biqtest(N):
               out.append( ( tm(N[0]), tm(N[1]) ) )
    return out



def bbrlist(M):
    ### find all blackboard biracks matching pattern ###
    L = []
    L.append( ( tm(M[0]),tm(M[1]) ) )
    out = []
    while len(L)>0:
       N = (tm(L[0][0]),tm(L[0][1]))
       L[0:1] = []
       q = bfindzero(N)
       if q:
           for i in pavail(getcolumn(N[q[0]],q[2])):
               M2 = [lm(N[0]),lm(N[1])]
               M2[q[0]][q[1]-1][q[2]-1] = i
               M3 = ybfill(M2)
               if M3: L.append( ( tm(M3[0]),tm(M3[1]) ) )
       else:
           out.append( ( tm(N[0]), tm(N[1]) ) )
    for x in out:
        if not bbrtest(x): out.remove(x)
    return out




def biqtest(M):
    U,L,n=M[0],M[1],len(M[0])
    Mb=barops(M)   
    if bfindzero(Mb): return False
    Ub,Lb = Mb[0],Mb[1]
    for i in range(1,n+1):
        ctp,ctm = 0,0
        for j in range(1,n+1):
            if U[j-1][i-1] == i and L[i-1][j-1] == j: ctp = ctp+1
            if Ub[j-1][i-1] == i and Lb[i-1][j-1] == j: ctm = ctm+1
            for k in range(1,n+1):
                if U[U[i-1][j-1]-1][k-1] != U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1]: return False
                if U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1] != L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1]: return False
                if L[L[k-1][j-1]-1][i-1] != L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1]: return False
                if Ub[Ub[i-1][j-1]-1][k-1] != Ub[Ub[i-1][Lb[k-1][j-1]-1]-1][Ub[j-1][k-1]-1]: return False
                if Ub[Lb[j-1][i-1]-1][Lb[k-1][Ub[i-1][j-1]-1]-1] != Lb[Ub[j-1][k-1]-1][Ub[i-1][Lb[k-1][j-1]-1]-1]: return False
                if Lb[Lb[k-1][j-1]-1][i-1] != Lb[Lb[k-1][Ub[i-1][j-1]-1]-1][Lb[j-1][i-1]-1]: return False
        if ctm != 1 or ctp != 1: return False
    return True



def diagonalmaps(M):
    ### get diagonal bijections; return False else ###
    S=sideops(M)
    S1,S2,n=S[0],S[1],len(S[0])
    x,y=[],[]
    for i in range(1,n+1):
        x.append(S1[i-1][i-1])
        y.append(S2[i-1][i-1])
    if reptest(x) and reptest(y):
        return [x,y]
    return False


def pi(M):
    ### type I move bijection positive twist ###
    D=diagonalmaps(M)
    return permcomp(perminv(D[0]),D[1])



def bbrtest(M):
    ### Test for blackboard birack axioms ###
    S0=sideops(M)
    S=[S0[0],mtranspose(S0[1])]
    if (not invops(M)) or (not invops(S)): return False
    if not diagonalmaps(M): return False
    if not ybtest(M): return False
    return True


def ybtest(M):
    ### Test whether M is a Yang-Baxter solution ###
    U,L,n=M[0],M[1],len(M[0])
    Mb=barops(M)
    Ub,Lb=Mb[0],Mb[1]
    for i in range(1,n+1):
        for j in range(1,n+1):
            for k in range(1,n+1):
                if U[U[i-1][j-1]-1][k-1] != U[U[i-1][L[k-1][j-1]-1]-1][U[j-1][k-1]-1]: return False
                if U[L[j-1][i-1]-1][L[k-1][U[i-1][j-1]-1]-1] != L[U[j-1][k-1]-1][U[i-1][L[k-1][j-1]-1]-1]: return False
                if L[L[k-1][j-1]-1][i-1] != L[L[k-1][U[i-1][j-1]-1]-1][L[j-1][i-1]-1]: return False
                if Ub[Ub[i-1][j-1]-1][k-1] != Ub[Ub[i-1][Lb[k-1][j-1]-1]-1][Ub[j-1][k-1]-1]: return False
                if Ub[Lb[j-1][i-1]-1][Lb[k-1][Ub[i-1][j-1]-1]-1] != Lb[Ub[j-1][k-1]-1][Ub[i-1][Lb[k-1][j-1]-1]-1]: return False
                if Lb[Lb[k-1][j-1]-1][i-1] != Lb[Lb[k-1][Ub[i-1][j-1]-1]-1][Lb[j-1][i-1]-1]: return False
    return True


def birackrank(M):
    ###find birack rank ###
    return permorder(pi(M))



def bisofill(M,N,f):
    ### fill isomorphism f:M -> N ###
    m,n = len(M[0]), len(N[0])
    if m != n: return False
    keepgoing = True
    while keepgoing:
        keepgoing = False
        for i in range(1,n+1):
            for j in range(1,n+1):
                if f[i-1] !=0 and f[j-1] != 0:
                    if f[M[0][i-1][j-1]-1] == 0: 
                        f[M[0][i-1][j-1]-1] = N[0][f[i-1]-1][f[j-1]-1]
                        keepgoing = True
                    if N[0][f[i-1]-1][f[j-1]-1] != f[M[0][i-1][j-1]-1]: return False
                    if f[M[1][i-1][j-1]-1] == 0: 
                        f[M[1][i-1][j-1]-1] = N[1][f[i-1]-1][f[j-1]-1]
                        keepgoing = True
                    if N[1][f[i-1]-1][f[j-1]-1] != f[M[1][i-1][j-1]-1]: return False
    if reptest(f): return f
    return False


def bhomfill(M,N,f):
    ### fill homomorphism f:M -> N ###
    m,n = len(M[0]), len(N[0])
    if m != n: return False
    keepgoing = True
    while keepgoing:
        keepgoing = False
        for i in range(1,n+1):
            for j in range(1,n+1):
                if f[i-1] !=0 and f[j-1] != 0:
                    if f[M[0][i-1][j-1]-1] == 0: 
                        f[M[0][i-1][j-1]-1] = N[0][f[i-1]-1][f[j-1]-1]
                        keepgoing = True
                    if N[0][f[i-1]-1][f[j-1]-1] != f[M[0][i-1][j-1]-1]: return False
                    if f[M[1][i-1][j-1]-1] == 0: 
                        f[M[1][i-1][j-1]-1] = N[1][f[i-1]-1][f[j-1]-1]
                        keepgoing = True
                    if N[1][f[i-1]-1][f[j-1]-1] != f[M[1][i-1][j-1]-1]: return False
    return f



def bhomlist(M,N):
    ###test for biquandle isomorphism###
    z = []
    for i in range(1,len(M[0])+1): z.append(0)
    L,out = [z],[]
    while len(L) > 0:
        w = L[0]
        L[0:1] = []  
        i = hfindzero(w)
        if not i: out.append(list(w))
        else:    
            for j in range(1,len(N[0])+1):
                phi = list(w)
                phi[i-1] = j
                v = bhomfill(M,N,phi)
                if v: L.append(tuple(v))
    return out

     

def bisotest(M,N):
    ###test for biquandle isomorphism###
    z = []
    for i in range(1,len(M[0])+1): z.append(0)
    L,out = [z],[]
    while len(L) != 0:
        w = L[0]
        L[0:1] = []  
        i = hfindzero(w)
        if (not i) and permtest(w): return True
        else:
            for j in pavail(w):
                phi = list(w)
                phi[i-1] = j
                v = bisofill(M,N,phi)
                if v: L.append(tuple(v))
    return False

     

def breducelist(L):
    ###Remove isomorphic copies from L###
    out = [tm(L[0])]
    W = L
    while len(W)>0:
        x = W[0]
        W[0:1] = []
        newbiq = True
        for y in out:
            if bisotest(x,y): newbiq = False
        if newbiq: out.append(tm(x))
    return out



def biqpdhomlist(PD,N):   
    ###Lists homomorphisms from knot biquandle/birack of planar diagram###
    z = []
    for i in range(1,2*len(PD)+1):
        z = z + [0]
    L = [z]
    out = []
    while len(L) != 0:
        w = L[0]
        L[0:1] = []  
        if w:
            i = hfindzero(w)
            if not i:
                out.append(w)
            else:
                for j in range(1,len(N[0])+1):
                    phi = list(w)
                    phi[i-1] = j
                    v = biqpdhomfill(PD,N,phi)
                    if v: L.append(tuple(v))
    return out     


def biqpdhomfill(PD,N,phi):   # fill in homomorphism
    ###Fills in entries in a biquandle homomorphism###
    U,L=N[0],N[1]
    Nb=barops(N)
    Ub,Lb=(Nb[0],Nb[1])
    f = phi
    c = True
    out = True
    while c == True:
        c = False
        for X in PD:
            if X[0] == 1:
                if f[X[1]-1] != 0 and f[X[4]-1] != 0:
                    if f[X[3]-1] == 0: 
                        f[X[3]-1] = U[f[X[1]-1]-1][f[X[4]-1]-1]
                        c = True
                    elif f[X[3]-1] != U[f[X[1]-1]-1][f[X[4]-1]-1]:
                        out = False
                    if f[X[2]-1] == 0: 
                        f[X[2]-1] = L[f[X[4]-1]-1][f[X[1]-1]-1]
                        c = True
                    elif f[X[2]-1] != L[f[X[4]-1]-1][f[X[1]-1]-1]:
                        out = False
            if X[0] == -1:
                if f[X[1]-1] != 0 and f[X[2]-1] != 0:
                    if f[X[3]-1] == 0: 
                        f[X[3]-1] = Ub[f[X[1]-1]-1][f[X[2]-1]-1]
                        c = True
                    elif f[X[3]-1] != Ub[f[X[1]-1]-1][f[X[2]-1]-1]:
                        out = False
                    if f[X[4]-1] == 0: 
                        f[X[4]-1] = Lb[f[X[2]-1]-1][f[X[1]-1]-1]
                        c = True
                    elif f[X[4]-1] != Lb[f[X[2]-1]-1][f[X[1]-1]-1]:
                        out = False                
    if out == True: 
       return f
    else:
       return out



def subbiq(L1,N):
    ###subbiquandle gen by L1###
    U,L = N[0],N[1]
    Nb = barops(N)
    Ub,Lb = (Nb[0],Nb[1])
    out = []
    for x in L1: 
        if not x in out:
            out.append(x)
    keepgoing = True
    while keepgoing:
        keepgoing = False
        for i in out: 
            for j in out:
                if not U[i-1][j-1] in out: 
                    out.append(U[i-1][j-1])
                    keepgoing = True
                if not L[i-1][j-1] in out: 
                    out.append(L[i-1][j-1])
                    keepgoing = True
                if not Ub[i-1][j-1] in out: 
                    out.append(Ub[i-1][j-1])
                    keepgoing = True
                if not Lb[i-1][j-1] in out: 
                    out.append(Lb[i-1][j-1])
                    keepgoing = True                
    return out




def abq(n,s,t):
    ###Alexander biquandle Z_n###
    if not invn(t,n): return False
    U,L=zeromatrix(n),zeromatrix(n)
    for i in range(1,n+1):
       for j in range(1,n+1):
           U[i-1][j-1] = (t*i+(1-s*t)*j) %n 
           L[i-1][j-1] = (s*i %n) 
    for i in range(1,n+1):
       for j in range(1,n+1):
           if U[i-1][j-1] == 0: U[i-1][j-1] = n
           if L[i-1][j-1] == 0: L[i-1][j-1] = n
    return([U,L])



def tsr(n,m,t,s,r):
    ### (t,s,r) birack over (Z_n)^m ###
    if (not invn(t,n) or not invn(r,n) or (s**2 %n) != (s*(1-t*r) %n)): 
        return False
    V=znm(n,m)
    U,L=zeromatrix(n**m),zeromatrix(n**m)
    for i in range(1,len(V)+1):
        for j in range(1,len(V)+1):
            u=vp(sm(t,V[i-1],n),sm(s,V[j-1],n),n)
            l=sm(r,V[j-1],n)
            for k in range(1,len(V)+1):
                if tuple(V[k-1])==tuple(u): U[i-1][j-1] = k
                if tuple(V[k-1])==tuple(l): L[j-1][i-1] = k
    return(U,L)


def gbbrcount1(g,M):
    ### blackboard birack counting polynomial ###
    N = birackrank(M)
    q = Symbol('q')
    x = Symbol('x')
    x = 0
    w = g
    for n in range(0,N):
        x = x + len(biqpdhomlist(gauss2pd(w),M))*q**(gausswrithe(w[0])% N)
        w = gaussinc(w,1)
    return x


def gaussbbrcount2(G,T):
    ### Gauss code birack counting invariant###
    comp = len(G)
    N = birackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    x = Symbol('x')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            x = x + len(biqpdhomlist(gauss2pd(w),T))*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x



def gaussbbrcount3(G,T):
    ###Gauss code birack counting invariant###
    comp = len(G)
    N = birackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    x = Symbol('x')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                x = x + len(biqpdhomlist(gauss2pd(w),T))*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x

def gaussbbrcount4(G,T):
    ###Gauss code rack counting invariant###
    comp = len(G)
    N = birackrank(T)
    q1 = Symbol('q1')
    q2 = Symbol('q2')
    q3 = Symbol('q3')
    q4 = Symbol('q4')
    x = Symbol('x')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for l in range(0,N):
                    x = x + len(biqpdhomlist(gauss2pd(w),T))*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*q4**(gausswrithe(w[3]) % N)
                    w = gaussinc(w,4)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x
    

def gbbrcount(G,T):
    ###Gauss code birack counting invariant###
    n = len(G)
    if n == 1:
        out = gbbrcount1(G,T)
    elif n == 2:
        out = gaussbbrcount2(G,T)
    elif n == 3:
        out = gaussbbrcount3(G,T)
    elif n == 4:
        out = gaussbbrcount4(G,T)
    return out



def subbbr(L1,N):
    ###subbbr gen by L1###
    U,L = N[0],N[1]
    Nb = barops(N)
    S = sideops(N)
    Ub,Lb,S1,S2 = Nb[0],Nb[1],S[0],S[1]
    out = []
    for x in L1: 
        if not x in out:
            out.append(x)
    keepgoing = True
    while keepgoing:
        keepgoing = False
        for i in out: 
            for j in out:
                if not U[i-1][j-1] in out: 
                    out.append(U[i-1][j-1])
                    keepgoing = True
                if not L[i-1][j-1] in out: 
                    out.append(L[i-1][j-1])
                    keepgoing = True
                if not Ub[i-1][j-1] in out: 
                    out.append(Ub[i-1][j-1])
                    keepgoing = True
                if not Lb[i-1][j-1] in out: 
                    out.append(Lb[i-1][j-1])
                    keepgoing = True 
                if not S1[i-1][j-1] in out: 
                    out.append(S1[i-1][j-1])
                    keepgoing = True
                if not S2[i-1][j-1] in out: 
                    out.append(S2[i-1][j-1])
                    keepgoing = True 
    return out


def bbrpoly(N):
    ### birack polynomial ###
    U,L=N[0],N[1]
    s1,s2,t1,t2=Symbol('s1'),Symbol('s2'),Symbol('t1'),Symbol('t2')
    out = 0
    for x in range(1,len(U)+1):
        c1,c2,r1,r2=0,0,0,0
        for y in range(1,len(U)+1):
            if U[y-1][x-1] == y: c1 = c1+1
            if L[y-1][x-1] == y: c2 = c2+1
            if U[x-1][y-1] == x: r1 = r1+1
            if L[x-1][y-1] == x: r2 = r2+1
        out = out + s1**c1*s2**c2*t1**r1*t2**r2;
    return out


def subbbrpoly(L1,N):
    ### subbirack polynomial ###
    U,L=N[0],N[1]
    s1,s2,t1,t2=Symbol('s1'),Symbol('s2'),Symbol('t1'),Symbol('t2')
    out = 0
    S=subbbr(L1,N)
    for x in S:
        c1,c2,r1,r2=0,0,0,0
        for y in range(1,len(U)+1):
            if U[y-1][x-1] == y: c1 = c1+1
            if L[y-1][x-1] == y: c2 = c2+1
            if U[x-1][y-1] == x: r1 = r1+1
            if L[x-1][y-1] == x: r2 = r2+1
        out = out + s1**c1*s2**c2*t1**r1*t2**r2;
    return out


def gbbrcountpoly1(g,T):
    ### blackboard birack counting polynomial enhanced polynomial ###
    N = birackrank(T)
    z,q,s1,s2,t1,t2 = Symbol('z'),Symbol('q'),Symbol('s1'),Symbol('s2'),Symbol('t1'),Symbol('t2')
    x = 0
    w = g
    for n in range(0,N):
        for f in biqpdhomlist(gauss2pd(w),T):
            x = x +z**(subbbrpoly(f,T))#*q**(gausswrithe(w[0])% N)
        w = gaussinc(w,1)
    return x


def gaussbbrcountpoly2(G,T):
    ### Gauss code birack counting polynomial enhanced invariant###
    comp = len(G)
    N = birackrank(T)
    z,q1,q2 = Symbol('z'),Symbol('q1'),Symbol('q2')
    s1,s2,t1,t2 = Symbol('s1'),Symbol('s2'),Symbol('t1'),Symbol('t2')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for f in biqpdhomlist(gauss2pd(w),T):
                x = x + z**(subbbrpoly(f,T))#*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x



def gaussbbrcountpoly3(G,T):
    ###Gauss code birack counting polynomial enhanced  invariant###
    comp = len(G)
    N = birackrank(T)
    q1,q2,q3 = Symbol('q1'),Symbol('q2'),Symbol('q3')
    s1,s2,t1,t2 = Symbol('s1'),Symbol('s2'),Symbol('t1'),Symbol('t2')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for f in biqpdhomlist(gauss2pd(w),T):
                    x = x + z**(subbbrpoly(f,T))#*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x

def gaussbbrcountpoly4(G,T):
    ###Gauss code rack counting polynomial enhanced invariant###
    comp = len(G)
    N = birackrank(T)
    q1 = Symbol('q1'),Symbol('q2'),Symbol('q3'),Symbol('q4')
    s1,s2,t1,t2 = Symbol('s1'),Symbol('s2'),Symbol('t1'),Symbol('t2')
    x = 0
    w = G
    for i in range(0,N):
        for j in range(0,N):
            for k in range(0,N):
                for l in range(0,N):
                    for f in biqpdhomlist(gauss2pd(w),T):
                        x = x + z**(subbbrpoly(f,T))#*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*q4**(gausswrithe(w[3]) % N)
                    w = gaussinc(w,4)
                w = gaussinc(w,3)
            w = gaussinc(w,2)
        w = gaussinc(w,1)
    return x


def gbbrpoly(G,T):
    ###Gauss code birack counting invariant###
    n = len(G)
    if n == 1:
        out = gbbrcountpoly1(G,T)
    elif n == 2:
        out = gaussbbrcountpoly2(G,T)
    elif n == 3:
        out = gaussbbrcountpoly3(G,T)
    elif n == 4:
        out = gaussbbrcountpoly4(G,T)
    return out





##############################
# R-Shadows
##############################


def rshadowtest(R,X):
    ### Test whether X is an R-shadow ###
    x,r=len(X),len(R)
    for i in range(1,x+1):
        if not reptest(getcolumn(X,i)): return False
        for j in range(1,r+1):
            for k in range(1,r+1):
                if X[X[i-1][j-1]-1][k-1] != X[X[i-1][k-1]-1][R[j-1][k-1]-1]: return False
    return True



def rshadowfill(R,Y):
    ### fill entries in an R-shadow matrix ###
    X=lm(Y)
    x,r=len(X),len(R)
    c = True
    while c:
        c = False
        for i in range(1,x+1):
            for j in range(1,r+1):
                if X[i-1][j-1] !=0:
                    if not reptest(getcolumn(X,j)): 
                        return False
                    for k in range(1,r+1):
                        if X[i-1][k-1]!= 0:
                            if X[X[i-1][j-1]-1][k-1] != 0 and X[X[i-1][k-1]-1][R[j-1][k-1]-1] == 0:
                                X[X[i-1][k-1]-1][R[j-1][k-1]-1] = X[X[i-1][j-1]-1][k-1]
                                c = True
                            if X[X[i-1][j-1]-1][k-1] == 0 and X[X[i-1][k-1]-1][R[j-1][k-1]-1] != 0:
                                X[X[i-1][j-1]-1][k-1] = X[X[i-1][k-1]-1][R[j-1][k-1]-1]
                                c = True
                            if X[X[i-1][j-1]-1][k-1] !=  X[X[i-1][k-1]-1][R[j-1][k-1]-1]:
                                return False
    return tm(X)  



def rshadowfind(R,n):
    ### Find all R-shadows of size n ###
    Z=nmzeromatrix(n,len(R))
    w = [tm(Z)]
    out = []
    while len(w) >0:
       X = w[0]
       w[0:1] = []
       f = findzero(X)
       if f: 
           for i in pavail(getcolumn(X,f[1])):
               N = lm(X)
               N[f[0]-1][f[1]-1] = i
               M = rshadowfill(R,N)
               if M: 
                   w.append(tm(M))
       else: out.append(tm(X))
    return out 





def gsrpinv1(g,R,X):
    ### shadow rack coloring invariant ###
    q,z,t,u=Symbol('q'),Symbol('z'),Symbol('t'),Symbol('u')
    N=rackrank(R)
    out=0
    w=g
    for k in range(0,N):
        L=pdhomlist(gauss2pd(w),R)
        for f in L:
            Im=subrack(f,R)
            for a in range(1,len(X)+1):
                Sr=[]
                for i in Im:
                    if not X[a-1][i-1] in Sr:
                        Sr.append(X[a-1][i-1])
                out = out+ z**srp(Sr,Im,X)*q**(gausswrithe(w[0]) % N)
        w=gaussinc(w,1)
    return(out)


def gsrpinv2(g,R,X):
    ### shadow rack coloring invariant ###
    q1,q2,z,t,u=Symbol('q1'),Symbol('q2'),Symbol('z'),Symbol('t'),Symbol('u')
    N=rackrank(R)
    out=0
    w=g
    for k1 in range(0,N):
        for k in range(0,N):
            L=pdhomlist(gauss2pd(w),R)
            for f in L:
                Im=subrack(f,R)
                for a in range(1,len(X)+1):
                    Sr=[]
                    for i in Im:
                        if not X[a-1][i-1] in Sr:
                            Sr.append(X[a-1][i-1])
                    out = out+ z**srp(Sr,Im,X)*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)
            w=gaussinc(w,1)
        w=gaussinc(w,2)
    return(out)


def gsrpinv3(g,R,X):
    ### shadow rack coloring invariant ###
    q1,q2,q3,z,t,u=Symbol('q1'),Symbol('q2'),Symbol('q3'),Symbol('z'),Symbol('t'),Symbol('u')
    N=rackrank(R)
    out=0
    w=g
    for k2 in range(0,N):
        for k1 in range(0,N):
            for k in range(0,N):
                L=pdhomlist(gauss2pd(w),R)
                for f in L:
                    Im=subrack(f,R)
                    for a in range(1,len(X)+1):
                        Sr=[]
                        for i in Im:
                            if not X[a-1][i-1] in Sr:
                                Sr.append(X[a-1][i-1])
                        out = out+ z**srp(Sr,Im,X)*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)
                w=gaussinc(w,1)
            w=gaussinc(w,2)
        w=gaussinc(w,3)
    return(out)


def gsrpinv4(g,R,X):
    ### shadow rack coloring invariant ###
    q1,q2,q3,q4,z,t,u=Symbol('q1'),Symbol('q2'),Symbol('q3'),Symbol('q4'),Symbol('z'),Symbol('t'),Symbol('u')
    N=rackrank(R)
    out=0
    w=g
    for k3 in range(0,N):
        for k2 in range(0,N):
            for k1 in range(0,N):
                for k in range(0,N):
                    L=pdhomlist(gauss2pd(w),R)
                    for f in L:
                        Im=subrack(f,R)
                        for a in range(1,len(X)+1):
                            Sr=[]
                            for i in Im:
                                if not X[a-1][i-1] in Sr:
                                    Sr.append(X[a-1][i-1])
                            out = out+ z**srp(Sr,Im,X)*q1**(gausswrithe(w[0]) % N)*q2**(gausswrithe(w[1]) % N)*q3**(gausswrithe(w[2]) % N)*q4**(gausswrithe(w[3]) % N)
                    w=gaussinc(w,1)
                w=gaussinc(w,2)
            w=gaussinc(w,3)
        w==gaussinc(w,4)
    return(out)




def gsrpinv(G,T,S):
    ### Gauss code shadow rack counting invariant ###
    n = len(G)
    if n == 1:
        out = gsrpinv1(G,T,S)
    elif n == 2:
        out = gsrpinv2(G,T,S)
    elif n == 3:
        out = gsrpinv3(G,T,S)
    elif n == 4:
        out = gsrpinv4(G,T,S)
    return out




def srp(S,R,X):
    ###shadow rack polynomial###
    out = 0
    t = Symbol('t')
    for x in S:
        c = 0
        for r in R:
            if X[x-1][r-1] == x:
                c = c+1
        out = out + t**c
    return out


############################################################
#  Rack Modules
############################################################


Z=Symbol('Z')

def rmodgens(X):
    ### Generators of rack algebra of R ###
    n,N=len(X),rackrank(X)
    temp=[]
    count=1
    for k in range(1,n+1):
        row=[]
        for j in range(1,n+1):
            row.append(count)
            count=count+1
        temp.append(row)
    s,t=tm(temp),lm(temp)
    for x in range(1,n+1):
        for y in range(1,n+1):
            for z in range(1,n+1):
                if X[X[x-1][y-1]-1][z-1]==X[X[x-1][z-1]-1][X[y-1][z-1]-1]:
                    if t[x-1][y-1]!=t[x-1][z-1]:
                        txz,txy=t[x-1][z-1],t[x-1][y-1]
                        if txy<txz:
                            t=lm(t)
                            t=genreplace(t,txz,txy)
                            t=tm(t)
                        elif txy>txz:
                            t=lm(t)
                            t=genreplace(t,txy,txz)
                            t=tm(t)
                if X[X[x-1][y-1]-1][z-1]==X[x-1][z-1]:
                    if t[x-1][y-1]!=t[X[x-1][z-1]-1][X[y-1][z-1]-1]:
                        txy,txzyz=t[x-1][y-1],t[X[x-1][z-1]-1][X[y-1][z-1]-1]
                        if txy<txzyz:
                            t=lm(t)
                            t=genreplace(t,txzyz,txy)
                            t=tm(t)
                        elif txy>txzyz:
                            t=lm(t)
                            t=genreplace(t,txy,txzyz)
                            t=tm(t)
                if X[x-1][y-1]==X[X[x-1][z-1]-1][X[y-1][z-1]-1]:
                    txyz,txz=t[X[x-1][y-1]-1][z-1],t[x-1][z-1]
                    if txyz!=txz:
                        if txyz<txz:
                            t=lm(t)
                            t=genreplace(t,txz,txyz)
                            t=tm(t)
                        elif txyz>txz:
                            t=lm(t)
                            t=genreplace(t,txyz,txz)
                            t=tm(t)
                if X[x-1][y-1]==X[x-1][z-1]:
                    if t[X[x-1][y-1]-1][z-1]!=t[X[x-1][z-1]-1][X[y-1][z-1]-1]:
                        txyz,txzyz=t[X[x-1][y-1]-1][z-1],t[X[x-1][z-1]-1][X[y-1][z-1]-1]
                        if txyz<txzyz:
                            t=lm(t)
                            t=genreplace(t,txzyz,txyz)
                            t=tm(t)
                        elif txyz>txzyz:
                            t=lm(t)
                            t=genreplace(t,txyz,txzyz)
                            t=tm(t)
    for x in range(1,n+1):
        for y in range(1,n+1):
            for z in range(1,n+1):
                if t[X[x-1][y-1]-1][z-1]==t[y-1][z-1]:
                    sxy,sxzyz=s[x-1][y-1],s[X[x-1][z-1]-1][X[y-1][z-1]-1]
                    if sxy!=sxzyz:
                        if sxy<sxzyz:
                            s=lm(s)
                            s=genreplace(s,sxzyz,sxy)
                            s=tm(s)
                        if sxy>sxzyz:
                            s=lm(s)
                            s=genreplace(s,sxy,sxzyz)
                            s=tm(s)
    return t,s




def genreplace(X,x,y):
    ### replace entry x with y ###
    out=[]
    for i in range(0,len(X)):
        row=[]
        for j in range(0,len(X[0])):
            if X[i][j]==x:
                row.append(y)
            else:
                row.append(X[i][j])
        out.append(row)
    return(out)
    

def zmatrix(n):
    ### return n by n Z matrix### 
    out = []
    for i in range(0,n): 
        r = []
        for j in range(0,n): r.append(Z)
        out.append(r)
    return out


def findz(M):
    ###Find position of first Z entry in a rack module matrix###
    for i in range(1,len(M[0])+1):
        for j in range(1,len(M[0])+1):
            if M[0][i-1][j-1] == Z:
                return (0,i,j)
            if M[1][i-1][j-1] == Z:
                 return (1,i,j)
    return False


def hfindz(w):
    ###Find position of first Z entry in a coloring###
    for i in range(0,len(w)):
        if w[i] == Z:
            return i
    return -1



def rackmodfind(X,n):
    ### Find rack module structures on Z_n ###
    m=len(X)
    working,out=[[tm(zmatrix(m)),tm(zmatrix(m))]],[]
    while len(working)!=0:
        c=working[0]
        working[0:1]=[]
        w=findz(c)
        if w:
            for k in range(0,n):
                wk=[lm(c[0]),lm(c[1])]
                wk[w[0]][w[1]-1][w[2]-1]=k
                wk2=False
                if (w[0]==0 and invn(k,n)!=False) or w[0]==1:
                    wk2=rackmodfill(wk,X,n)
                if wk2: 
                    working.append((tm(wk2[0]),tm(wk2[1])))
        else: 
            if Ntest(c,X,n):
                out.append(c)
    return out


def Ntest(M,X,n):
    ### Test for rank-N condition ###
    m=len(X)
    t,s=M[0],M[1]
    N=rackrank(X)
    for x in range(1,m+1):
        p,y=1,x
        for k in range(0,N):
            p=(p*(t[y-1][y-1]+s[y-1][y-1]))%n
            y=X[y-1][y-1]
        if p!=1: 
            return False
    return True


def rackmodfill(M,X,n):
    ### Fill rack module entries ###
    m=len(X)
    t,s=lm(M[0]),lm(M[1])
    for x in range(1,m+1):
        for y in range(1,m+1):
            for z in range(1,m+1):
# t-fill
                if t[X[x-1][y-1]-1][z-1]!=Z and t[x-1][y-1]!=Z and t[X[x-1][z-1]-1][X[y-1][z-1]-1]!=Z and t[x-1][z-1]!=Z:
                    if ((t[X[x-1][y-1]-1][z-1]*t[x-1][y-1]) %n) != ((t[X[x-1][z-1]-1][X[y-1][z-1]-1]*t[x-1][z-1]) %n):
                        return False
                if t[X[x-1][y-1]-1][z-1]==Z:
                    if t[x-1][y-1]!=Z and t[X[x-1][z-1]-1][X[y-1][z-1]-1]!=Z and t[x-1][z-1]!=Z:
   ing[0]
        working[0:1]=[]
        z=hfindz(w)
        if z==-1:
            out.append(w)
        else:
            for k in range(0,n):
                w1=list(w)
                w1[z]=k
                w2=rmodhfill(w1,v,PD,X,M,n)
                if w2:
                    working.append(tuple(w2))
    return out



def rmodhfill(w,v,PD,X,M,n):
    ### fill rmod coloring ###
    t,s=M[0],M[1]
    w1,c=list(w),True
    while c:
        c=False
        for x in PD:
            if x[0]==1:
                if w1[x[2]-1]==Z and w1[x[4]-1]!=Z:
                    w1[x[2]-1]=w1[x[4]-1]
                    c=True
                if w1[x[4]-1]==Z and w1[x[2]-1]!=Z:
                    w1[x[4]-1]=w1[x[2]-1]
                    c=True
                if w1[x[2]-1]!=Z and w1[x[4]-1]!=Z:
                    if w1[x[2]-1]!=w1[x[4]-1]: 
                        return False
                if w1[x[3]-1]==Z and w1[x[1]-1]!=Z and w1[x[4]-1]!=Z:
                    w1[x[3]-1]=(t[v[x[1]-1]-1][v[x[4]-1]-1]*w1[x[1]-1]+s[v[x[1]-1]-1][v[x[4]-1]-1]*w1[x[4]-1])%n
                    c=True
                if w1[x[3]-1]!=Z and w1[x[1]-1]!=Z and w1[x[4]-1]!=Z:
                    if w1[x[3]-1]!=(t[v[x[1]-1]-1][v[x[4]-1]-1]*w1[x[1]-1]+s[v[x[1]-1]-1][v[x[4]-1]-1]*w1[x[4]-1])%n:
                        return False
            if x[0]==-1:
                if w1[x[2]-1]==Z and w1[x[4]-1]!=Z:
                    w1[xx[2]-1]=w1[x[4]-1]
                    c=True
                if w1[x[4]-1]==Z and w1[x[2]-1]!=Z:
                    w1[x[4]-1]=w1[x[2]-1]
                    c=True
                if w1[x[2]-1]!=Z and w1[x[4]-1]!=Z:
                    if w1[x[2]-1]!=w1[x[4]-1]: 
                        return False
                if w1[x[1]-1]==Z and w1[x[3]-1]!=Z and w1[x[4]-1]!=Z:
                    w1[x[1]-1]=(t[v[x[3]-1]-1][v[x[4]-1]-1]*w1[x[3]-1]+s[v[x[3]-1]-1][v[x[4]-1]-1]*w1[x[4]-1])%n
                    c=True
                if w1[x[1]-1]!=Z and w1[x[3]-1]!=Z and w1[x[4]-1]!=Z:
                    if w1[x[1]-1]!=(t[v[x[3]-1]-1][v[x[4]-1]-1]*w1[x[3]-1]+s[v[x[3]-1]-1][v[x[4]-1]-1]*w1[x[4]-1])%n:
                        return False
    return w1



def grackmodinv1(G,X,M,n):
    ### rack module invariant ###
    N=rackrank(X)
    g,u=G,Symbol('u')
    out=0
    for k in range(0,N):
        p=gauss2pd(g)
        L=pdhomlist(p,X)
        for x in L:
            out=out+u**(len(rmodcolor(x,p,X,M,n)))
        g=gaussinc(g,1)
    return out


def grackmodinv12(G,X,M,n):
    ### rack module invariant ###
    N=rackrank(X)
    g,u=G,Symbol('u')
    out=0
    for k1 in range(0,N):
        for k in range(0,N):
            p=gauss2pd(g)
            L=pdhomlist(p,X)
            for x in L:
                out=out+u**(len(rmodcolor(x,p,X,M,n)))
            g=gaussinc(g,1)
        g=gaussinc(g,2)
    return out


def grackmodinv13(G,X,M,n):
    ### rack module invariant ###
    N=rackrank(X)
    g,u=G,Symbol('u')
    out=0
    for k2 in range(0,N):
        for k3 in range(0,N):
            for k in range(0,N):
                p=gauss2pd(g)
                L=pdhomlist(p,X)
                for x in L:
                    out=out+u**(len(rmodcolor(x,p,X,M,n)))
                g=gaussinc(g,1)
            g=gaussinc(g,2)
        g=gaussinc(g,3)
    return out

def grackmodinv14(G,X,M,n):
    ### rack module invariant ###
    N=rackrank(X)
    g,u=G,Symbol('u')
    out=0
    for k4 in range(0,N):
        for k2 in range(0,N):
            for k3 in range(0,N):
                for k in range(0,N):
                    p=gauss2pd(g)
                    L=pdhomlist(p,X)
                    for x in L:
                        out=out+u**(len(rmodcolor(x,p,X,M,n)))
                    g=gaussinc(g,1)
                g=gaussinc(g,2)
            g=gaussinc(g,3)
        g=gaussinc(g,4)
    return out



def grackmodinv(G,X,M1,n):
    ### rack module invariant ###
    if len(G)==1:
        return grackmodinv1(G,X,M1,n)
    if len(G)==2:
        return grackmodinv12(G,X,M1,n)
    if len(G)==3:
        return grackmodinv13(G,X,M1,n)
    if len(G)==4:
        return grackmodinv14(G,X,M1,n)