printf("To see the procedure in this file execute: Info(); \n");

Info:=proc()

printf("This file loads the packages: combinat, combstruct, group and LinearAlgebra and contains the following additional procedures: \n\n");

printf("Aut(M) [also AutFast(M)] and Inn(M) which returns the automorphism group and inner automorphism group given quandle  matrix M.  \n\n");

printf("AreIsomorphic(M,N) returns true if quandle M is isomorphic to quandle N, and false if not.\n"); printf("  Excecute n:=quandle order and P:=permute(n) before using AreIsomorphic! \n\n");

printf("OrderSequence(n,G) gives the order sequence of a permgroup G of degree n. \n\n");

printf("ConvertDisjCyc(M) converts the columns of M to disjcyc notation. \n\n");

printf("IsQuandle(A) tells whether or not A is a quandle matrix. \n\n");

printf("CycleType(x) gives the cycle type of a permutation given in disjcyc form. \n\n");

printf("CycleType2(x,n) gives the partition that describes the cycle type of x as element of Sym(n). \n\n");

printf("IsAbelianQuandle(A) obvious.\n\n");

printf("IsLeftDistributive(A) \n\n");

printf("IsHomo(f,A) \n");

printf("IsHomo2(f,A,B) \n");

printf("PrintMatrix(A) prints the matrix A \n");

printf("IsCrossedSet(A) returns true if A is crossed set, false if not \n");

printf("IsLatin(A) \n\n");

printf("Inv3(A) is an invariant of the quandle A. If Inv3(A) <> Inv3(B) then the quandles A and B are not isomorphic.\n\n");

printf("ConvertToLoop(A,e) converts idempotent quasigroup to a loop with identity e. \n\n");

printf("PRODUCT(A,B) gives the matrix of the quandle which is a product of the quandles whose matrices are A and B.\n\n");

printf("IsKei(A) returns true if A is involutory, false otherwise.\n\n");

printf("IsConnected(A) returns true if A is connected, false otherwise.\n\n");

printf("Trans(A) returns the transvection group of A as a permgroup.\n\n");

printf("DisjUnion(A,B) returns the trivial disjoint union of quandles A and B.\n\n");

printf("IsAlexander(A) tells whether or not a *connected* quandle is Alexander.\n\n");

printf("IsFaithful(A) tells whether or not a quandle is faithful.\n\n");

printf("DualQuandle(A) returns the dual quandle of A.\n\n");

printf("AbelianGroupTypes(n) returns the set of all factorizations of n as a product of prime powers.\n\n");

printf("NumberOfColorings:=proc(A::Matrix,B::list), and NumberOfColoring1s:=proc(A::Matrix,B::list), the latter fixed one color so is faster.\n");

end proc:

interface(rtablesize=50):

with(combinat):
with(numtheory):
with(combstruct):
with(group):
with(LinearAlgebra):

DualQuandle:=proc(A)
local i,j,k,B,n;
n:=RowDimension(A);
B:=Matrix(n,n);
for i from 1 to n do
for j from 1 to n do
k:=A[i,j];
B[k,j]:=i;
od:
od:
B;
end proc:

AbelianGroupTypes:=proc(n)
local x,y,i,a,T,AG,p,S;
description "returns the set of all factorizations of n as a product of prime powers";
  x:=ifactors(n);
  y:=x[2];
  for i from 1 to nops(y) do
    a:=y[i][1];
    p:=combinat[partition](y[i][2]);
    S[i]:=map(t->[seq(a^j,j=t)],p);
  od;
  T:=combinat[cartprod]([seq(S[i],i=1..nops(y))]):
  AG[n]:={}:
  while not(T[finished]) do
    AG[n]:=AG[n] union {T[nextvalue]()};
  od:
  map(ListTools:-FlattenOnce,AG[n]);
 end proc:



IsAlexander:=proc(A)
description "this works only if A is known to be connected."; 
evalb(IsLatin(A) and IsAbelianQuandle(A)); 
end proc:

IsFaithful:=proc(A)
local n,s;
n:=RowDimension(A);
s:=convert(Transpose(A),listlist);
s:=convert(s,set);
evalb(nops(s)=n);
end proc:




DisjUnion:=proc(A,B)
local n,k,i,j,C;
n:=RowDimension(A):
k:=RowDimension(B):
C:=Matrix(n+k,n+k):
  for i from 1 to n do
  for j from 1 to n do
    C[i,j]:=A[i,j]:
  od:
  od:
for i from 1 to k do
for j from 1 to n do
    C[n+i,j]:=n+i;
od:
od:
 for i from 1 to n do
 for j from 1 to k do
  C[i,n+j]:=i;
 od:
 od:
for i from 1 to k do
for j from 1 to k do
 C[n+i,n+j]:=B[i,j]+n;
od:
od:
  C;
end proc:

Trans:=proc(A)
  local T,T2,T3,x,y;
  T:=convert(Transpose(A),listlist);
  T2:=map(t->convert(t,disjcyc),T);
  T3:={}:
  for x in T2 do
   for y in T2 do
    T3:=T3 union {mulperms(x,invperm(y))};
   od:
 od:
 permgroup(RowDimension(A),T3);
end proc:

IsLatin:=proc(A)
local i,j,n;
n:=RowDimension(A);
for i from 1 to n do
if nops({seq(A[i,j],j=1..n)})< n then return false; fi;
od:
true;
end proc:

IsConnected:=proc(A)
local G,Orb;
G:=Inn(A):
Orb:=orbit(G,1);
evalb(nops(Orb)=RowDimension(A))
end proc:


PrintMatrix:=proc(A)
local i,j,n,m;
n:=RowDimension(A);
m:=ColumnDimension(A);
for i from 1 to n do
 for j from 1 to m do
  printf("%2d ",A[i,j]);
 od;
 printf("\n");
od:
NULL;
end proc:

IsCrossedSet:=proc(A)
local n,a,b;
n:=RowDimension(A);
for a from 1 to n do
for b from 1 to n do
if evalb(A[a,b] = a) <> evalb(A[b,a]=b) then return false; fi;
od:
od:
true;
end proc:


IsAbelianQuandle:=proc(A)
local n,a,b,c,d;
n:=RowDimension(A);
for a from 1 to n do
for b from 1 to n do
for c from 1 to n do
for d from 1 to n do
 if A[A[a,b],A[c,d]]<>A[A[a,c],A[b,d]] then return false; fi;
od:
od:
od:
od:
true;
end proc:

IsLeftDistributive:=proc(M)
local i,j,k,n;
n:=RowDimension(M);
for i from 1 to n do
for j from 1 to n do
for k from 1 to n do
if M[k,M[i,j]]<>M[M[k,i],M[k,j]] then return false; fi:
od:
od:
od:
true:
end proc:


IsQuandle:=proc(M)
local i,j,k,n;
n:=RowDimension(M);
for i from 1 to n do
  if M[i,i]<>i then return false; fi;
od:
for i from 1 to n do
for j from 1 to n do
for k from 1 to n do
if M[M[i,j],k]<>M[M[i,k],M[j,k]] then return false; fi:
od:
od:
od:
for j from 1 to n do
if nops({seq(M[i,j],i=1..n)}) < n then return false; fi;
od:
true:
end proc:

CycleType2:=proc(x,n)
  local y,k,j;
  if x = [] then return [seq(1,i=1..n)]; fi;
  y:=[op(map(nops,x))];
  k:=add(j,j=y):
  sort([seq(1,i=1..n-k),op(y)]);
end proc:


OrderSequence:=proc(n,G)
local g,E;
description "if G is a permgroup of degree n return the sequence of orders of elements of G";
E:=elements(G);
sort([seq(grouporder(permgroup(n,{g})),g=E)]);
end proc:



ConvertDisjCyc:=proc(M)
local i,L;
description "this returns the list of columns in disjoint cycle notation of quandle matrix M";
Transpose(M);
convert(%,listlist);
map(convert,%,disjcyc);
end proc:

CycleType:=proc(x)
if x = [] then return 1; fi;
sort(map(nops,x));
end proc:

IsHomo:=proc(p,M)
local i,j,n;
description "IsHomo(p,M) = true if permutation p is homomorphism of quandle with matrix M, false if not";
 n:=LinearAlgebra:-RowDimension(M);
 for i from 1 to n do
    for j from 1 to n do
     if p[M[i,j]] <> M[p[i],p[j]]then return false; fi;
    od:
 od:
true;
end proc:

IsHomo2:=proc(f,B,A)
local i,j;
for i from 1 to RowDimension(B) do
for j from 1 to RowDimension(A) do
 if f(B[i,j])<>A[f(i),f(j)] then return false; fi;
od:
od:
true;
end proc:

Aut:=proc(M)
 local i,p,A,AA,allp,n;
 description "Aut(M) is the automorphism group of quandle with matrix M";
 n:=RowDimension(M);
 A:={}:
 allp := iterstructs(Permutation([seq(i,i=1..n)])):
 while `not`(finished(allp)) do
    p:=nextstruct(allp);
    if IsHomo(p,M) then A:=A union {p}; fi;
 od:
AA:=map(convert,A,disjcyc);
permgroup(n,AA);
end proc:

PartitionRowsAndCols:=proc(A::Matrix)
local n,x,y,i,S,j;
n:=RowDimension(A);
convert(A,listlist):
map(convert,%,set);
x:=map(nops,%);
y:=map(CycleType,ConvertDisjCyc(A));
for i from 1 to n do
  S[i]:={}:
  for j from 1 to n do
   if [x[i],y[i]]=[x[j],y[j]] then S[i]:=S[i] union {j}; fi:
  od:
od:
{seq(S[i],i=1..n)};
convert(%,list);
map(convert,%,list);
end proc:

AutFast:=proc(A)
local x,n,Parts,i,PP,GS,T,q,j,u,v,h,f;
n:=RowDimension(A);
Parts:=PartitionRowsAndCols(A);
#if nops(Parts)=1 then return FAIL; fi;
for i from 1 to nops(Parts) do
  PP[i]:=permute(Parts[i]);
od:
GS:={};
T := cartprod([seq(PP[t],t=1..nops(Parts))]):
while `not`(T[finished]) do
  q:=T[nextvalue]();
  for j from 1 to nops(Parts) do
    for u from 1 to nops(Parts[j]) do
      h[Parts[j][u]]:=q[j][u];
    od:
  od:
  f:=[seq(h[v],v=1..n)];
  if IsHomo(f,A) then GS:=GS union {convert(f,'disjcyc')}; fi;
end do:
permgroup(n,GS);
end proc:



Inn:=proc(M)
 local S,i,x,n;
 description "Inn(M) is the inner automorphism group of quandle with matrix M";
 n:=RowDimension(M);
 convert(Transpose(M),listlist);
 convert(%,set);
 map(convert,%,disjcyc);
 permgroup(n, %):
end proc:

n:=3:
P:=permute(n):
AreIsomorphic:=proc(A,B)
local i,j,q,flag;
global n,P;
for q in P do
  flag:=true;
  for i from 1 to n do
   for j from 1 to n do
    if q[A[i,j]]<>B[q[i],q[j]] then flag:=false; break; fi;
   od:
   if flag=false then break; fi;
  od:
  if flag=true then return true; fi;
od:
false;
end proc:

Floyd:=proc(f,x0)
local tortoise,hare,mu,lam;
tortoise,hare:=f(x0),f(f(x0));
while tortoise <> hare do
  tortoise:=f(tortoise);
  hare:=f(f(hare));
od:
mu:=0;
tortoise,hare:=x0,tortoise;
while tortoise<>hare do
  tortoise:=f(tortoise);
  hare:=f(hare);
  mu := mu+1;
od:
lam:=1;
hare:=f(tortoise);
while tortoise<>hare do
  hare:=f(hare);
  lam:=lam+1;
od;
return [mu,lam];
end proc:
FloydRowInvariant:=proc(A,a)
local f,x,S;
f:=x->A[a,x];
S:=[seq](Floyd(f,x),x=1..RowDimension(A));
convert(%,multiset);
convert(%,set);
end proc:

IndegreeRow:=proc(A,a)
local x,L;
L:=convert(Row(A,a),list);
convert(L,multiset);
sort(map(x->x[2],%));
end proc:

ColumnCycleType:=proc(A,a)
convert(Column(A,a),list);
convert(%,disjcyc);
map(nops,%);
sort(%);
end proc:


Inv3:=proc(A)
local i;
[seq]([IndegreeRow(A,i),FloydRowInvariant(A,i),ColumnCycleType(A,i)],i=1..RowDimension(A));
convert(%,multiset);
convert(%,set);
end proc:

ConvertToLoop:=proc(A,e)
local n,invR,invL,x,L,y;
description "this converts an idempotent quasigroup into a loop with identity e";
n:=RowDimension(A);
for x from 1 to n do
invR[e,A[x,e]]:=x;
invL[e,A[e,x]]:=x;
od:
L:=Matrix(n,n):
for x from 1 to n do
for y from 1 to n do
L[x,y]:=A[invR[e,x],invL[e,y]];
od:
od:
L;
end proc:

PRODUCT:=proc(A,B)
local i,f,T,L,inv,n,m,x;
description "this constructs the matrix of the product of the quandles whose matrices are A and B";
n:=LinearAlgebra:-RowDimension(A);
m:=LinearAlgebra:-RowDimension(B);
L:=[[seq(i,i=1..n)],[seq(j,j=1..m)]];
T:=combinat[cartprod](L):
i:=0:
while `not`(T[finished]) do
x:=T[nextvalue]();
i:=i+1;
f(i):=x;
inv(x):=i;
od:
Matrix(n*m,n*m,(i,j)->inv([A[f(i)[1],f(j)[1]],B[f(i)[2],f(j)[2]]]));
end proc:

IsKei:=proc(A)
local i,j,n;
n:=RowDimension(A);
for i from 1 to n do
for j from 1 to n do
if A[A[i,j],j]<>i then return false; fi;
od:
od:
true;
end proc:

Conj:=proc(Y,X)
description "Conj(Y,X) = X.Y.X^(-1) where X,Y are in S_n in disjcyc notation"; 
mulperms(invperm(X),mulperms(Y,X)); 
end proc:

ConjQuandle:=proc(S)
local M,i,j,n;
description "converts subquandle S of conj(S_n) to a Matrix";
n:=nops(S):
M:=Matrix(n,n,0);
for i from 1 to n do
f(S[i]):=i;
od:
for i from 1 to n do
for j from 1 to n do
  M[i,j]:=f(Conj(S[i],S[j]));
od:
od:
M;
end proc:

IsColoring:=proc(a::list,B::list,f::Matrix,g::Matrix)
local n,i,v,k,S,N,j;
  n,N:=max(op(map(abs,B)))+1,nops(B);
   for i to n do 
     v[0, i] := a[i];
   end do;
  #lprint(seq(v[0,i],i=1..n));
   for k to N do 
    i := B[k]; 
    for j to n do 
     if (j <> abs(i)) or (j <> abs(i)+1) then v[k, j] := v[k-1, j] end if;
    end do; 
    if i > 0 then 
       v[k, i], v[k, i+1] := v[k-1, i+1], f[v[k-1, i], v[k-1, i+1]];
    elif i < 0 then 
      i := -i; 
      v[k, i], v[k, i+1] := g[v[k-1, i+1], v[k-1, i]], v[k-1, i];
    end if;
   # lprint(seq(v[k,i],i=1..n));
  end do; 
  for i to n do 
    if v[0, i] <> v[N, i] then return false; fi;
  end do;
  true;
end proc:

NumberOfColorings:=proc(A::Matrix,B::list)
local A2,i,j,k,T,count,xx,q,b;
  q:=RowDimension(A);
  b:=max(op(map(abs,B)))+1;
  A2:=DualQuandle(A);
  T:=cartprod([seq([seq(i,i=1..q)],t=1..b)]):
  count:=0;
  while not T[finished] do
   xx:=T[nextvalue]();
   if  IsColoring(xx,B,A,A2) then count:=count+1; fi;
  od:
 return count;
end proc:

NumberOfColorings1:=proc(A::Matrix,B::list)
local A2,i,j,k,T,count,xx,q,b;
  q:=RowDimension(A);
  b:=max(op(map(abs,B)))+1;
  A2:=DualQuandle(A);
  T:=cartprod([seq([seq(i,i=1..q)],t=1..b-1)]):
  count:=0;
  while not T[finished] do
   xx:=[1,op(T[nextvalue]())];
   if  IsColoring(xx,B,A,A2) then count:=count+1; fi;
  od:
 return count*q;
end proc: