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

qpnm:=proc(M,n,m)
#
# (m,n) quandle polynomial
#
f:=0;
for x from 1 to coldim(M) do
 cnx:=0;
 rmx:=0;
 for y from 1 to coldim(M) do
  if qit(M,x,y,n)=x then
   cnx:=cnx+1;
  fi;
 od;  
 for y from 1 to coldim(M) do
  if qit(M,y,x,m)=y then
   rmx:=rmx+1;
  fi;
 od;
 f:=f+s^cnx*t^rmx;
od;
eval(f);
end;

qit:=proc(M,x,y,n)
#
# iterate quandle operation
#
Md:=qdual2(M);
z:=x;
if n>0 then
 for i from 1 to n do
  z:=M[z,y];
 od;
 elif n<0 then
  for i from 1 to -1*n do
  z:=Md[z,y];
 od;
fi;
eval(z);
end;

permorder:=proc(x)
#
# find order of permuatation
#
count:=1;
y:=eval(x);
c:=true;
while c do
 c:=false;
 z:=eval(y);
 for i from 1 to nops(x) do
  y[i]:=x[z[i]];
 od;
 if eval(y)<>eval(x) then 
  c:=true;
  count:=count+1;
 fi;
od;
eval(count);
end;

qdual2 := 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);
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);
end;

qpmatrix:=proc(M)
#
# compute quandle polynomial matrix
#
l:=[];
for i from 1 to coldim(M) do
 x:=[];
 for j from 1 to coldim(M) do
  x:=[op(x),M[j,i]];
 od;
 l:=[op(l),permorder(x)];
od;
N:=1;
for i from 1 to nops(l) do
 N:=lcm(N,l[i]);
od;
out:=[];
for i from 1 to N do
 temp:=[];
 for j from 1 to N do
  temp:=[op(temp),qpnm(M,j-1,i-1)];
 od;
 out:=[op(out),temp];
od;
eval(matrix(out));
end;


bqpnm:=proc(B,n,m)
#
# biquandle polynomial
#
f:=0;
for x from 1 to coldim(B[1]) do
 cnx1:=0;
 rmx1:=0;
 cnx2:=0;
 rmx2:=0;
 cnx3:=0;
 rmx3:=0;
 cnx4:=0;
 rmx4:=0;
 for y from 1 to coldim(B[1]) do
  if qit(B[1],x,y,n)=x then
   cnx1:=cnx1+1;
  fi;
  if qit(B[2],x,y,n)=x then
   cnx2:=cnx2+1;
  fi;
  if qit(B[3],x,y,n)=x then
   cnx3:=cnx3+1;
  fi;
  if qit(B[4],x,y,n)=x then
   cnx4:=cnx4+1;
  fi;
 od;  
 for y from 1 to coldim(B[1]) do
  if qit(B[1],y,x,m)=y then
   rmx1:=rmx1+1;
  fi;
  if qit(B[2],y,x,m)=y then
   rmx2:=rmx2+1;
  fi;
  if qit(B[3],y,x,m)=y then
   rmx3:=rmx3+1;
  fi;
  if qit(B[4],y,x,m)=y then
   rmx4:=rmx4+1;
  fi;
 od;
 f:=f+s1^cnx1*t1^rmx1*s2^cnx2*t2^rmx2*s3^cnx3*t3^rmx3*s4^cnx4*t4^rmx4;
od;
eval(f);
end;

bpmatrix:=proc(B)
#
# compute biquandle polynomial matrix
#
l:=[];
for z from 1 to 4 do
 for i from 1 to coldim(B[1]) do
  x:=[];
  for j from 1 to coldim(B[1]) do
   x:=[op(x),B[z][j,i]];
  od;
  l:=[op(l),permorder(x)];
 od;
od;
N:=1;
for i from 1 to nops(l) do
 N:=lcm(N,l[i]);
od;
out:=[];
for i from 1 to N do
 temp:=[];
 for j from 1 to N do
  temp:=[op(temp),bqpnm(B,j-1,i-1)];
 od;
 out:=[op(out),temp];
od;
eval(matrix(out));
end;

