##########################################################
#                                                        #
#  Updated on July 24, 2006:                             #
#                                                        #
#  The multi-component Gauss code to quandle             #
#  presentation converter had a bug. Rather than fixing  #
#  the quandle version, I've added a new program which   #
#  uses a faster counting algorithm which converts a     #
#  multicomponent Gauss code to a reduced biquandle      #
#  presentation and computes homomorphisms into the      # 
#  target quandle considered as a biquandle.             #
#                                                        #
#  To use this code you'll need both quandles-maple.txt  #
#  and biquandles-maple.txt saved in the same directory  #
#                                                        #
##########################################################

read("quandles-maple.txt");
read("biquandles-maple.txt");

gfix := proc (v)
#
# Fix Gauss code by doing type I moves
#
l:=coldim(matrix([v]))/2;
c:=true;
w:=[];
for i from 1 to 2*l do
 w:=[op(w),v[i]];
od;
while c do 
 c:=false;
 f:=array(1..l);   # position in v of undercrossing label of ith arc
 m:=array(1..2*l); # arc labels at crossings
 b:=array(1..l);   # overcrossing arc of ith arc
 j:=0;
 for i from 1 to 2*l do
  if Re(w[i])<0 then
   j:=j+1;
   f[j]:=i;
  fi;
   m[i]:=eval(j);
 od;
 for i from 1 to 2*l do
  if m[i]=0 then
   m[i]:=eval(l);
  fi;
 od;
 for i from 1 to l do
  j:=1;
  while w[j]+w[f[i]]<> 0 do
   j:=j+1;
  od;
  b[i]:=eval(m[j]);
 od;
 for i from 1 to l-1 do
  if not c and Im(w[f[i]])<>0  and Im(w[f[i+1]])=0 
   and b[m[f[i]]]=b[m[f[i+1]]] then
    z:=[];
    for k from 1 to f[i] do
     z:=[op(z),eval(w[k])];
    od;
    z:=[op(z),eval(l)+1];
    z:=[op(z),-1*eval(l)-1];
   if f[i]<>2*l then
    for k from f[i]+1 to 2*l do
     z:=[op(z),eval(w[k])];
    od;
   fi;
   c:=true;
   l:=l+1;
  fi;
 od;
 if not c and Im(w[f[l]])<>0 and Im(w[f[1]])=0 
  and b[m[f[l]]]=b[m[f[1]]] then
   z:=[];
   for k from 1 to f[l] do
    z:=[op(z),eval(w[k])];
   od;
   z:=[op(z),eval(l+1)];
   z:=[op(z),eval(-l-1)];
   if f[i]<>2*l then
    for k from f[i]+1 to 2*l do
     z:=[op(z),eval(w[k])];
    od;
   fi;
   c:=true;
   l:=l+1;
 fi;
 if c then
  w:=[];
  for i from 1 to 2*l do
   w:=[op(w),z[i]];
  od;
 fi;
od;
eval(w);
end;

gauss2presold := proc(v)
#
# converts Gauss code to upper quandle presentation matrix
#
n:=coldim(matrix([v]))/2;
f:=array(1..n);   # position in v of undercrossing label of ith arc
m:=array(1..2*n+1); # arc labels at crossings
b:=array(1..n);   # overcrossing arc of ith arc
j:=0;
for i from 1 to 2*n do
 if Re(v[i])<0 then
  j:=j+1;
  f[j]:=i;
 fi;
 if j>n then
  j:=1;
 fi;   
 m[i]:=eval(j);
od;
for i from 1 to 2*n do
 if m[i]= 0 then
  m[i]:=eval(j);
 fi;
od;
m[2*n+1]:=eval(m[1]);
for i from 1 to n do
 J:=1;
 while v[J]+v[f[i]]<> 0 do
  J:=J+1;
 od;
 b[i]:=eval(m[J]);
od;
temp:=[];
for i from 1 to n do
 temp:=[op(temp),0];
