###############################################
#
# This maple code uses biquandles-maple.txt
#
# Permission to redistribute and reimplement
# granted. Bugfixes will be made as necessary.
#
###############################################

read("biquandles-maple.txt");

pairop:=proc(x,y,A)
#
# (x,y) |--> alpha x + (x M y^t) y
# where A = [alpha, M]
#
M:=evalm(matrix([x])&*A[2]&*transpose(matrix([y])));
evalm(A[1]*x+M[1,1]*y);
end;

pairopn:=proc(x,y,A,n)
#
# (x,y) |--> alpha x + (x M y^t) y  mod n
# where A = [alpha, M]
#
M:=evalm(matrix([x])&*A[2]&*transpose(matrix([y])));
X:=A[1]*x+M[1,1]*y;
out:=[];
for i from 1 to coldim(matrix([X])) do
 out:=[op(out),X[i] mod n];
od;
eval(out);
end;


bv:=proc(j,n)
#
# jth stadard basis vector in znm
#
out:=[];
for i from 1 to j-1 do
 out:=[op(out),0];
od;
out:=[op(out),1];
for i from j+1 to n do
 out:=[op(out),0];
od;
eval(out);
end;

znm:=proc(n,m)
#
# Create list of (Zn)^m elements
#
t1:=[];
for i from 0 to n-1 do
 t1:=[op(t1),eval(i)];
od;
t2:=[];
for i from 1 to m do
 t2:=[op(t2),eval(t1)];
od;
V:=[];
T:=cartprod(t2): while not T[finished] do 
 J:=T[nextvalue]();
 V:=[op(V),J]; 
od;
eval(V);
end;


bfind:=proc(A,n)
#
# find all bilinear biquandle matrices with first matrix
# determines by the A=[alpha,M]
#
m:=coldim(A[2]);
L:=znm(n,m);
br:=[];
for i from 1 to n^m do
 br:=[op(br),0];
od;
bm:=[];
for i from 1 to n^m do
 bm:=[op(bm),br];
od;
M:=matrix(bm);
P:=[evalm(M),evalm(M),evalm(M),evalm(M)];
for i from 1 to n^m do
 for j from 1 to n^m do
  en:=pairopn(L[i],L[j],A,n);
  for k from 1 to n^m do
   if en=L[k] then
    P[1][i,j]:=k;
   fi;
  od;
 od;
od;
eval(biqlist(P));
end;

bdecode:=proc(M,n,m)
#
# recover pair operator from biquandle matrix assuming
# znm order of module elements
#
br:=[];
for i from 1 to m do 
 br:=[op(br),0];
od;
bm:=[];
for i from 1 to m do
 bm:=[op(bm),br];
od;
X:=matrix(bm);
L:=znm(n,m);
a:=L[M[bvfind(1,m,L),bvfind(2,m,L)]][1];
for i from 1 to m do
 for j from 1 to m do
  if i<> j then
   X[i,j]:=L[M[bvfind(i,m,L),bvfind(j,m,L)]][j];
  else
   X[i,j]:=L[M[bvfind(i,m,L),bvfind(j,m,L)]][j]-a mod n;
  fi;
 od;
od;
eval([a,evalm(X)]);
end;

bvfind:=proc(i,n,L)
#
# find which row/column correspond to ith basis vector
#
v:=bv(i,n);
out:=0;
for j from 1 to coldim(matrix([L])) do
 if v=L[j] then
  out:=j;
 fi;
od;
eval(out);
end;

bilbiq:=proc(n,A,B,C,DD)
#
# make operation matrix for (Z_n)^m with biquandle pair operations
#
m:=coldim(A[2]);
temp:=[];
temp2:=[];
for i from 1 to n^m do
 temp:=[op(temp),0];
od;
for i from 1 to n^m do
 temp2:=[op(temp2),eval(temp)];