sbqpnm:=proc(S,B,n,m)
#
# subbiquandle polynomial
#
f:=0;
for x in S do
 cnx1:=0;
 rmx1:=0;
 cnx2:=0;
 rmx2:=0;
 cnx3:=0;
 rmx3:=0;
 cnx4:=0;
 rmx4:=0;
 for y from 1 to coldim(B[1]) do
  if qit(B[1],x,y,n)=x then
   cnx1:=cnx1+1;
  fi;
  if qit(B[2],x,y,n)=x then
   cnx2:=cnx2+1;
  fi;
  if qit(B[3],x,y,n)=x then
   cnx3:=cnx3+1;
  fi;
  if qit(B[4],x,y,n)=x then
   cnx4:=cnx4+1;
  fi;
 od;  
 for y from 1 to coldim(B[1]) do
  if qit(B[1],y,x,m)=y then
   rmx1:=rmx1+1;
  fi;
  if qit(B[2],y,x,m)=y then
   rmx2:=rmx2+1;
  fi;
  if qit(B[3],y,x,m)=y then
   rmx3:=rmx3+1;
  fi;
  if qit(B[4],y,x,m)=y then
   rmx4:=rmx4+1;
  fi;
 od;
 f:=f+s1^cnx1*t1^rmx1*s2^cnx2*t2^rmx2*s3^cnx3*t3^rmx3*s4^cnx4*t4^rmx4;
od;
eval(f);
end;


subbiq:=proc(L1,B)
#
# sub-biquandle of B spanned by L1
#
L:=convert(L1,set);
c:=true;
while c do
 c:=false;
 for x in L do
  for y in L do
   if not B[1][x,y] in L then
    c:=true;
    L:= {eval(B[1][x,y])} union L;
   fi;
   if not B[2][x,y] in L then
    c:=true;
    L:= {eval(B[2][x,y])} union L;
   fi;
   if not B[3][x,y] in L then
    c:=true;
    L:= {eval(B[3][x,y])} union L;
   fi;
   if not B[4][x,y] in L then
    c:=true;
    L:= {eval(B[4][x,y])} union L;
   fi;
  od;
 od;
od;
eval(L);
end;

subq:=proc(L1,B)
#
# subquandle of B spanned by L1
#
L:=convert(L1,set);
c:=true;
while c do
 c:=false;
 for x in L do
  for y in L do
   if not B[x,y] in L then
    c:=true;
    L:= {eval(B[x,y])} union L;
   fi;
  od;
 od;
od;
eval(L);
end;

subbpinv:=proc(g,B,n,m)
#
# multiset version
#
out:=[];
H:=bhomlist2(g,B);
for f in H do 
 S:=[];
 for i from 1 to nops(f) do
  if f[i]<>0 then 
   S:=[op(S),f[i]];
  fi;
 od;
 IM:=subbiq2(S,B);
 out:=[op(out),sbqpnm(IM,B,n,m)];
od;
eval(out);
end;

specsubbpinv:=proc(g,B,n,m,x1,x2,x3,x4,y1,y2,y3,y4)
#
# specialized polynomial version
#
out:=0;
H:=bhomlist2(g,B);
for f in H do 
 S:=[];
 for i from 1 to nops(f) do
  if f[i]<>0 then 
   S:=[op(S),f[i]];
  fi;
 od;
 IM:=subbiq2(S,B);
 out:=out+
  q^(subs(s1=x1,s2=x2,s3=x3,s4=x4,t1=y1,t2=y2,t3=y3,t4=y4,sbqpnm(IM,B,n,m)));
od;
eval(out);
end;

mssubbpinv:=proc(g,B,x1,x2,x3,x4,y1,y2,y3,y4)
#
# matrix of specialized subbiqinvs
#
l:=[];
for z from 1 to 4 do
 for i from 1 to coldim(B[1]) do
  x:=[];
  for j from 1 to coldim(B[1]) do
   x:=[op(x),B[z][j,i]];
  od;
  l:=[op(l),permorder(x)];
 od;
od;
N:=1;
for i from 1 to nops(l) do
 N:=lcm(N,l[i]);
od;
out:=[];
for i from 1 to N do
 temp:=[];
 for j from 1 to N do
  temp:=[op(temp),specsubbpinv(g,B,j-1,i-1,x1,x2,x3,x4,y1,y2,y3,y4)];
 od;
 out:=[op(out),temp];
od;
eval(matrix(out));
end;

q2biq:=proc(M)
#
# convert quandle to biquandle
#
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:=[evalm(Z),evalm(M),evalm(T2),evalm(T2)];
end;