#################################################################
#################################################################
#                                                               #
#  Written by Sam Nelson. Permission granted to use, adapt,     #
#  reimplement and redistribute as desired.                     #
#                                                               #
#  Tested on Maple versions 7 and 9.5.                          #
#                                                               #
#  Report any bugs to knots -at- esotericka.org                 #
#                                                               #
#################################################################
#################################################################

with(linalg);
with(combinat);

permtest := proc (v)
#
# # Used in qtest
# #
# # tests whether a vector represents a permutation
#
ret:=true;
for i from 1 to coldim(matrix([v])) do
 for j from 1 to i-1 do
  if v[i]=v[j] then 
   ret:=false; 
  fi;
 od;
od;
#
ret;
end;


nextperm := proc (v)
#
# # Returns next permutation in dictionary order
#
if permtest(v) then
 R:=permute(coldim(matrix([v])));
 ret:=R[coldim(matrix([v]))!];
 for i from 1 to coldim(matrix([v]))!-1 do
  if R[i]=eval(v) then
   ret:=R[i+1];
   break;
  fi;
 od;
else ret:=false;
fi;
eval(ret);
end;

nextmap := proc (v,n)
#
# # returns next m-component vector in dictionary order 
# # with entries 1..n
#
m:=coldim(matrix([v]));
R:=listmaps(m,n);
for i from 1 to n^m do
 if v=R[i] then 
  w:=R[i+1];
  break;
 fi;
od;
eval(w);
end;


listmaps := proc(a,b)
#
# # list all a-component vectors with entries in 1..b
#
w:=initperm(b);
v:=[];
for i from 1 to a do
 v:=[op(v),w];
od;
ret:=[];
T:=cartprod(v);
while not T[finished] do
 t:=T[nextvalue]();
 ret:=[op(ret),t];
od;
eval(ret);
end;


homtest := proc(A, B, v)
#
# # detemines whether v:A->B is a quandle homomorphism
#
n:=coldim(A);
m:=coldim(B);
ret:=true;
#
if n<>coldim(matrix([v])) then ret:=false; 
else
 for i from 1 to n do
  for j from 1 to n do
    if A[i,j]<>0 and v[A[i,j]]<> B[v[i],v[j]] then 
     ret:=false; 
     break; 
    fi;
  od;
 od;
fi;
#
ret;
end;


homlist := proc(A,B)
#
# # lists all quandle homomorphisms from A to B
#
n:=coldim(A);
m:=coldim(B);
w:=initmap(n);
ret:=[]; 
ret:=[eval(w)];
#
for i from 2 to m^n do
 v:=nextmap(w,m);
 if homtest(A,B,v) then 
  ret:=[op(ret),v]; 
 fi;
 w:=v;
od;
#
eval(ret);
end;


homcount := proc(A,B)
#
# # Counts homomorphisms from A to B
#
n:=coldim(matrix([homlist2(A,B)]));
#
eval(n);
end;


initperm:= proc (n)
#
# # outputs the permutation vector [1,2,...,n]
#
w:=[];
for i from 1 to n do 
 w:=[op(w),i]; 
od;
#
eval(w);
end;


initmap:= proc (n)
#
# # outputs an n-component vector of 1s [1,1,...,1]
#
w:=[];
for i from 1 to n do 
 w:=[op(w),1]; 
od;
#
eval(w);
end;

autlist := proc(A)
#
# # Lists automorphisms of a quandle A
#
n:=coldim(A);
w:=initperm(n);
ret:=[initperm(n)];
#
for i from 2 to n! do
 v:=nextperm(w);
 if homtest(A,A,v) then 
  ret:=[op(ret),eval(v)]; 
 fi;
 w:=eval(v);
od;
#
eval(ret);
end;

qtest := proc(A)
#
# # Tests whether A is a quandle matrix in standard form
#
ret:=true;
#
for i from 1 to coldim(A) do
 if A[i,i] <> i then 
  ret:=false; 
  break; 
 fi;
od;
#
if ret then
 for i from 1 to coldim(A) do
  if not permtest(col(A,i)) then 
   ret:=false; 
   break; 
  fi;
 od;
fi;
# 
if ret then 
 ret:=q3test(A);
fi;
#
eval(ret);
end;


q3test := proc (A)
#
# # Tests a matrix for quandle axiom (iii)
#
v:= true;
#
for a from 1 to rowdim(A) do if v then
 for b from 1 to rowdim(A) do if v then
  for c from 1 to rowdim(A) do
   if A[A[a,b],c] <> A[A[a,c],A[b,c]]
    then v:=false;
    break;
   fi;
  od; fi;
 od; fi;
od;
v;
end;


permat := proc (A,B)
#
# # Used in stndiso
# #
# # applies permutation B to the entries of A
#
C:=array(1..coldim(matrix([B])),1..coldim(matrix([B])));
#
for i from 1 to rowdim(A) do
 for j from 1 to rowdim(A) do
  C[i,j]:=B[A[i,j]];
 od;
od;
eval(C);
end;


pm := proc (B)
#
# # Used in stndiso
# # 
# # computes permutation matrix of permutation B
#
C:=array(1..coldim(matrix([B])),1..coldim(matrix([B])));
#
for i from 1 to coldim(matrix([B])) do
 for j from 1 to coldim(matrix([B])) do
  C[i,j]:=0;
 od;
od;
#
for i from 1 to coldim(matrix([B])) do
 for j from 1 to coldim(matrix([B])) do
  if j=B[i] then 
   C[i,j]:=1; 
  fi;
 od;
od;
eval(C);
end;


stndiso := proc (A,v)
#
# # applies permutation v to A, then converts to standard form
#
multiply(inverse(pm(v)),multiply(permat(A,v),pm(v)));
end;


compare := proc (A,B)
#
# # compares two matrices
#
ret:=true;
for i from 1 to rowdim(A) do if ret then
 for j from 1 to rowdim(A) do
  if A[i,j] <> B[i,j] then 
   ret:=false; 
   break; 
  fi;
 od; fi;
od;
#
ret;
end;


isotest := proc (A,B)
#
# # determines whether quandles A and B are isomorphic
#
ret:=false;
R:=permute(coldim(A));
#
if coldim(A)=coldim(B) then
 if coldim(A)=1 then 
  ret:=true; 
 else
  for i from 1 to coldim(B)! do
   if homtest(A,B,R[i]) then 
    ret:=true; 
    break; 
   fi;
  od;
 fi;
fi;
#
ret;
end;


listpermsi := proc(n,i)
#
# # Used in quandleslist
# # 
# # lists all permutations of n letters with an i in the
# # ith position
#
G:=[];
R:=permute(n);
#
for j from 1 to coldim(matrix([R])) do 
 if R[j][i]=i then 
  G:=[op(G),eval(R[j])]; 
 fi; 
od;
#
eval(G);
end;


quandleslist := proc(n)
#
# # lists all quandles of order n
# # 
# # not very useful for n>5, unless your machine is very fast
# #
# # quandleslist5 below seems to work better on machines
# # with less memory
# 
R:=[];
clist:=[];
#
for i from 1 to n do
 clist:=[op(clist),eval(listpermsi(n,i))];
od;
#
l:=initmap(n);
w:=[];
for k from 1 to ((n-1)!)^n-1 do
 for m from 1 to n do
  w:=[op(w),clist[m][l[m]]];
 od; 
 M:=transpose(matrix(w));
  if q3test(M) then
   R:=[op(R),eval(M)];
  fi;
 l:=nextmap(l,(n-1)!);
 w:=[];
od;
#
eval(R);
end;

quandleslist5 := proc()
#
# # alternate algorithm for quandles of order 5
# # 
# # uses hard-coded nested loops in place of control vector
#
R:=[];
C1:=listpermsi(5,1);
C2:=listpermsi(5,2);
C3:=listpermsi(5,3);
C4:=listpermsi(5,4);
C5:=listpermsi(5,5);
#
for a in C1 do
 for b in C2 do
  for c in C3 do
   for e in C4 do
    for f in C5 do
      M:=transpose(matrix([a,b,c,e,f]));
      if q3test(M) then 
       R:=[op(R),eval(M)]; 
      fi;
    od;
   od;
  od;
 od;
od;
#
eval(R);
end;

reducelist := proc(R)
#
# # removes isomorphic quandles, giving a list with a
# # a single representative from each isomorphism class
#
r:=R;
n:= coldim(R[1]);
m:=coldim(matrix([R]));
#
for i from 2 to m do
 for j from 1 to i do
  if type(r[j],matrix) and 
     type(r[i],matrix) and 
     i<>j and 
     isotest2(r[i],r[j]) 
  then r[i]:=0; fi;
 od;
od;
#
F:=[];
for i from 1 to m do 
 if type(r[i],matrix) then 
  F:=[op(F),eval(r[i])]; 
 fi;
od;
#
eval(F);
end;


mplus:= proc(A,m)
#
# # used in cprod
# # 
# # adds m to every entry of A
#
B:=matrix(rowdim(A),coldim(A));
#
for i from 1 to rowdim(A) do
 for j from 1 to coldim(A) do
  B[i,j]:=A[i,j]+m;
 od;
od;
#
eval(B);
end;


cprod := proc(A,B)
#
# # finds a matrix for the matrix of the cartesian product of 
# # two quandle matrices (or other binary operation tables)
#
tmp:=[];
#
for i from 1 to rowdim(B) do
 for j from 1 to coldim(B) do
  tmp:=[op(tmp),mplus(A,coldim(A)*(B[i,j]-1))];
 od;
