Print("Loaded the following functions: IsFaithful(Q), ConjugationQuandleN_fold(C,n),ConjugationQuandle(C),IsQuandle(Q), Type(Q),IsConjugationQuandle(Q),FiniteEnvelopingGroup(Q), IdConjugationQuandle(Q), Inn(Q), IsConnected(Q), IsMedialQuandle(Q), IsLatinQuandle(Q), IsKei(Q), DualQuandle(Q), IsColoring(a,B,f,g), NumberOfColorings(A,Braid), NumberOfColorings3(A,Braid),IsSimpleQuandle(Q),IsCrossedSet(Q),  TripleData(Q1),NumberOfColoringsT(A,PT, Braid ), multiset(L),Orbits2(Q),Orbits2_1(Q),Orbits3(Q,orbs2), Orbits4(Q,orbs3), MakeGAlex(f,G),MakeGAlex2(x,G) (where <x^G> = G), MakeGAlex3(Q) (Q is GAlex),BraidFunction(a,B,f), Psi(B,Q) = vector invariant, ConnectedSum(B1,B2),IsGeneralizedAlexander(Q), ImagePhi(Q), Fix(G,f) \n");
#note the GAP command  Collected works for multiset
#also it has some sorting properties.

RTrans:=function ( Q )
    local  B, e, n, m, i, j;
    B := Set( TransposedMat( Q ) );
    ;
    e := List( B, function ( x )
            return PermList( x );
        end );
    return e;
end;;

Type:=function(Q)
return Order( RTrans(Q)[1]);
end;;

FiniteEnvelopingGroup:=function(Q)
local i,j,f,RelList,g,H,n,t;
n:=Size(Q);
f:=FreeGroup(n);
RelList:=[];
for i in [1..n] do
  for j in [1..n] do
   Add(RelList,f.(j)*f.(Q[i][j])*f.(j)^(-1)*f.(i)^(-1));
  od;
od;
for t in [1..n] do Add(RelList,f.(t)^Type(Q)); od;
g:=f/RelList;
return g;
end;;

#If Q is a conjugation quandle then this procedure returns the finite
#enveloping group (for which Q is a conjugation quandle), otherwise
#it returns fail. 


IsConjugationQuandle:=function(Q)
local G,U;
G:=FiniteEnvelopingGroup(Q);
if (Size(G) in [1024, 1152, 1536,  1920]) or Size(G) > 2000 then U:=Size(G); else U:=IdGroup(G); fi;
return [Size(Q) = Size(Set(GeneratorsOfGroup(G))),U ];
end;;


multiset:=function(L)
local S,MS,i,j,x;
S:=Set(L);
MS:=[];
for j in S do
 x:=0;
 for i in L do
  if i = j then x:=x+1; fi;
 od;
 Add(MS,[j,x]);
od;
return MS;
end;; 

ImagePhi:=function(A)
local B,e,n,m,i,j;
B:=Set(TransposedMat(A));;
e:=List(B,x->PermList(x));
n:=Length(e);
m := NullMat(n,n);
for i in [1..n] do
  for j in [1..n] do
     m[i][j]:=Position(e,e[i]^e[j]);
  od;
od;
return m;
end;;

#Fix(G,f) for f in Aut(G) gives the subgroup of x in G such that f(x) = x.

Fix:=function(G,f)
local x,H;
H:=[];
for x in Elements(G) do
 if Image(f,x) = x then Add(H,x); fi;
od;
return Subgroup(G, H);;
end;;



IsCrossedSet:=function(A)
local n,a,b;
n:=Length(A);
for a in [1..n] do
for b in [1..n] do
if (A[a][b] = a) <> (A[b][a]=b) then return false; fi;
od;
od;
return true;
end;;


Inn:=function(A)
local B,LL;
B:=TransposedMat(A);
   LL:=List(B,x->PermList(x));
   return Group(LL);
end;;

IsGeneralizedAlexander:=function(Q)
local g,h;
g:=Inn(Q);
h:=DerivedSubgroup(g);
return Size(h) = Size(Q);
end;;

