read("quandles-maple.txt");

assoctest := proc(A)
#
# # tests a Cayley matrix for associativity
#
ret:=true;
for i from 1 to coldim(A) do if ret then
 for j from 1 to coldim(A) do if ret then
  for k from 1 to coldim(A) do
   if A[A[i,j],k]<> A[i,A[j,k]] then 
    ret:=false; 
    break; 
   fi;
  od; fi;
 od; fi;
od;
#
eval(ret);
end;

commtest := proc(A)
#
# # tests a Cayley matrix for commutativity
#
ret:=true;
for i from 1 to coldim(A) do if ret then
 for j from 1 to coldim(A) do 
  if A[i,j]<>A[j,i] then 
   ret:=false; 
   break; 
  fi;
 od; fi;
od;
#
eval(ret);
end;

invtest := proc(A)
#
# # test possibly incomplete Cayley matrix for presence of 
# # row or column without inverses. assumes standard from
#
R:=true;
#
for i from 1 to coldim(A) do
 RT:=false;
 RZ:=false;
 CT:=false;
 CZ:=false;
 for j from 1 to coldim(A) do
  if A[i,j]=1 or A[i,j]=0 then RT:=true; fi;
  if A[j,i]=1 or A[j,i]=0 then CT:=true; fi;
 od;
 if not CT or not RT then R:=false; fi;
od;
#
eval(R);
end;

abqtest :=proc(A)
#
# test whether a quandle is abelian
#
n:=coldim(A);
C:=false;
for i from 1 to n do
 if not C then
  for j from 1 to n do
   if not C then
    for k from 1 to n do
     if not C then
      for l from 1 to n do
       if A[A[i,j],A[k,l]]<>A[A[i,k],A[j,l]] then
        C:=true;
       fi;
      od;
     fi;
    od;
   fi;
  od;
 fi;
od;
eval(not C)
end;

cafill := proc(A)
#
# # uses commutativity and associativity to fill in zeros
# # in a Cayley matrix
#
C:=false;
ch:=true;
R:=evalm(A);
while ch do
 ch:=false;
 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 R[i,j]<>0 and R[j,k]<>0 then
     if R[R[i,j],k] <> R[i,R[j,k]] and
        R[R[i,j],k] <> 0 and R[i,R[j,k]]<>0 then
      C:=true; 
     elif R[R[i,j],k]=0 and R[i,R[j,k]]<>0 then
      R[R[i,j],k]:=eval(R[i,R[j,k]]); 
      ch:=true;
     elif R[i,R[j,k]]=0 and R[R[i,j],k]<>0 then
      R[i,R[j,k]]:=eval(R[R[i,j],k]); 
      ch:=true;
     fi;
    fi;
   od;
  od;
 od;
 for i from 1 to coldim(A) do
  for j from 1 to coldim(A) do   
   if R[i,j]<>0 and R[j,i]<>0 and R[i,j]<>R[j,i]
    then C:=true;
   fi;
   if R[i,j]=0 and R[j,i]<>0 then
      R[i,j]:=eval(R[j,i]);
    ch:=true;
   fi;    
   if R[j,i]=0 and R[i,j]<>0 then
      R[j,i]:=eval(R[i,j]);
    ch:=true;
   fi;    
  od;
 od;
od;
if C then 
 eval(not C);
else
 evalm(R);
fi;
end;

findzero :=proc(A)
#
# # finds first zero position
#
ret:=false;
C:=true;
for i from 1 to coldim(A) do 
 if C then
  for j from 1 to coldim(A) do
   if A[i,j]=0 then
    ret:=[i,j];
    C:=false;
    break;
   fi;
  od; fi;
od;
eval(ret);
end;

zerofill := proc(A)
#
# # returns all standard-form Cayley matrices
# # matching the given partly-filled in matrix
#
cont:=true;
ret:=[evalm(A)];
while cont do
 cont:=false;
 for j from 1 to coldim(matrix([ret])) do
  X:=evalm(ret[1]);
  f:=findzero(X);
  ret:=subsop(1=NULL,ret);
  if type(f,list) and invtest(X) then
   B:=evalm(X);
   for i from 1 to coldim(A) do
    B[f[1],f[2]]:=i;
    Z:=cafill(B);
    if type(Z,matrix) and invtest(Z) then
     ret:=[op(ret),evalm(Z)];
     if type(findzero(Z),list) then
      cont:=true;
     fi;
    fi;
   od;
  else 
   ret:=[op(ret),evalm(X)];
  fi;
 od;
od;
eval(ret);
end;

abgroupfill := proc(Q)
#
# # fills in standard form Cayley matrix from a quandle 
# # matrix 
#
R:=matrix(coldim(Q),coldim(Q));
for i from 1 to coldim(Q) do
 R[i,1]:=i;
 R[1,i]:=i;
 for j from 2 to coldim(Q) do
  R[i,j]:=0;
 od;
od;
#
C:=false;  # contradiction counter
ch:=true;  # changed counter 
#
while ch do
 ch:=false;
 for i from 1 to coldim(Q) do 
  if not C then
   for j from 1 to coldim(Q) do
    if not C then 
     for k from 1 to coldim(Q) do
      if R[Q[i,j],Q[j,k]]<>0 and R[Q[i,k],j]<>0 and
         R[Q[i,j],Q[j,k]]<>R[Q[i,k],j] then
        C:=true;
        break;
      fi;
      if R[Q[i,j],Q[j,k]]=0 and R[Q[i,k],j]<>0 then
       R[Q[i,j],Q[j,k]]:=eval(R[Q[i,k],j]);
       ch:=true;
      elif R[Q[i,j],Q[j,k]]<>0 and R[Q[i,k],j]=0 then
       R[Q[i,k],j]:=eval(R[Q[i,j],Q[j,k]]);
       ch:=true;
      fi;
     od;
    fi;
   od;
  fi;
 od;
 for i from 1 to coldim(Q) do
  if not C then
   for j from 1 to coldim(Q) do
    if R[Q[i,j],Q[j,i]]<>0 and R[i,j]<>0 and
       R[Q[i,j],Q[j,i]]<>R[i,j]  then
     C:=true;
     break;
    fi;
    if R[Q[i,j],Q[j,i]]=0 and R[i,j]<>0 then
     R[Q[i,j],Q[j,i]]:=eval(Q[i,j]);
     ch:=true;
    elif R[Q[i,j],Q[j,i]]<>0 and R[i,j]=0 then
     R[i,j]:=eval(R[Q[i,j],Q[j,i]]);
     ch:=true;
    fi;
   od;
  fi;
 od;
 T:=cafill(R);
 if type(T,matrix) then
  R:=evalm(T);
 else 
  C:=true;
 fi;
od;
if not C then
 eval(R);
else
 eval(not C);
fi;
end;

alextest:=proc(Q)
#
# # determines whether a quandle is Alexander
#
if abqtest(Q) then
 A:=abgroupfill(Q);
 ret:=[];
 if type(A,matrix) then
  X:=zerofill(A);
  for Z in X do
   if assoctest(Z) and commtest(Z) and 
      invtest(Z) and homtest(Z,Z,eval(col(Q,1))) then
    ret:=[op(ret),evalm(Z)];
   fi;
  od;
 fi;
 if coldim(matrix([ret])) =0 then 
  ret:=false;
 fi;
 else ret:=false;
fi;
eval(ret);
end;

alexquandle := proc(Q,v);
#
# # computes the quandle matrix for the Alexander quandle 
# # with Cayley matrix Q and automorphism v
#
if homtest(Q,Q,v) then
 tmp1:=[]; 
 for i from 1 to coldim(Q) do 
  tmp1:=[op(tmp1),0]; 
 od;
 n:=eval(tmp1);
 tmp2:=[]; 
 for i from 1 to coldim(Q) do 
  tmp2:=[op(tmp2),tmp1]; 
 od;
 M:=matrix(tmp2);
 for i from 1 to coldim(Q) do
  for j from 1 to coldim(Q) do
   if Q[i,j]=1 then  
    n[i]:=j; 
   fi;
  od;
 od;
#
 for i from 1 to coldim(Q) do
  for j from 1 to coldim(Q) do
   M[i,j]:=Q[v[i],Q[j,n[v[j]]]];
  od;
 od;
 evalm(M);
else 
 print("Error, second argument must be an automorphism"); 
fi;
end;

cayley := proc (n)
#
# # computes the cayley matrix of Z mod n
#
tmp1:=[]; 
for i from 1 to n do 
 tmp1:=[op(tmp1),0]; 
od;
tmp2:=[]; 
for i from 1 to n do 
 tmp2:=[op(tmp2),tmp1]; 
od;
M:=matrix(tmp2);
for i from 1 to n do
 for j from 1 to n do
  M[i,j]:= i+j-2 mod n +1;
 od;
od;
eval(M);
end;

conjq := proc(G)
#
# # Returns the conjugation quandle of the group with Cayley
# # matrix G
#
tmp1:=[]; 
for i from 1 to coldim(G) do 
 tmp1:=[op(tmp1),0]; 
od;
n:=eval(tmp1);
tmp2:=[]; 
for i from 1 to coldim(G) do 
 tmp2:=[op(tmp2),tmp1]; 
od;
M:=matrix(tmp2);
for i from 1 to coldim(G) do
 for j from 1 to coldim(G) do
  if G[i,j]=1 then  
   n[i]:=j; fi;
 od;
od;
for i from 1 to coldim(G) do
 for j from 1 to coldim(G) do
  M[i,j]:= G[n[j],G[i,j]];
 od;
od;
eval(M);
end;