od;
#
ret:=blockmatrix(rowdim(B),coldim(B),tmp);
#
eval(ret);
end;


ractiontest:=proc(A,B,M)
#
# # tests whether M is a rack action of A on B
#
ret:=true;
#
if coldim(M)<> coldim(A) or rowdim(M)<> coldim(B) then ret:=false; 
else 
for j from 1 to coldim(M) do if ret then
 for k from 1 to coldim(M) do if ret then
  for i from 1 to rowdim(M) do
   if M[M[i,j],k]<>M[M[i,k],A[j,k]] then 
    ret:=false; 
    break; 
   fi;
  od; fi;
 od; fi;
od;
fi;
#
eval(ret);
end;


ractionlist := proc(A,B)
#
# # output a list of all rack actions of A on B by quandle 
# # automorphisms of B, M:A->aut(B)
#
R:=[];
autb:=autlist(B);
l:=initmap(coldim(A));
an:=coldim(matrix([autb]));
a:=coldim(A);
b:=coldim(B);
#
w:=[];
for k from 1 to (an)^a-1 do
 for m from 1 to coldim(A) do
  w:=[op(w),autb[l[m]]];
 od; 
 M:=transpose(matrix(w));
  if ractiontest(A,B,M) then
   R:=[op(R),eval(M)];
  fi;
 l:=nextmap(l,an);
 w:=[];
od;
#
eval(R);
end;


compattest:=proc(A,B,M,N)
#
# # test rack actions M:B->aut(A) and N:A->aut(B) for compatibility
#
Q:=blockmatrix(2,2,[A,M,mplus(N,coldim(A)),mplus(B,coldim(A))]);
#
ret:=true;
for i from 1 to coldim(A) do if ret then
 for j from 1+coldim(B) to coldim(A)+coldim(B) do if ret then
  for k from 1 to coldim(A) do
   if Q[Q[i,j],k]<>Q[Q[i,k],Q[j,k]] then 
    ret:=false; 
    break; 
   fi;
  od; fi;
 od; fi;
od;
for i from coldim(B)+1 to coldim(A)+coldim(B) do if ret then
 for j from 1 to coldim(A) do if ret then
  for k from coldim(B)+1 to coldim(A)+coldim(B) do
   if Q[Q[i,j],k]<>Q[Q[i,k],Q[j,k]] then 
    ret:=false; break; 
   fi;
  od; fi;
 od; fi;
od;
#
eval(ret);
end;


qunion:=proc(A,B)
#
# # find all quandle structures on A union B
#
rab:=ractionlist(A,B);
rba:=ractionlist(B,A);
R:=[];
#
for i from 1 to coldim(matrix([rab])) do
 for j from 1 to coldim(matrix([rba])) do
  if compattest(A,B,rba[j],rab[i]) then
   R:=[op(R),blockmatrix(2,2,
   [A,rba[j],mplus(rab[i],coldim(A)),mplus(B,coldim(A))])]; 
  fi;
 od;
od;
#
eval(R);
end;


subquandletest := proc(A,v)
#
# # tests whether a subset v of 1..n is a subquandle of A
#
ret:=true;
#
for i in v do if ret then
 for j in v do
  if not member(A[i,j],v) then 
   ret:=false; 
   break; 
  fi;
 od; fi;
od; 
#
eval(ret);
end;


qcomptest := proc(A,v)
#
# # determines whether a subset v of 1..n is an A-complemented
# # subquandle
#
S:={};
for i from 1 to coldim(A) do S:={op(S),i}; od;
C:=S minus eval(v);
ret:=false;
if subquandletest(A,v) and subquandletest(A,C) then 
 ret:=true; 
fi;
#
eval(ret);
end;


orbdec := proc(A)
#
# # computes the orbit decomposition of quandle with matrix A
# # given as subsets
#
ret:={};
S:={};
for i from 1 to coldim(A) do S:={op(S),i}; od;
C:=combinat[powerset](S);
Cred:={};
for v in C do 
 if qcomptest(A,eval(v)) then 
  Cred:={op(Cred),v}; 
 fi; 
od;
#
for i in S do 
 u:=S;
 for w in Cred do 
  if member(i,w) then
   u:=u intersect w;
  fi;
 od;
 ret:={op(ret),u};
od;
#
eval(ret);
end;