IsConnected:=function(A)
return IsTransitive(Inn(A),[1..Length(A)]);
end;; 

IsFaithful:=function(A)
local B,LL;
B:=TransposedMat(A);
return Size(B) = Size(Set(B));
end;;

IdConjugationQuandle:=function(QM)
local g;
g:=Inn(QM);
return [IdSmallGroup(g), StructureDescription(g)];
end;;

ConjugationQuandleN_fold:=function(C, n)
  local m,i,j,e;
  e := Elements(C);
  m := NullMat(Size(e),Size(e));
  for i in [1..Size(e)] do
    for j in [1..Size(e)] do
      m[i][j] := Position(e, e[i]^(e[j]^n));
    od;
  od;
  return m;
end;;

ConjugationQuandle:=function(C)
local e,i,j,n,m;
e:=Elements(C);
n:=Length(e);
m := NullMat(n,n);
for i in [1..n] do
for j in [1..n] do
 m[i][j]:=Position(e,e[i]^e[j]);
od;
od;
return m;
end;;

IsQuandle:=function(x)
local n,i,j,k,y;
n:=Length(x);
for i in [1..n] do 
  if x[i][i]<>i then return false; fi;
od;
y:=TransposedMat(x);
for i in [1..n] do
  if PermList(y[i]) = fail then return false; fi;
od;
for i in [1..n] do
for j in [1..n] do
for k in [1..n] do
 if x[x[i][j]][k] <> x[x[i][k]][x[j][k]] then return false; fi;
od;
od;
od;
return true;
end;;