od;
T:=[];
for i from 1 to n do
 T:=[op(T),eval(temp)];
od;
M:=matrix(T);
for i from 2 to n do
 if Im(v[f[i]])<>0 then
  M[i,b[i]]:=i-1;
 else
  M[i-1,b[i]]:=i;
 fi;
od;
if Im(v[f[1]])<>0 then
 M[1,b[1]]:=n;
else 
 M[n,b[1]]:=1;
fi;
evalm(M);
end;


q2chom :=proc(v,T)
#
# compute 2-component quandle counting invariant for Gauss 
# code v and target quandle T
#
U:=homcount(gauss2pres(gfix(v)),T);
w:=[];;
for i from 1 to coldim(matrix([v])) do
 w:=[op(w),-1*eval(v[i])];
od;
x:=0;
while Re(w[1])<>-1 do
 x:=w[1];
 w:=subsop(1=NULL,w);
 w:=[op(w),eval(x)];
od;  
L:= homcount(gauss2pres(gfix(w)),T);
eval([U,L]);
end;

qdiff :=proc(v,T)
#
# compute quandle counting difference for Gauss 
# code v and target quandle T
#
U:=homcount(gauss2pres(gfix(v)),T);
w:=[];;
for i from 1 to coldim(matrix([v])) do
 w:=[op(w),-1*eval(v[i])];
od;
x:=0;
while Re(w[1])<>-1 do
 x:=w[1];
 w:=subsop(1=NULL,w);
 w:=[op(w),eval(x)];
od;  
L:= homcount(gauss2pres(gfix(w)),T);
eval(U-L);
end;

glist := proc(n)
#
# list single-component Gauss codes with n crossings
#
T:=[];
U:=[];
for i from 1 to n do
 T:=[op(T),[0,I]];
od;
S:=cartprod(T);
while not S[finished] do # S[nextvalue]()
 U:=[op(U),S[nextvalue]()];
od;
A:=permute(2*n);
B:=[];
for v in U do
 for i from 1 to (2*n)! do
  w:=A[i];
  t:=[];
  if w[1]=1 then
   for j from 1 to 2*n do
    if w[j] mod 2 =1 then 
     z:=(eval(w[j])+1)/2;
     t:=[op(t),z+eval(v[z])];
    elif w[j] mod 2 =0 then
     z:=eval(w[j])/2; 
     t:=[op(t),-1*z-eval(v[z])];
    fi;
   od;
   B:=[op(B),eval(t)];
  fi;
 od;
od;
eval(B);
end;

redI := proc(A)
#
# return the sublist of Gauss codes with no type-I reducibility
#
B:=[];
n:=coldim(matrix([A[1]]));
for a in A do
 c:=false;
 if a[1]+a[n]=0 then 
  c:=true;
 fi;
 for i from 1 to n-1 do
  if a[i]+a[i+1]=0 then
   c:=true; 
   break;
  fi;
 od;
 if not c then
  B:=[op(B),eval(a)];
 fi;
od;
eval(B);
end;

redII := proc(A)
#
# return the sublist of Gauss codes with no type-II reducibility
#
B:=[];
n:=coldim(matrix([A[1]]));
for v in A do
 c:=false;
 a:=[];
 for i from 1 to n do
  a:=[op(a),eval(v[i])];
 od;
 a:=[op(a),eval(v[1])];
 for i from 1 to n do
  if ((Re(a[i])>0 and Re(a[i+1])>0) or (Re(a[i])<0 and Re(a[i+1])<0)) 
   and (Im(a[i])<> Im(a[i+1])) then
   for j from 1 to n do
    if (a[j] = -1*eval(a[i]) and a[j+1] = -1*a[i+1]) or
       (a[j] = -1*eval(a[i+1]) and a[j+1] = -1*a[i]) then
     c:=true;
    fi;
   od;
  fi;
 od;
 if not c then
  a:=subsop(n+1=NULL,a);
  B:=[op(B),eval(a)];
 fi;
od;
eval(B);
end;

noevenint := proc(A)
#
# return the sublist of Gauss codes which are not evenly intersticed
#
B:=[];
n:=coldim(matrix([A[1]]));
for a in A do
 c:=true;
 for i from 1 to n-1 do 
  if c then
   j:=1;
   while a[i]+a[j]<>0 do
    j:=j+1;
   od;
   if (j-1) mod 2 =0 then
    c:=false;
   fi; 
  fi;
 od; 
 if not c then
  B:=[op(B),eval(a)];
 fi;
od;
eval(B);
end;

rglist:= proc(n)
#
# generates list of n-crossing Gauss codes with
# no obvious classicality or reducibility 
#
eval(noevenint(redI(redII(glist(n)))));
end;

qdcount:=proc(L,T)
#
# counts how many codes have non-classicality detected by QD_T
#
j:=0;
for i from 1 to coldim(matrix([L])) do
 if qdiff(L[i],T)<>0 then
  j:=j+1;
 fi;
od;
eval(j);
end;

gauss2list2:=proc(v)
#
# biquandle presentation list from Gauss code
#
v2:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  v2:=[op(v2),eval(v[i])];
 fi;
od;
b:=[]; # position in v2 of corresponding crossing label
for i from 1 to coldim(matrix([v2])) do
 j:=1;
 while v2[i]+v2[j]<>0 do
  j:=j+1;
 od;
 b:=[op(b),eval(j)];
od;
nx:=[];
lt:=1;
c:=0;
for i from 1 to coldim(matrix([v]))-1 do
 if v[i+1]=0 then
  nx:=[op(nx),eval(lt)];
  lt:=eval(i+1-c);
  c:=c+1;
 elif v[i]<>0 then
  nx:=[op(nx),eval(i+1-c)];
 fi;
od;
out:=[];
for i from 1 to coldim([v2]) do
 if Re(v2[i]) >0 and Im(v2[i])=0 then
  out:=[op(out),[M[4][eval(x[i]),eval(x[b[i]])], eval(x[eval(nx[i])])]];
 elif Re(v2[i]) >0 and Im(v2[i])<>0 then
  out:=[op(out),[M[3][eval(x[i]),eval(x[b[i]])], eval(x[eval(nx[i])])]];
 elif Re(v2[i]) <0 and Im(v2[i])=0 then
  out:=[op(out),[M[2][eval(x[i]),eval(x[b[i]])], eval(x[eval(nx[i])])]];
 elif Re(v2[i]) <0 and Im(v2[i])<>0 then
  out:=[op(out),[M[1][eval(x[i]),eval(x[b[i]])], eval(x[eval(nx[i])])]];
 fi;
od;
eval(out);
end;

reducepreslist:=proc(A)
#
# reduce biquandle presentation list
#
L:=[];
for i from 1 to coldim(matrix([A])) do
 L:=[op(L),eval(A[i])];
od;
c:=true;
while c do
 c:=false;
 for i from 1 to coldim(matrix([L])) do
  w:=eval(L[1]);
  L:=subsop(1=NULL, L);
  if subs(eval(w[2])=z,eval(w[1]))=eval(w[1]) then
   c:=true;
   L:=subs(eval(w[2])=eval(w[1]),L);
  else L:=[op(L),eval(w)];
  fi;
 od;
od;
eval(L);
end;


bhomlist2:=proc(v,A)
#
# List homomorphisms from the biquandle of the Gauss code v
# into the target biquandle A by reducing the presentation
# then reconstruct full coloring
#
out:=[];
L1:=gauss2list2(v);
L2:=reducepreslist(L1);
L:=eval(L2);
for i from 1 to coldim(matrix([L])) do
 L:=subs(eval(L[i][2])=f[i],L);