orbdecmat := proc(A)
#
# # returns the matrices of the orbit subquandles of A
#
ret:=[];
L:=convert(orbdec(A),list);
for i from 1 to coldim([L]) do
 M:=evalm(A);
 for j from 1 to coldim(A) do
  if not member(j,L[i]) then 
   for k from 1 to coldim(A) do
    M[k,j]:=0;
    M[j,k]:=0;
   od; 
  fi;
 od;
  for k from 1 to coldim(A) do
   if M[coldim(A)-k+1,coldim(A)-k+1] =0 then 
    M:=minor(eval(M),coldim(A)-k+1,coldim(A)-k+1);
   fi;
  od;
  ret:=[op(ret),evalm(M)];
od;
#
eval(ret);
end;


convert2stnd := proc(A)
#
# # convert quandle matrix to standard form
# # useful for subquandles
#
M:=eval(A);
N:=eval(A);
max:=1;
for i from 1 to coldim(A) do
 if A[1,i] >max 
  then max:=A[1,i];
 fi;
od;
for i from 1 to coldim(A) do
 for j from 1 to coldim(A) do
  M[i,j]:=A[i,j]+max+1;
 od;
od;
for i from 1 to coldim(A) do
 v[i]:=M[i,i]; 
od;
for i from 1 to coldim(A) do
 for j from 1 to coldim(A) do
  for k from 1 to coldim(A) do
   if M[i,j] = v[k] then N[i,j]:=k; fi;
  od;
 od;
od;
eval(N);
end;


isotest2 := proc(A,B)
#
# # Uses orbit decomposition to speed up testing for 
# # isomorphism
#
ret:=true;
if coldim(matrix([convert(orbdec(A),list)])) <> 
   coldim(matrix([convert(orbdec(B),list)])) then
 ret:=false;
fi;
if ret then 
 PA:=orbdecmat(A);
 PB:=orbdecmat(B);
 for i from 1 to coldim(matrix([PA])) do
  for j from 1 to coldim(matrix([PB])) do
   if type(evalm(PA[i]),matrix) and type (evalm(PB[j]),matrix)
      and isotest(convert2stnd(evalm(PA[i])),convert2stnd(evalm(PB[j]))) then
      PA[i]:=0; 
      PB[j]:=0;
   fi;
  od;
 od;
 for i from 1 to coldim(matrix([PA])) do
  if type(PA[i],matrix) then 
   ret:=false; 
  fi; 
 od;
fi;
if ret then 
 ret:=isotest(A,B); 
fi;
#
eval(ret);
end;

qdual := proc(A)
#
# # finds the matrix of the dual quandle of A
#
r:=[]; 
for i from 1 to coldim(A) do 
 r:=[op(r),0]; 
od;
c:=[];
for i from 1 to coldim(A) do
 c:=[op(c),r];
od;
B:=matrix(c);
if qtest(A) then
 for j from 1 to coldim(A) do
  for i from 1 to coldim(A) do
   B[A[i,j],j]:=i;
  od;
 od;
 eval(B);
else print("Error, A must be a quandle matrix") 
fi;
end;

homlist2 := proc(A,B)
#
# # lists homomorphisms from A to B
#
n:=coldim(A);
m:=coldim(B);
out:=[];
st:=[];
for i from 1 to n do
 st:=[op(st),0];
od;
w:=[eval(st)];
c:=true;
while c do
 c:=false;
 for j from 1 to coldim(matrix([w])) do
  v:=eval(w[1]);
  w:=subsop(1=NULL,w);
  l:=0;
  for i from 1 to n do
   if v[i]=0 then
    l:=i;
    break;
   fi;
  od;
  x:=eval(v);
  if l=0 then
   out:=[op(out),eval(x)];
  else
   for i from 1 to m do
    x[l]:=i;
    z:=homfill(A,B,x);
    if type(z,list) then
     w:=[op(w),eval(x)];
     c:=true;
    fi;
   od;
  fi;
 od;
od;
eval(out);
end;

homfill:=proc(A,B,v)
#
# v:A -> B
#
ret:=true;
c:=true;
w:=eval(v);
while c do
 c:=false;
 for i from 1 to coldim(A) do
  for j from 1 to coldim(A) do
   if w[i]<>0 and w[j]<>0 and A[i,j]<>0 and
      w[A[i,j]]<> B[w[i],w[j]] then
    if w[A[i,j]]=0 then
     w[A[i,j]]:=eval(B[w[i],w[j]]);
     c:=true;
    else
     ret:=false;
    fi;
   fi;
  od;
 od;
od;
if ret then
 eval(w);
else
 eval(ret);
fi;
end;

reducelist2 := proc(R)
#
# # removes isomorphic quandles, giving a list with a
# # a single representative from each isomorphism class
#
r:=[eval(R[1])];
#
for i from 2 to coldim(matrix([R])) do
 c:=false;
 for j in r do
  if isotest2(R[i],j) then 
   c:=true;
  fi;   
 od;
 if not c then
  r:=[op(r),evalm(j)];
 fi;
od;
eval(r);
end;