od;
M:=matrix(temp2);
out:=[evalm(M),evalm(M),evalm(M),evalm(M)];
V:=znm(n,m);
for i from 1 to n^m do
 for j from 1 to n^m do
  x1:=pairopn(V[i],V[j],A,n);
  x2:=pairopn(V[i],V[j],B,n);
  x3:=pairopn(V[i],V[j],C,n);
  x4:=pairopn(V[i],V[j],DD,n);
  for k from 1 to n^m do
   if x1=V[k] then
    out[1][i,j]:=eval(k);
   fi;
   if x2=V[k] then
    out[2][i,j]:=eval(k);
   fi;
   if x3=V[k] then
    out[3][i,j]:=eval(k);
   fi;
   if x4=V[k] then
    out[4][i,j]:=eval(k);
   fi;
  od;
 od;
od;
eval(out);
end;


subbiq:=proc(L1,B)
#
# sub-biquandle of B spanned by L1
#
L:=eval(L1);
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;

vssubbiq:=proc(L1,B,n)
#
# submodule of B spanned by L1
#
L2:=subbiq(L1,B);
L:=znm(n,log[n](coldim(B[1])));
S:={};
for x in L2 do 
 S:=S union {eval(L[x])};
od;
c:=true;
while c do
 c:=false;
 for x in S do
  for y in S do
   w:=x+y mod n;
   if not w in S then
    c:=true;
    S:=S union {eval(w)};
   fi;
  od;
  for i from 0 to n-1 do
   w:=i*x mod n;
   if not w in S then
    c:=true;
    S:=S union {eval(w)};
   fi;
  od;
 od;
od;
eval(S);
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;


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;

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;


bbinv:=proc(g,B,n)
#
# bilinear biquandle invariant
#
H:=bhomlist2(g,B);
out:=0;
for f in H do
 S:=convert(f,set);
 sb:=coldim(matrix([convert(subbiq(S,B),list)]));
 sv:=coldim(matrix([convert(vssubbiq(S,B,n),list)]));
 out:=out+q^sb*z^sv;
od;
eval(out);
end;

bilfind2:=proc(n)
#
# bilinear biquandle structures on (Z_n)^2
#
InvEl:={};
for i from 1 to n-1 do
 for j from 1 to i do
  if i*j mod n =1 then
    InvEl:=InvEl union {i};
  fi;
 od;
od;
out:=[];
for alpha in InvEl do
 for beta in InvEl do
  ai:=1/alpha mod n;
  bi:=1/beta mod n;
  alist:=[];
  for f from 0 to n-1 do
   if alpha*(1-beta^2)*f mod n = 0 and
      beta*(1-beta^2)*f mod n = 0 then 
    alist:=[op(alist),eval(f)];
   fi;
  od;
  omega:=ai^2-beta*ai-(ai*bi)^2 mod n;
  for f1 in alist do
   for f2 in alist do
    M:=matrix([[bi-alpha mod n,f1],[f2,bi-alpha mod n]]);
    M2:=matrix([[omega*(bi-alpha) mod n,omega*f1 mod n],
                [omega*f2 mod n,omega*(bi-alpha) mod n]]);
    Z:=matrix([[0,0],[0,0]]);
    if biqtest2(bilbiq(n,[ai,M2],[alpha,M],[bi,Z],[beta,Z])) then 
     out:=[op(out),[alpha,evalm(M),beta]];
    fi;
   od;
  od;
 od;
od;
eval(out);
end;

bilfind3:=proc(n)
#
# bilinear biquandle structures on (Z_n)^3
#
InvEl:={};
for i from 1 to n-1 do
 for j from 1 to i do
  if i*j mod n =1 then
    InvEl:=InvEl union {i};
  fi;
 od;
od;
out:=[];
for alpha in InvEl do
 for beta in InvEl do
  ai:=1/alpha mod n;
  bi:=1/beta mod n;
  alist:=[];
  for f from 0 to n-1 do
   if alpha*(1-beta^2)*f mod n = 0 and
      beta*(1-beta^2)*f mod n = 0 then 
    alist:=[op(alist),eval(f)];
   fi;
  od;
  omega:=ai^2-beta*ai-(ai*bi)^2 mod n;
  for f1 in alist do
   for f2 in alist do
    for f3 in alist do
     for f4 in alist do
      for f5 in alist do
       for f6 in alist do
        M:=matrix([[bi-alpha mod n,f1,f2],
                   [f3,bi-alpha mod n,f4],
                   [f5,f6,bi-alpha mod n]]);
        M2:=matrix([[omega*(bi-alpha) mod n,omega*f1 mod n,omega*f2 mod n],
                    [omega*f3 mod n,omega*(bi-alpha) mod n,omega*f4 mod n],
                    [omega*f5 mod n,omega*f6 mod n,omega*(bi-alpha) mod n]]);
        Z:=matrix([[0,0,0],[0,0,0],[0,0,0]]);
        if biqtest2(bilbiq(n,[ai,M2],[alpha,M],[bi,Z],[beta,Z])) then 
             out:=[op(out),[alpha,evalm(M),beta]];
        fi;
       od; 
      od;
     od;
    od;
   od;
  od;
 od;
od;
eval(out);
end;

bilbiq2:=proc(alpha,beta,M,n)
#
# Biquandle operation matrix from minimal information
#
ai:=1/alpha mod n;
bi:=1/beta mod n;
omega:=ai^2-beta*ai-(ai*bi)^2 mod n;
m:=[];
z:=[];
for i from 1 to coldim(M) do
 r:=[];
 zr:=[];
 for j from 1 to coldim(M) do
  r:=[op(r),eval(omega*M[i,j] mod n)];
  zr:=[op(zr),0];
 od;
 m:=[op(m),eval(r)];
 z:=[op(z),eval(zr)];
od;
M2:=matrix(m);
Z:=matrix(z);
bilbiq(n,[ai,M2],[alpha,M],[bi,Z],[beta,Z]);
end;


bisolist2 := proc(A,B)
#
# lists isomorphisms from A to B
#
n:=coldim(A[1]);
m:=coldim(B[1]);
out:=[];
st:=[];
for i from 1 to n do
 st:=[op(st),0];
od;
w:=[eval(st)];
c:=true;
while c do
 c:=false;
 for j from 1 to coldim(matrix([w])) do
  v:=eval(w[1]);
  w:=subsop(1=NULL,w);
  l:=0;
  for i from 1 to n do
   if v[i]=0 then
    l:=i;
    break;
   fi;
  od;
  x:=eval(v);
  if l=0 then
   out:=[op(out),eval(x)];
  else
   for i from 1 to m do
    x[l]:=i;
    z:=bisofill(A,B,x);
    if type(z,list) then
     w:=[op(w),eval(z)];
     c:=true;
    fi;
   od;
  fi;
 od;
od;
eval(out);
end;


bisofill:=proc(A,B,v)
#
# v:A -> B
#
c:=true;
ret:=true;
w:=eval(v);
while c do
 c:=false;
 for i from 1 to eval(coldim(A[1])) do
  for j from 1 to eval(coldim(A[1])) do
   if w[i]<>0 and w[j]<>0 and A[1][i,j]<>0 and
      w[A[1][i,j]]<> B[1][w[i],w[j]] then
    if w[A[1][i,j]]=0 then
     w[A[1][i,j]]:=eval(B[1][w[i],w[j]]);
     c:=true;
    else
     ret:=false;
    fi;
   fi;
   if w[i]<>0 and w[j]<>0 and A[2][i,j]<>0 and
      w[A[2][i,j]]<> B[2][w[i],w[j]] then
    if w[A[2][i,j]]=0 then
     w[A[2][i,j]]:=eval(B[2][w[i],w[j]]);
     c:=true;
    else
     ret:=false;
    fi;
   fi;
   if w[i]<>0 and w[j]<>0 and A[3][i,j]<>0 and
      w[A[3][i,j]]<> B[3][w[i],w[j]] then
    if w[A[3][i,j]]=0 then
     w[A[3][i,j]]:=eval(B[3][w[i],w[j]]);
     c:=true;
    else
     ret:=false;
    fi;
   fi;
   if w[i]<>0 and w[j]<>0 and A[4][i,j]<>0 and
      w[A[4][i,j]]<> B[4][w[i],w[j]] then
    if w[A[4][i,j]]=0 then
     w[A[4][i,j]]:=eval(B[4][w[i],w[j]]);
     c:=true;
    else
     ret:=false;
    fi;
   fi;
  od;  
 od;
 if not permtest(w) then
  c:=false;
  ret:=false;
 fi;