IsMedialQuandle:=function(A)
local n,a,b,c,d;
n:=Length(A);
for a in [1..n] do
for b in [1..n] do
for c in [1..n] do
for d in [1..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;
return true;
end;;

IsLatinQuandle:=function(A) 
local r,n; 
n := Length(A); 
for r in A do 
if Length(Set(r)) < n then return false; fi;
od;
return true;
end;;

IsKei:=function(A) 
local i, j, n; 
n := Length(A); 
for i in [1..n] do
  for j in [1..n] do
    if A[A[i][ j]][ j] <> i then return false; fi;
 od;
od;
return true;
end;;

DualQuandle:=function(A)
local i,j,k,B,n;
n:=Length(A);
B:=NullMat(n,n);
for i in [1..n] do
for j in [1..n] do
k:=A[i][j];
B[k][j]:=i;
od;
od;
return B;
end;;

IsColoring:=function(a,B,f,g)
local n,i,v,k,S,N,j;
  n:=Maximum(List(B,AbsoluteValue))+1;
  N:=Length(B)+1;
  v:=NullMat(N,n);
   for i in [1..n] do 
     v[1][i] := a[i];
   od;
   for k in [2..N] do 
    i := B[k-1]; 
    for j in [1..n] do 
     if (j <> AbsoluteValue(i)) or (j <> AbsoluteValue(i)+1) then v[k][ j] := v[k-1][ j]; fi;
    od; 
    if i > 0 then 
       v[k][ i]:=  v[k-1][ i+1];
       v[k][ i+1] := f[v[k-1][ i]][ v[k-1][ i+1]];
    elif i < 0 then 
      i := -i; 
      v[k][ i]:= g[v[k-1][ i+1]][ v[k-1][ i]];
      v[k][ i+1] := v[k-1][ i];    
    fi; 
  od; 
  for i in [1..n] do 
    if v[1][i] <> v[N][i] then return false; fi;
  od;
 return true;
 end;;

##in the following A is a quandle matrix and Braid is the
##braid description of a knot as a sequence of integers

NumberOfColorings:=function(A,Braid)
local A2,i,enum,xx,yy,q,b,count;
  q:=Length(A);
  b:=Maximum(List(Braid,AbsoluteValue))+1;
  A2:=DualQuandle(A);
  enum:=EnumeratorOfTuples([1..q],b-1);
  count:=0;
  i:=0;
  while i < q^(b-1) do
   i:=i+1;
   yy:=enum[i];
   xx:=ShallowCopy(yy);
   Add(xx,1);
   if  IsColoring(xx,Braid,A,A2) then count:=count+1; fi;
  od;
 # return count*q; #for all colorings
  return count*q - q; #for non-trivial colorings
end;;

SetOfColorings:=function(A,Braid)
local A2,i,enum,xx,yy,q,b,SS;
  SS:=[];
  q:=Length(A);
  b:=Maximum(List(Braid,AbsoluteValue))+1;
  A2:=DualQuandle(A);
  enum:=EnumeratorOfTuples([1..q],b);
  i:=0;
  while i < q^b do
   i:=i+1;
   yy:=enum[i];
   xx:=ShallowCopy(yy);
   if  IsColoring(xx,Braid,A,A2) then Add(SS,xx); fi;
  od;
  return SS;
end;;


##The following is for braid index 3 knots only
NumberOfColorings3:=function( A, Braid )
    local a1,a2, A2, i, enum, xx,  q, b, count;
    q := Length( A );
    b := Maximum( List( Braid, AbsoluteValue ) ) + 1;
    A2 := DualQuandle( A );   
    count := 0;
    for a1 in [1..q] do
     for a2 in [1..q] do
       xx:=[1,a1,a2];       
        if IsColoring( xx, Braid, A, A2 )  then
            count := count + 1;
        fi;
    od;
   od;
    return count * q - q;
end;;

IsSimpleQuandle:=function(q)
local g,N,gg,n;
if IsFaithful(q) = false then return false; fi;
#if IsConnected2(q) = false then return false; fi;
g:=Inn(q);;
if Size(Center(g))>1 then return false; fi;
N:=NormalSubgroups(g);;
gg:=DerivedSubgroup(g);;
for n in N do
  if Size(n) = 1 then continue; fi;
  if IsSubset(gg,n) and Size(n)<Size(gg) then return false;  fi;
od;
return true; 
end;;

#TripleData(Q1) list of [[i1,i2,i3],N([i1,i2,i3])] 
#where [i1,i2,i3] is representative of a triple orbit of 
#Inn(Q1) and N([i1,i2,i3]) is the size of this orbit.

TripleData:=function(Q1)
  local z,g,n,orbs1,orbs2,cn,x,y,h,S,TD;
  TD:=[];
  g:=Inn(Q1);
  n:=Size(Q1);
  orbs2:=OrbitsDomain(g,Tuples([1..Size(Q1)],2),OnTuples); 
  cn:=0;
  for x in orbs2 do
   y:=x[1];
   h:=Stabilizer(g,y,OnPairs);
   S:=[1..Size(Q1)];
    orbs1:=Orbits(h,S);
    for z in orbs1 do
     Add( TD, [[y[1],y[2],z[1]], Size(x)*Size(z)]);
     cn:=cn+Size(x)*Size(z);
    od;
  od;
  if cn <> n^3 then Print("Error! \n"); fi;
  return TD;
end;; 

#in the following A = quandle matrix
#PT = output of TripleData(A)
#Braid = braid description of knot

NumberOfColoringsT:=function (A,PT, Braid )
    local  A2, i,cc,j,jj,zz,pt, enum, xx, yy, q, b, count;
    q := Length( A );
    b := Maximum( List( Braid, AbsoluteValue ) ) + 1;
    A2 := DualQuandle( A );
    cc:=[];
    for j in [1..Size(PT)] do cc[j]:=0; od;
    enum := EnumeratorOfTuples([ 1 .. q ], b - 3 );    
      i := 0;
      while i < q ^ (b - 3) do
        i := i + 1;
        yy := enum[i];
        #xx := ShallowCopy( yy );    
        for j in [1..Size(PT)] do
          pt:=PT[j];    
          zz:=Concatenation(pt[1], yy);
          if IsColoring( zz, Braid, A, A2 )  then
             cc[j]:=cc[j]+1;
          fi;
       od;
    od;
   count:=0;
   for jj in [1..Size(PT)] do
      count:=count+cc[jj]*PT[jj][2];
   od;
   return count-q;
end;;

Orbits2:=function(Q)
local n,orbs,g,h,y, OrbsData;
 n:=Size(Q);
 g:=Inn(Q);
 OrbsData:=[];
 h:=Stabilizer(g,1);
 orbs:=Orbits(h,[1..n]);
 for y in orbs do
  Add(OrbsData,[[1,y[1]],n*Size(y)]);
 od;
 return OrbsData;
end;; 

Orbits2_1:=function(Q)
local n,orbs,g,h,y, OrbsData;
 n:=Size(Q);
 g:=Inn(Q);
 OrbsData:=[];
 h:=Stabilizer(g,1);
 orbs:=Orbits(h,[1..n]);
 for y in orbs do
  Add(OrbsData,[1,y[1],n*Size(y)]);
 od;
 return OrbsData;
end;; 


Orbits3:=function(Q,orbs2 )
    local  z, g, n,sorbs, x, y, h, orbs3;
    orbs3 := [  ];
    g := Inn( Q );
    n := Size( Q );
    for x  in orbs2  do
        y := x[1];
        h := Stabilizer( g, y, OnPairs );
        sorbs := Orbits(h, [1..n]);
        for z  in sorbs  do
            Add( orbs3, [ [ y[1], y[2], z[1] ], x[2] * Size(z) ] );
        od;
    od;
    return orbs3;
end;;

Orbits4:=function(Q,orbs3 )
    local  z, g, n,sorbs, x, y, h, orbs4;
    orbs4 := [  ];
    g := Inn( Q );
    n := Size( Q );
    for x  in orbs3  do
        y := x[1];
        h := Stabilizer( g, y, OnTuples );
        sorbs := Orbits(h, [1..n]);
        for z  in sorbs  do
            Add( orbs4, [ [ y[1], y[2],y[3], z[1] ], x[2] * Size(z) ] );
        od;
    od;
    return orbs4;
end;;

MakeGAlex:=function(f,G)
local e,n,QM,i,j;
  e:=Elements(G);
  n:=Length(e);
  QM:=List([1..n],t->[1..n]);
    for i in [1..n] do
      for j in [1..n] do
       QM[i][j]:=Position(e,Image(f,e[i]*e[j]^(-1))*e[j]);
      od;
    od;
  return QM;
end;;

ConnectedSum:=function(B1,B2)
local t,B3,m;
m:=Maximum(List(B1,AbsoluteValue)) + 1;
return Concatenation(B1,List(B2,t->t + SignInt(t)*(m-1)));
end;

#in following need <x^G> = G.

MakeGAlex2:=function ( x, G )
    local  e, n, QM, i, j,f;
    f:=InnerAutomorphism(G,x);
    e := Elements(DerivedSubgroup(G) );
    n := Length( e );
    QM := NullMat(n,n);
    for i  in [ 1 .. n ]  do
        for j  in [ 1 .. n ]  do
            QM[i][j] := Position( e, Image( f, e[i] * e[j] ^ -1 ) * e[j] );
        od;
    od;
    return QM;
end;;

MakeGAlex3:=function(Q)
local B,LL,t,G,x,QQ,dG;
B:=TransposedMat(Q);
LL:=List(B,t->PermList(t));
G:=Group(LL);
x:=PermList(B[1]);
QQ:=MakeGAlex2(x,G);
dG:=DerivedSubgroup(G);
return [QQ,StructureDescription(Intersection(dG,Centralizer(G,x)))];
end;;

BraidFunction:=function(a,B,f)  #a = color list,B = braid, f = quandle, g = dual f
local n,i,v,k,S,N,j,g;
  g:=DualQuandle(f);
  n:=Maximum(List(B,AbsoluteValue))+1;
  N:=Length(B)+1;
  v:=NullMat(N,n);
   for i in [1..n] do 
     v[1][i] := a[i];
   od;
   for k in [2..N] do 
    i := B[k-1]; 
    for j in [1..n] do 
     if (j <> AbsoluteValue(i)) or (j <> AbsoluteValue(i)+1) then v[k][ j] := v[k-1][ j]; fi;
    od; 
    if i > 0 then 
       v[k][ i]:=  v[k-1][ i+1];
       v[k][ i+1] := f[v[k-1][ i]][ v[k-1][ i+1]];
    elif i < 0 then 
      i := -i; 
      v[k][ i]:= g[v[k-1][ i+1]][ v[k-1][ i]];
      v[k][ i+1] := v[k-1][ i];    
    fi; 
  od;  
 return v[N];
 end;;

Gindices2:=function(Q)
local T,W,i,j,g;
g:=Inn(Q);
T:=TransposedMat(Q);
W:=[];
for i in [1..Size(T)] do
 if T[1] = T[i] then Add(W,i); fi;
od;;
for i in [2..Size(W)] do
for j in [i+1..Size(W)] do
  if Set(Orbit(g,[1,W[i]],OnPairs)) = Set(Orbit(g,[1,W[j]],OnPairs)) then Remove(W,j); fi;
od;
od;
return W;
end;;

GindicesPerm:=function(Q)
 local  T, W, i,j,p,G,s,g; 
    g:=Inn(Q);  
    T := TransposedMat( Q );
    W := [  ];
    for i  in [ 1 .. Size( T ) ]  do
        if T[1] = T[i]  then
            Add( W, i );
        fi;
    od;
    p:=[];
    for i in [1..Size(W)] do
     s:=W[i];
     for j in [1..Size(W)] do
      if Set( Orbit( g, [ s,1 ], OnPairs ) ) = Set( Orbit( g, [ 1, W[j] ], OnPairs ) ) then p[i]:=Position(W,W[j]); fi;
     od;
   od;
   return [W,PermList(p)];      
end;;              

Gindices:=function ( Q )
    local  T, W, i;   
    T := TransposedMat( Q );
    W := [  ];
    for i  in [ 1 .. Size( T ) ]  do
        if T[1] = T[i]  then
            Add( W, i );
        fi;
    od;
    return W;
end;;

#the following works for the knot 12a_0690 = Knot[1491] only
NumCol:=function(n,f,g,z)
local N,y;
N:=0;
for y in [1..n] do
  if 
f[f[f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]][g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]]][g[g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]][f[f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]][g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]]]] = f[g[g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]][f[f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]][g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]]]][g[z][f[g[g[1][y]][f[y][g[1][y]]]][z]]]