od;
F:=listmaps(coldim(matrix([L])),coldim(A[1]));
for i from 1 to coldim(A[1])^coldim(matrix([L])) do
 w:=eval(F[i]);
 c:=true;
 for j from 1 to coldim(matrix([L])) do
  X:=L[j][1];
  Y:=L[j][2];
  if eval(subs(f=w,M=A,eval(X)))<>eval(subs(f=w,M=A,eval(Y))) then
   c:=false;
  fi;
 od;
 if c then
  out:=[op(out),eval(w)];
 fi;
od;
out2:=[];
blk:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  blk:=[op(blk),0];
 fi;
od;
B:=gauss2biq(v);
for C in out do
 t:=eval(blk);
 for i from 1 to coldim(matrix([L2])) do
  for j from 1 to coldim(matrix([blk])) do
   if L2[i][2]=x[j] then
    t[j]:=eval(C[i]);
   fi;
  od;
 od;
 out2:=[op(out2),eval(bhomcomplete(B,A,t))];
od;
eval(out2);
end;

bhomcomplete:=proc(B,A,t)
#
# t:B-> A
#
v:=eval(t);
for i from 1 to coldim(B[1]) do
 if v[i]=0 then
  j:=1;
  w:= eval(v);
  w[i]:=j;
  while type(bhomfill(B,A,w),boolean) do
   j:=j+1;
   w:= eval(v);
   w[i]:=j;
  od;
  v:=eval(w);
 fi;
od;
eval(v);
end;

gauss2biq:=proc(v)
#
# convert gauss code to biquandle presentation matrix
#
w:=[];
n:=coldim(matrix([v]));
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  w:=[op(w),0];
 fi;
od;
W:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  W:=[op(W),eval(w)];
 fi;
od;
v2:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  v2:=[op(v2),eval(v[i])];
 fi;
od;
A:=[];
for i from 1 to 4 do
 A:=[op(A),evalm(W)];
od;
b:=[]; # position in v2 of corresponding crossing label
for i from 1 to coldim(matrix([v2])) do
 j:=1;
 while v2[i]+v2[j]<>0 do
  j:=j+1;
 od;
 b:=[op(b),eval(j)];
od;
nx:=[];
lt:=1;
c:=0;
for i from 1 to coldim(matrix([v]))-1 do
 if v[i+1]=0 then
  nx:=[op(nx),eval(lt)];
  lt:=eval(i+1-c);
  c:=c+1;
 elif v[i]<>0 then
  nx:=[op(nx),eval(i+1-c)];
 fi;
od;
for i from 1 to coldim(matrix([v2])) do
 if Re(v2[i]) >0 and Im(v2[i])=0 then
  A[4][i,b[i]]:=eval(nx[i]);
 elif Re(v2[i]) >0 and Im(v2[i])<>0 then
  A[3][i,b[i]]:=eval(nx[i]);
 elif Re(v2[i]) <0 and Im(v2[i])=0 then
  A[2][i,b[i]]:=eval(nx[i]);
 elif Re(v2[i]) <0 and Im(v2[i])<>0 then
  A[1][i,b[i]]:=eval(nx[i]);
 fi;
od;
eval(A);
end;


qhomcount:=proc(v,M)
#
# count homomorphisms from knot quandle of Gauss code v 
# to finite quandle M by converting to biquandles
#
t:=[];
z:=[];
T:=[];
Z:=[];
for i from 1 to coldim(M) do
 t:=[op(t),eval(i)];
 z:=[op(z),0];
od;
for i from 1 to coldim(M) do
 T:=[op(T),eval(t)];
 Z:=[op(Z),eval(z)];
od;
T2:=transpose(matrix(T));
for i from 1 to coldim(M) do
 for j from 1 to coldim(M) do
  for k from 1 to coldim(M) do
   if M[i,j]=k then
    Z[k,j]:=i;
   fi;
  od;
 od;
od;
out:=coldim(matrix([bhomlist2(v,[evalm(Z),evalm(M),evalm(T2),evalm(T2)])]));
eval(out);
end;