od;
if ret then
 eval(w);
else
 eval(ret);
fi;
end;


bisotest2 := proc(A,B)
#
# tests whether A is isomorphic to B
#
n:=coldim(A[1]);
m:=coldim(B[1]);
out:=[];
iso:=false;
st:=[];
for i from 1 to n do
 st:=[op(st),0];
od;
w:=[eval(st)];
c:=true;
while c do
 c:=false;
 for j from 1 to coldim(matrix([w])) do
  v:=eval(w[1]);
  w:=subsop(1=NULL,w);
  l:=0;
  for i from 1 to n do
   if v[i]=0 then
    l:=i;
    break;
   fi;
  od;
  x:=eval(v);
  if l=0 then
   out:=[op(out),eval(x)];
  else
   for i from 1 to m do
    x[l]:=i;
    z:=bisofill(A,B,x);
    if type(z,list) then
     w:=[op(w),eval(z)];
     c:=true;
    fi;
   od;
  fi;
 od;
 if coldim(matrix([out])) <>0 then
  c:=false;
  iso:=true;
 fi;
od;
eval(iso);
end;

bilfind4:=proc(n)
#
# bilinear biquandle structures on (Z_n)^4
#
InvEl:={};
for i from 1 to n-1 do
 for j from 1 to i do
  if i*j mod n =1 then
    InvEl:=InvEl union {i};
  fi;
 od;
od;
out:=[];
for alpha in InvEl do
 for beta in InvEl do
  ai:=1/alpha mod n;
  bi:=1/beta mod n;
  alist:=[];
  for f from 0 to n-1 do
   if alpha*(1-beta^2)*f mod n = 0 and
      beta*(1-beta^2)*f mod n = 0 then 
    alist:=[op(alist),eval(f)];
   fi;
  od;
  omega:=ai^2-beta*ai-(ai*bi)^2 mod n;
  for f1 in alist do
   for f2 in alist do
    for f3 in alist do
     for f4 in alist do
      for f5 in alist do
       for f6 in alist do
        for f7 in alist do
         for f8 in alist do
          for f9 in alist do
           for f10 in alist do
            for f11 in alist do
             for f12 in alist do
        M:=matrix([[bi-alpha mod n,f1,f2,f3],
                   [f4,bi-alpha mod n,f5,f6],
                   [f7,f8,bi-alpha mod n,f9],
                   [f10,f11,f12,bi-alpha mod n]]);
       M2:=matrix([[omega*(bi-alpha) mod n,omega*f1 mod n,omega*f2 mod n,omega*f3 mod n],
                   [omega*f4 mod n,omega*(bi-alpha) mod n,omega*f5 mod n,omega*f6 mod n],
                   [omega*f7 mod n,omega*f8 mod n,omega*(bi-alpha) mod n,omega*f9 mod n],
                   [omega*f10 mod n,omega*f11 mod n,omega*f12 mod n,omega*(bi-alpha) mod n]]);
        Z:=matrix([[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]);
        if biqtest2(bilbiq(n,[ai,M2],[alpha,M],[bi,Z],[beta,Z])) then 
             out:=[op(out),[alpha,evalm(M),beta]];
        fi;
             od;
            od;
           od;
          od;
         od;
        od;
       od; 
      od;
     od;
    od;
   od;
  od;
 od;
od;
eval(out);
end;