and 

g[z][f[g[g[1][y]][f[y][g[1][y]]]][z]] = f[y][f[f[f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]][g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]]][g[g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]][f[f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]][g[f[g[g[1][y]][f[y][g[1][y]]]][z]][f[f[y][g[1][y]]][g[g[1][y]][f[y][g[1][y]]]]]]]]]
then 
   N:=N+1; 
fi;
od;
return N;
end;;

#the following works for the knot 2-bridge knot S(7,2) only
#the quandle has order n.
NColS72:=function(n,f,g)
local N,x,u0,u1,u2,u3,u4,u5,u6;
N:=0;
for x in [1..n] do
 u0:=1;
 u1:=x;
 u2:=g[u0][u1];
 u5:=g[u2][u0];
 u4:=f[u5][u1];
 u3:=f[u4][u0];
 u6:=g[u3][u1];
 if f[u1][u0] = u6 then N:=N+1; fi;  
od;
return N*n-n;
end;;



Psi:=function ( B, f )
    local  W, n, q, b, N, t, y, enum, i, yy, xx, w, z, flag, j;
    W := Gindices( f );
    n := Length( B );
    q := Length( f );
    b := Maximum( List( B, AbsoluteValue ) ) + 1;
    N := [  ];
    for t  in [ 1 .. Size( W ) ]  do
        N[t] := 0;
    od;
    enum := EnumeratorOfTuples( [ 1 .. q ], b - 1 );
    i := 0;
    while i < q ^ (b - 1)  do
        i := i + 1;
        yy := enum[i];
        xx := ShallowCopy( yy );
        w := Concatenation( [ 1 ], xx );
        z := BraidFunction( w, B, f );
        flag := true;
        for j  in [ 2 .. b ]  do
            if z[j] <> w[j]  then
                flag := false;
                break;
            fi;
        od;

        if flag = true  then
         N[Position( W, z[1] )] := N[Position( W, z[1] )] + 1;
        fi;
    od;
    return N;
end;;