r := function(biquandle, x, y)
  return [biquandle.lperms[x][y], biquandle.rperms[y][x]];
end;



paste:=function(clases)
local i,j,c,cl,n;

cl:=clases;

for i in [1..Size(cl)] do
  cl[i]:=UnionSet(cl[i],[]);
 od; 

cl:=UnionSet(cl,[]);
n:=Size(cl);

for i in [1..n] do
  for j in [i+1..n] do
   if j>Size(cl) then continue;
   else
    c:=IntersectionSet(cl[i],cl[j]);
    if Size(c)>0 then
     cl[i]:=UnionSet(cl[i],cl[j]);
     if cl[i]=cl[j] then continue; 
     else cl:=Difference(cl, [cl[j]] );
     fi;
    fi;
   fi;
  od;
 od;

 return cl;
 end;

relationgenerated:=function(clase)
local i, cl;
 cl:=clase;
 for i in [1..Size(clase)] do
  cl:=paste(cl);
 od;
 return cl;
 end;





#A biquandle is given by rperms, lperms, but also contains the information
# about size, labels (just in case), and function s. 


# In case there is given only rperms and lperms, this functions
# add the information about size, s, and label.
oficialbiquandle:=function(bi)
local x,y,ss,N, labels;
ss:=[];
N:=Size(bi.lperms);
labels:=UnionSet(bi.lperms[1],[]);

for x in [1..N] do
 for y in [1..N] do
  if r(bi,x,y)=[x,y] then
   Add(ss,y);
  fi;
 od;
od;
return rec( lperms := bi.lperms, rperms := bi.rperms, size := N, labels := labels, s :=  ss);
end;



#Bialexander biquandle
bialexander := function(n, s, t)
  local e, lperms, rperms, x, y, ss;
  
  s := s*One(ZmodnZ(n));
  t := t*One(ZmodnZ(n));
  e := Enumerator(ZmodnZ(n));

  lperms := NullMat(Size(e), Size(e));
  rperms := NullMat(Size(e), Size(e));

  for x in e do
    for y in e do
      lperms[Position(e, x)][Position(e, y)] := Position(e, s*y);
      rperms[Position(e, y)][Position(e, x)] := Position(e, t*x+(1-s*t)*y);
    od;
  od;

  ss := [];
  for x in e do
    Add(ss, Position(e, Inverse(s)*x));
  od;

  return rec( lperms := lperms, rperms := rperms, size := n, labels := e, s :=  ss);
end;


#Bialexander in a field
bialexanderField := function(q, s,t) #q=p^n
  local e, lperms, rperms, x, y, ss;
  

  e := Enumerator(Field(Z(q)));

  lperms := NullMat(Size(e), Size(e));
  rperms := NullMat(Size(e), Size(e));

  for x in e do
    for y in e do
      lperms[Position(e, x)][Position(e, y)] := Position(e, s*y);
      rperms[Position(e, y)][Position(e, x)] := Position(e, t*x+(1-s*t)*y);
    od;
  od;

  ss := [];
  for x in e do
    Add(ss, Position(e, Inverse(s)*x));
  od;

  return rec( lperms := lperms, rperms := rperms, size := Size(Field(Z(q))), labels := e, s :=  ss);
end;

#flip coproduct
flipcop:=function(p,q)
local n,m,lperms,rperms,i,j,biq;

n:=Size(p.rperms);
m:=Size(q.rperms);


lperms := NullMat(n+m,n+m);
rperms := NullMat(n+m,n+m);

for i in [1..n] do
for j in [1..n] do
lperms[i][j]:= p.lperms[i][j];
rperms[i][j]:= p.rperms[i][j];
od;od;

for i in [1..m] do
for j in [1..m] do
lperms[i+n][j+n]:= n+q.lperms[i][j];
rperms[i+n][j+n]:= n+q.rperms[i][j];
od;od;

for i in [1..n] do
for j in [1..m] do
lperms[i][j+n]:= j+n;
lperms[j+n][i]:= i;
rperms[i][j+n]:= j+n;
rperms[j+n][i]:= i;
od;od;

biq:=rec(rperms:=rperms,lperms:=lperms);
return oficialbiquandle(biq);
end;


#One point quandle extension
1_extension:=function(q,f) #q es un quandle y f es un automorfismo de quandles tal que f(x*y)=f(x)*y
local n,m,lperms,rperms,i,j,biq;

n:=Size(q.rperms);


lperms := NullMat(n+1,n+1);
rperms := NullMat(n+1,n+1);

for i in [1..n] do
for j in [1..n] do
lperms[i][j]:= q.lperms[i][j];
rperms[i][j]:= q.rperms[i][j];
od;od;

lperms[n+1][n+1]:= n+1;
rperms[n+1][n+1]:= n+1;

for i in [1..n] do
lperms[i][n+1]:= n+1;
lperms[n+1][i]:= i;
rperms[i][n+1]:= n+1;
rperms[n+1][i]:= f[i];
od;

biq:=rec(rperms:=rperms,lperms:=lperms);
return biq;
end;



#Wadda's solution
wada := function(group)
  local x, y, e, s, lperms, rperms;

  e := Elements(group);
  s := [];
  lperms := NullMat(Size(group), Size(group));
  rperms := NullMat(Size(group), Size(group));
  for x in group do
    for y in group do
      lperms[Position(e, x)][Position(e, y)] := Position(e, x*Inverse(y)*Inverse(x));
      rperms[Position(e, y)][Position(e, x)] := Position(e, x*y^2);
    od;
  od;

  for x in e do
    Add(s, Position(e, Inverse(x)));
  od;

  return rec( lperms := lperms, rperms := rperms, size := Order(group), labels := e, s :=  s);
end;

#This biquandle is the quandle solution associated to a subset of a group
# that is stable under conjugation
conj := function(conjclass)
  local x, y, e, s, lperms, rperms;

  e := Elements(conjclass);
  s := [];
  lperms := NullMat(Size(conjclass), Size(conjclass));
  rperms := NullMat(Size(conjclass), Size(conjclass));
  for x in conjclass do
    for y in conjclass do
      lperms[Position(e, x)][Position(e, y)] := Position(e, y);
      rperms[Position(e, y)][Position(e, x)] := Position(e, Inverse(y)*x*y);
    od;
  od;

  for x in e do
    Add(s, Position(e, x));
  od;

  return rec( lperms := lperms, rperms := rperms, size := Size(conjclass), labels := e, s :=  s);
end;

##involutive Zn-biquandle
invo := function(n)
  local x, y, e, s, lperms, rperms;
  e := Elements(ZmodnZ(n));
  s := [];
  lperms := NullMat(n,n);
  rperms := NullMat(n,n);
  for x in ZmodnZ(n) do
    for y in ZmodnZ(n) do
      lperms[Position(e, x)][Position(e, y)] := Position(e, y+1);
      rperms[Position(e, y)][Position(e, x)] := Position(e, x-1);
    od;
  od;

  for x in e do
    Add(s, Position(e, x-1));#Chequear si s(x) es x-1 o x+1!
  od;

  return rec( lperms := lperms, rperms := rperms, size := n, labels := e, s :=  s);
end;

#The inverse solution of a given biquandle
inversebiquandle:=function(biquandle)

local x,y,u,v,n,L,R;

n:=biquandle.size;
L := NullMat(n,n);
R := NullMat(n,n);

for x in [1..n] do
  for y in [1..n] do
    for u in [1..n] do
      for v in [1..n] do
        if r(biquandle, x,y)=[u,v] then
         L[u][v] := x;
         R[v][u] := y;
        fi;
      od;
    od;
  od;
od;

return rec( lperms := L, rperms := R, size := n, labels := biquandle.labels, s := biquandle.s);
end;

# This biquandle is invo(2) \coprod {3} 
aflip := function(n)
  local e, lperms, rperms, N,x, y, ss;
  
  e := Enumerator(ZmodnZ(n));

  lperms := NullMat(Size(e), Size(e));
  rperms := NullMat(Size(e), Size(e));

  for x in e do
    for y in e do
      lperms[Position(e, x)][Position(e, y)] := Position(e, y+x*x*y);
      rperms[Position(e, y)][Position(e, x)] := Position(e, x+x*y*y);
    od;
  od;

ss:=[];
N:=3;

for x in [1..N] do
 for y in [1..N] do
  if 
[lperms[x][y], rperms[y][x] ]=[x,y] then
   Add(ss,y);
  fi;
 od;
od;

  return rec( lperms := lperms, rperms := rperms, size := 3, labels := e, s :=  ss);
end;



#Derived solution: given a (non deg) solution of the YBeq this functions give the associated rack (or quandle)
derived:=function(q)
local lperms,rperms,x,y,g,f,gg,N;

N:=Size(q.rperms);

lperms := NullMat(N,N);
rperms := NullMat(N,N);
gg := NullMat(N,N);
g := q.lperms;
f := q.rperms;

for x in [1..N] do
for y in [1..N] do
gg[x][g[x][y]]:=y;
od;od;

for x in [1..N] do
for y in [1..N] do
lperms[x][y]:=y;
rperms[y][x]:=g[y][ f[gg[x][y]][x] ];
od;od;


return rec(lperms:=lperms,rperms:=rperms);
end;





# for a crosssing given by  c=[\pm 1,[a,b,c,d] = [c[1] , [ c[2][1],c[2][2],c[2][3],c[2][4] ] 
# this function checks if a coloring is well done
check_equation := function(c, set, biquandle)
  if c[1] = 1 then
    if not r(biquandle, set[c[2][1]], set[c[2][4]]) = [set[c[2][2]], set[c[2][3]]] then
      return false;
    else
      return true;
    fi;
  else
    if not r(biquandle, set[c[2][3]], set[c[2][4]]) = [set[c[2][2]], set[c[2][1]]] then
      return false;
    else
      return true;
    fi;
  fi;
end;


### this function computes all biquandle-colorings of a planar diagram (that includes the sign of the crossing) 
colorings := function(pd, biquandle)
  local e, i, crossing, p, candidato, max, coloreos;

  e := [1..biquandle.size];
  i := 0;
  max := Maximum(Flat(pd));
  coloreos := [];
  crossing := true;

  for candidato in Iterator(Tuples(e, max)) do
    if ForAny(pd, x->check_equation(x, candidato, biquandle)=false) then
      continue;
    else
      Add(coloreos, candidato);
    fi;
  od;
  return coloreos;
end;


### this function computes all biquandle-colorings in a more efficient way
colormejor := function(pd, biquandle)
  local  lista, n, c, colorprevio, colorviejo, candidatos, precandidatos, candidatoacumulado2,e, colnuevo, solocolor, ppd, coloreos, semiarcospasados, x, y, i, j, cn, candi_acum;

  semiarcospasados:=[];
  colorprevio:=[];
  colnuevo:=[]; 
  colorviejo:=[];
  candi_acum:=[];
  candidatoacumulado2:=[];
  e := [1..biquandle.size];
 
  #paints the first crossing

  candidatos:=[];
  precandidatos:= [];
  i:=0;
  for x in e do
    for y in e do
      if pd[1][1]=1 then
        precandidatos:=[
[pd[1][2][1],x],
[pd[1][2][2],r(biquandle, x,y)[1] ],
[pd[1][2][3],r(biquandle, x,y)[2] ],
[pd[1][2][4],y]
]; 
      else
        precandidatos:=[
[pd[1][2][1],r(biquandle, x,y)[2] ],
[pd[1][2][2],r(biquandle, x,y)[1] ], 
[pd[1][2][3],x],
[pd[1][2][4],y]
]; 
      fi;
      i:=i+1;
      candidatos[i]:=precandidatos;
      precandidatos:=[];
    od;
  od;
  semiarcospasados:=pd[1][2];
  lista:=candidatos;

  ppd:=Difference(pd,[pd[1]]);

  for c in ppd do;
  #paints every crossing
  j:=0;
  candidatos:=[];
  precandidatos:= [];
    for x in e do
      for y in e do
  #colors every crossing satisfying r
  #without checking if a semiarc appears more than once
        if c[1]=1 then
        precandidatos:=[
[c[2][1],x],
[c[2][2],r(biquandle, x,y)[1] ],
[c[2][3],r(biquandle, x,y)[2] ],
[c[2][4],y]
]; 
      else
        precandidatos:=[
[c[2][1],r(biquandle, x,y)[2] ],
[c[2][2],r(biquandle, x,y)[1] ], 
[c[2][3],x],
[c[2][4],y]
]; 
        fi;
      j:=j+1;
      colnuevo[j]:=precandidatos;
      precandidatos:=[];
      od;
    od;

  #finishes using "colnuevo"

    n:=Size(IntersectionSet(semiarcospasados,c[2]));
    i:=0;
    candidatoacumulado2:=[];
    for colorviejo in lista do
      for cn in colnuevo do
        if Size(IntersectionSet(cn,colorviejo))=n then
          i:=i+1;
          candidatoacumulado2[i]:=UnionSet(cn,colorviejo);
        fi;
      od;
    od;
    lista:=candidatoacumulado2;
#Print(Size(lista));
#Print("\n");
    candidatoacumulado2:=[];
    semiarcospasados:=UnionSet(semiarcospasados,c[2]);

  od;

  #transforms the list [[semicarco,color]..] in [[color],..]
  solocolor:=[];
  coloreos:=[];
for x in lista do
  solocolor:=[];
  for j in [1..Size(semiarcospasados)] do
    for i in [1..Size(semiarcospasados)] do
      if x[i][1]=j then
        solocolor[j]:=x[i][2];
       fi;
    od;
  od;
  Add(coloreos,solocolor);
od;

return(coloreos);

end;


# give the list of colorings in terms of the labels of the biquandle
# (instead of in terms of the numbering of the elements of the biquandle)
colorings_withlabels := function(pd, biquandle)
  local c;
  c := colormejor(pd, biquandle);
  return List(c, x->List(x, y->biquandle.labels[y]));
end;

#add 1 in a pd (jndependent from sign) but for coloring in flip 
#it doesn't mater
add1:=function(pd)
local 1pd,i;
1pd:=[];
for i in [1..Size(pd)] do
Add(1pd,[1,pd[i]]);
od;
return 1pd;
end;

### connected components
componentesconexas := function(pd)
  local n;

if Size(pd[1])=2 then
n:=Log2Int(Size(colormejor(pd,bialexander(2,1,1))));
else
n:=Log2Int(Size(colormejor(add1(pd),bialexander(2,1,1))));
fi;
return n;

end;

###first semiarc of the second component
n1 := function(pd)
  local n, i, A, c;
if Size(pd[1])=2 then
  c:=colormejor(pd,bialexander(2,1,1))[2];
else
  c:=colormejor(add1(pd),bialexander(2,1,1))[2];
fi;

  A:= Maximum(Flat(pd));
  for i in [1..A-1] do
    if c[i]<> c[i+1] then
    n:=i+1;
    fi;
  od;
  return n;
end;


### gives [n2,n3]= [first semiarc of the second component, first semiarc of the third component]
ns := function(pd)
  local n, i, M, c,Ns, spd;
if Size(pd[1])=4 then
spd:=add1(pd);
else
spd:=pd;
fi;
  
  Ns:=[];
for c in colormejor(spd,bialexander(2,1,1)) do;
  M:= Maximum(Flat(spd));
  for i in [1..M-1] do
    if c[i]<> c[i+1] then
    Add(Ns,i+1);
    fi;
  od;
od;
  Ns:=UnionSet(Ns,[]);
  return Ns;
end;



#Componentes conexas de un biquandle
T:=function(Q)

local N,x,y,clases;

N:=oficialbiquandle(Q).size;
clases:=[];

for x in [1..N] do
for y in [1..N] do
clases:=UnionSet(clases,[[x,r(Q,x,y)[2]],[y,r(Q,x,y)[1]]]);
od;
clases:=relationgenerated(clases);
od;
return clases;

end;

#sistem of representatives in a biquandle module
repclase:=function(biquandle)
local x, reps, igualdades;
igualdades:=[];

for x in [1..biquandle.size] do
 Add( igualdades,[x,biquandle.s[x]] );
od;
reps:=relationgenerated(igualdades);
return reps;
end;

# this function has as input some equalites, some set of trivial
# elements S, and a set of cocycle equations. The function
# evaluates the elements of S in the cocycle equaion, identify elements
# that are equal, and get -if possible- new trivial elements
# (comming from equalities or from cocycle condition), and get
# new equalities (from cocycle condition of type ab=ac, ba=ca, or sa=bs with s in S, etc)
#returns an enlarged set S, a new set of equalities, and shorter list of cocycle conditions
New:=function(lista,clases)
local i, l, c, cl, lista2,lista3;
 cl:=clases;
 #find equalities cocycle eq.
 lista2:=[];
 for l in lista do #try to deduce new equalities
   if l[1]=l[3] then
     if l[2]=l[4] then continue;
     else 
      cl:=relationgenerated(    UnionSet(cl,[[l[2],l[4]]]) );#use the new equalitye to generate the new classes
     fi; 
   elif l[2]=l[4] then 
      cl:=relationgenerated(UnionSet(cl,[[l[1],l[3]]]) );
   else
     lista2:=UnionSet(lista2,[l]);
   fi;
 od;
 #write the equations choosing 1 representative for each class
 for l in lista2 do
  for i in [1..4] do
    for c in cl do
     if l[i] in c then l[i]:=c[1]; fi;
    od;
  od;
 od;
 lista2:=UnionSet(lista2,[]);
 #write in normal form whenever 1.a or a.1 in on the equations
 for l in lista2 do
   if l[2]=[] then
   l[2]:=l[1];
   l[1]:=[];
   fi;
   if l[4]=[] then
   l[4]:=l[3];
   l[3]:=[];
   fi;
 od;

lista3:=[];
for l in lista2 do
Add(lista3, l);
lista3:=Difference(lista3,[ [l[3],l[4],l[1],l[2] ]]);
od;
 lista2:=UnionSet(lista3,[]);
return rec( list:=lista2, clases:=cl); 
end;

#this computes a set of generators and relations for Unc. Gives as answer the set
# of trivial elements S, the equivalence classes (where any cocycle takes the same value) 
# and the equations (in terms of representatives of the classes)
# from the 2-cocycle condition.
# As input has a biquandle and a set S (so it can compute te reduced Unc)
# the procedure is to iterate the above function
unc:=function(biquandle,S)
local c, clases, i, j, x, y, z, A, d1 , d2, lista, n, lis, cob, cl, S2, gam1,gam2,usados, News;

lista:=[];
A:=[];

#generate clases using condition 2
clases:=[];
n:=biquandle.size;
for x in [1..n] do
  for y in [1..n] do
    for z in [1..n] do
clases:=
UnionSet(
clases,
 [[
 [y,z],[r(biquandle,x,y)[1],r(biquandle,r(biquandle,x,y)[2],z)[1] ] 
]]);
    od;
  od;
od;
#Print(clases);
clases:=relationgenerated(clases);

# Add [] and S = { (x,s(x))} to the set of classes 
S2:=UnionSet(S,[[]]);
n:=biquandle.size;
for x in [1..n] do
S2:=
UnionSet(S2, [[x,biquandle.s[x] ]]);
od;
clases:=UnionSet(clases,[S2]);
clases:=relationgenerated(clases);


#generates the list of equations of condition 1 (cocycle condition)
n:=biquandle.size;
for x in [1..n] do
 for y in [1..n] do
  for z in [1..n] do
   A[1]:=[x,y];
   A[2]:=[r(biquandle,x,y)[2],z];
   A[3]:=[x,r(biquandle,y,z)[1]];
   d1:=r(biquandle,x,r(biquandle,y,z)[1])[2];
   d2:=r(biquandle,y,z)[2];
   A[4]:=[d1,d2];
   for i in [1..4] do
    for c in clases do
     if A[i] in c then A[i]:=c[1]; fi; #choose representatives
    od;
   od;
   lista:=UnionSet(lista,[A]);    
   A:=[];
  od;
 od;  
od;


c:=Size(clases)+1;
lis:=Size(lista)-1;
for i in [1..60] do #repetimos un procedimiento varias veces
 if c=Size(clases) and lis=Size(lista) then continue;
 else
  c:=Size(clases);
  lis:=Size(lista);
 
  News:=New(lista,clases);
  lista:=News.list;
  clases:=News.clases;
 fi;
od;


#otra vuelta por las dudas
  News:=New(lista,clases);
  lista:=News.list;
  clases:=News.clases;


cob:=[];
usados:=[];

for c in clases do
if [] in c then S2:=c; fi;
od;

cl:=Difference(clases,[S2]);

n:=Size(repclase(biquandle));
for i in [1..Size(cl)] do
x:=cl[i][1][1];
y:=cl[i][1][2];
z:=r(biquandle,x,y)[2]; #sigma2(x,y)
for j in [1..n] do
 if x in repclase(biquandle)[j]  then
  gam1:= repclase(biquandle)[j][1];
 fi;
 if z in repclase(biquandle)[j]  then
  gam2:= repclase(biquandle)[j][1];
 fi;
od;

if gam1=gam2 then 
  continue;
 elif [gam1,gam2] in usados then continue;
 else
  usados:=UnionSet(usados,[[gam1,gam2]]);
  Add(cob,
  [
  #clases[i][1],
  #=
  gam1,cl[i][1],gam2#^-1"
  ]);
 fi;
od;

return rec( S:=S2, clases := cl, equations := lista, cocycle := cob);
end;














imprimir_ecuaciones:=function(biquandle,S)
local i, j, clases, Unc, n, co,lista,l, trivial;

Unc:=unc(biquandle,S);

trivial:=Unc.S;
clases:=Unc.clases;
lista:=Unc.equations;

n:=Size(clases);
if n=1 then  
 Print("Unc tiene ");Print(n);Print(" generador:");
elif n=0 then
 Print("Unc ={1}");
else
 Print("Unc tiene ");Print(n);Print(" generadores:");
fi;
Print("\n");
for i in [1..Size(clases)] do
 Print("f_{"); Print(i);Print("}");
 for j in clases[i] do
  Print("=");
Print("(");Print(j[1]);Print(",");Print(j[2]);Print(")");
 od;
 Print(", ");
od;
Print("\n");


Print("elementos triviales:");
Print("\n");
#Print("1");
for i in trivial do
  Print("=");
#Print("(");Print(i[1]);Print(",");Print(i[2]);Print(")");
Print(i);
od;

Print("\n");


Print("equations: ");
Print("\n");

n:=Size(clases);
for l in lista do
  for i in [1..4] do;
  if i=3 then Print("="); fi;  
  if l[i] in trivial then continue;
    else
      for j in [1..n] do
        if l[i] in clases[j] then
          Print("f_{"); Print(j);Print("}");
        fi;
      od; 
    fi; 
  
  od;
    Print(", ");
od;

Print("\n");
Print("\n");

Print("cobordant conditions to eventually add elements to S:");
Print("\n");
for co in Unc.cocycle do
Print("\\gamma_");
Print(co[1]);
Print(co[2]);
Print("\\gamma_");
Print(co[3]);
Print("^{-1}, ");
od;

Print("\n");

end;


#Computes the Bolzman Weight of a crossing, for a given coloring, a set
# of equivalent classes of pairs (were the cocycle takes the same values) 
# and a group (Fnc) given by generators
BWgen := function(cruce, coloreo, clases,Fnc)
local bwg,i;

bwg:=One(Fnc);
  if cruce[1]=1 then 
    for i in [1..Size(clases)] do
      if [coloreo[cruce[2][1]],coloreo[cruce[2][4]]] in clases[i] then
        bwg := GeneratorsOfGroup(Fnc)[i];
      fi;
    od;

  else 
    for i in [1..Size(clases)] do
      if [coloreo[cruce[2][3]],coloreo[cruce[2][4]]] in clases[i] then
        bwg := Inverse(GeneratorsOfGroup(Fnc)[i]);
      fi;
    od;
  fi;
return bwg;
end;

# over all colorings, this function computes the invariant (of each connected component)
# up to 3 connected components. The group is the one computed by "unc", that have
# an S as input, so it can (nearly) compute the reduced Unc, or also a quotient of Unc
invariantgen:=function(pd,biquandle,S)
  local generators,
inv,max, W, W1,W2, W3, BWg, coloreo,   cruce, Fnc,
n, rk, N1, N2,i, j,k,l;

max := Maximum(Flat(pd));
inv:=[];

#calculates U_nc equations and generators
unc(biquandle,S);
generators:=unc(biquandle,S).clases;
rk:=Size(generators);
Fnc := FreeGroup(rk);

#and given a coloring calculates the invariant

n:=componentesconexas(pd);

if n=1 then
  for coloreo in colormejor(pd, biquandle) do
    W := One(Fnc);
    for i in [1..max] do
      for cruce in pd do
        if cruce[2][1]=i then
          W:=W*BWgen(cruce, coloreo, generators,Fnc);
        fi;
      od;
    od;
    Add(inv,[coloreo,W]);
  od;
  return(inv);
elif n=2 then
  N1:=n1(pd);
  for coloreo in colormejor(pd, biquandle) do
    W1 := One(Fnc);
    W2 := One(Fnc);
    for i in [1..N1-1] do
      for cruce in pd do
        if cruce[2][1]=i then
          W1:=W1*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    for j in [N1..max] do
      for cruce in pd do
        if cruce[2][1]=j then
          W2:=W2*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    Add(inv,[coloreo,[W1,W2]]);
  od;
return(inv);
elif n=3 then
  N1:=ns(pd)[1];
  N2:=ns(pd)[2];
  for coloreo in colormejor(pd, biquandle) do
    W1 := One(Fnc);
    W2 := One(Fnc);
    W3 := One(Fnc);
    for i in [1..N1-1] do
      for cruce in pd do
        if cruce[2][1]=i then
          W1:=W1*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    for j in [N1..N2-1] do
      for cruce in pd do
        if cruce[2][1]=j then
          W2:=W2*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    for k in [N2..max] do
      for cruce in pd do
        if cruce[2][1]=k then
          W3:=W3*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    Add(inv,[coloreo,[W1,W2,W3]]);
  od;
return inv;
fi;


end;

#devuelve true si todas las igualdades de la ec de trenzas son ciertas, devuelve false si no.
checkYB:=function(bi)
local x,y,z, N,V;
V:=true;
N:=Size(bi.rperms);

for x in [1..N] do;
for y in [1..N] do;
for z in [1..N] do;
if [r(bi,r(bi,x,y)[1],r(bi,r(bi,x,y)[2],z)[1])[1],
r(bi,r(bi,x,y)[1],r(bi,r(bi,x,y)[2],z)[1])[2]
,  r(bi,r(bi,x,y)[2],z)[2]]
=
[r(bi,x,r(bi,y,z)[1])[1],
r(bi,r(bi,x,r(bi,y,z)[1])[2],r(bi,y,z)[2])[1],
r(bi,r(bi,x,r(bi,y,z)[1])[2],r(bi,y,z)[2])[2]] then continue;
else #Print("no cumple YBeq!!!! ");
V:=false;
z:=N;x:=N;y:=N;
fi;
od;od;od;
return V;
end;

checks:=function(bQ)
local x,y, N,V,ss;
N:=Size(bQ.rperms);
V:=true;

for x in [1..N] do
ss:=[];
 for y in [1..N] do
  if r(bQ,x,y)=[x,y] then
   Add(ss,[x,y]);
  fi;
 od;
 if Size(ss)=1 then continue;
 else V:=false; x:=N;
 fi;
od;

return V;
end;




#testYB a una lista
testYB :=function(listaq)
local i,N;
N:=Size(listaq[1].lperms[1]);
for i in [1..Size(listaq)] do
 if checkYB(listaq[i],N)=true then continue;
 else 
  Print(i);
  Print(" no cumple YB!");
 fi;
od;
end;


#coloreo
testcolor :=function(pd,listaq)
local i;

for i in [1..Size(listaq)] do
Print("  ");
Print(i);
Print(" c=");
Print(Size(colormejor(pd,oficialbiquandle(listaq[i]))));
#Print("\n");
od;

end;

r2:=function(bi,x,y)
local rr;
rr:= r(bi,r(bi,x,y)[1],r(bi,x,y)[2]);
return rr;
end;

#
r3:=function(bi,x,y)
local rr;
rr:= r(bi,r2(bi,x,y)[1],r2(bi,x,y)[2]);
return rr;
end;

#
r4:=function(bi,x,y)
local rr;
rr:= r(bi,r3(bi,x,y)[1],r3(bi,x,y)[2]);
return rr;
end;

#
r5:=function(bi,x,y)
local rr;
rr:= r(bi,r4(bi,x,y)[1],r4(bi,x,y)[2]);
return rr;
end;

#
r6:=function(bi,x,y)
local rr;
rr:= r(bi,r5(bi,x,y)[1],r5(bi,x,y)[2]);
return rr;
end;



#gives the number of (x,y) such that r^N(x,y)=(x,y), for each 
#member of a list of biquandles
orden:=function(lista,N)

local x,y,n,rr,bolsa,i;

rr:=[r2,r3,r4,r5,r6];

n:=Size(lista[1].rperms);

for i in [1..Size(lista)] do
 bolsa:=[];
 for x in [1..n] do
  for y in [1..n] do
   if rr[N](lista[i],x,y)=[x,y] then
    Add(bolsa,[[x,y]]);
   fi;
  od;
 od;
#if Size(bolsa)=16 then
    Print("bi-");
   Print(i);
   Print(" ");
  Print(Size(bolsa));
   Print(" ");
   Print(" ");
#  Print("\n");
#fi; 
#if Size(bolsa)=16 then 
 Print(" ");
 Print(Size(T(oficialbiquandle(lista[i])))); 
   if Size(T(oficialbiquandle(lista[i])))=1 then Print("AAACA HAY UNO INDESCOMPONIBLE!!"); fi; 
 Print("\n");
#fi;
od;
end;
 
#puntos fijos de la diagonal
diagonal:=function(lista)

local i,j,N,bi,St;
N:=Size(lista[1].rperms);

for j in [1..Size(lista)] do
St:=[];
 for i in [1..N] do
  if oficialbiquandle(lista[j]).s[i]=i then
  Add(St,i);
  fi;
 od;
Print("bi-");Print(j);
Print(" ");Print(Size(St));
Print("\n");
od;
end;


biquandleiso:=function(biquandle)
local x,y,sa,sam,i,L,R,bi,n;

n:=biquandle.size;
sa:=[];
sam:=[];
bi:=[];

sa[1]:=[2,1,3];
sa[2]:=[1,3,2];
sa[3]:=[3,2,1];
sa[4]:=[2,3,1];
sa[5]:=[3,1,2];
sa[6]:=[1,2,3];

sam[1]:=[2,1,3];
sam[2]:=[1,3,2];
sam[3]:=[3,2,1];
sam[4]:=[3,1,2];
sam[5]:=[2,3,1];
sam[6]:=[1,2,3];


for i in [1..6] do
L := NullMat(n,n);
R := NullMat(n,n);
for x in [1..n] do
for y in [1..n] do
R[x][y]:=sam[i][biquandle.rperms[sa[i][x]][sa[i][y]]];
L[x][y]:=sam[i][biquandle.lperms[sa[i][x]][sa[i][y]]];
od;
od;
bi[i]:=oficialbiquandle(rec(rperms:=R,lperms:=L));
od;

return bi;
end;



#le das uno famoso y la lista y lo busca
buscar:=function(biquandle,lista)

local i,j,bi;

for i in [1..Size(lista)] do
for j in [1..6] do;
bi:=biquandleiso(biquandle)[j];

  if bi.rperms= lista[i].rperms then
  if bi.lperms=lista[i].lperms then
Print("\n");
Print("encontrado!!!! es bi3-");
Print(i);
Print(" con el isomorfismo ");
Print(j);
  else continue;
  fi;
else continue;
fi;
od;
od;
end;

testeoinverso:=function(biquandle)

local x,y,n,bi,ibi;
bi:=biquandle;
ibi:=inversebiquandle(bi);
n:=bi.size;

for x in [1..n] do
for y in [1..n] do
if [x,y] =r(ibi,r(bi,x,y)[1],r(bi,x,y)[2]) then continue;
else
Print("ops");
fi;
od;
od;

end;


### This function returns the signed planar diagram (spd) of a KNOT
spd := function(pd)## here pd means a not signed planar diagram
  local c, s, max;
  s := [];
  max := Maximum(Flat(pd));

  for c in pd do
    if c[2] = max and c[4]=max-1 then 
      Add(s,[1,c]);
      continue;
    elif c[2] = max and c[4]=1 then 
      Add(s,[-1,c]);
    continue;
    elif c[1] = max and c[2]>c[4] then
      Add(s, [1,c]);
      continue;
    elif c[1] = max and c[2]<c[4] then
      Add(s, [-1,c]);
      continue;
    elif c[4] = max and c[2]=max-1 then
      Add(s, [-1,c]);
      continue;
    elif c[4] = max and c[2]=1 then
      Add(s, [1,c]);
      continue;
    elif c[4]<c[2]  and c[1]<c[3]then
      Add(s, [1,c]);
      continue;
    elif c[4]>c[2]  and c[1]<c[3]then
      Add(s, [-1,c]);
      continue;
  else
      Print("Nothing?\n");
      return fail;
    fi;
  od;
  return s;
end;


#given a knot, calculates the signed pd (for a given and not signed pd)
signocruce:=function(c)
local  max;
max := Maximum(Flat(c));

     if c[2] = max and c[1]<c[3] then 
      return [-1,c];
      continue;
    elif c[1] = max and c[2]>c[4] then
      return [1,c];
      continue;
    elif c[1] = max and c[2]<c[4] then
      return [-1,c];
      continue;
    elif c[4] = max and c[1]<c[3] then
      return [1,c];
      continue;
    elif c[4]<c[2]  and c[1]<c[3]then
      return [1,c];
      continue;
    elif c[4]>c[2]  and c[1]<c[3]then
      return [-1,c];
    else 
      Print("Huston!\n");
      continue;
    fi;
end;

# calculates the sing of a crossing for a 2-component 
#link knowing N1=N and  max=maximun(Flat(sspd))
signon1:=function(c,N,max)
     if c[1] = max and c[4]=N-1 then 
      return [1,c];
    elif c[1] = N-1 and c[4]=max then
      return [1,c];
    elif c[1] = max and c[2]=N-1 then
      return [-1,c];
    elif c[1] = N-1 and c[2]=max then
    return [-1,c];
    elif c[1]=max  and c[4]+1=c[2]then
      return [1,c];
    elif c[4]=max  and c[1]+1=c[3]then
      return [1,c];
    elif c[1]=N-1  and c[2]=c[4]+1 then
      return [1,c];
    elif c[4]=N-1  and c[1]+1=c[3] then
      return [1,c];
    elif c[1]=max  and c[3]=1+c[2]then
      return [-1,c];
    elif c[2]=max  and c[1]+1=c[3]then
      return [-1,c];
    elif c[1]=N-1  and c[2]+1=c[4]then
      return [-1,c];
    elif c[2]=N-1  and c[1]+1=c[3]then
      return [-1,c];
    elif c[4]+1=c[2]  and c[1]+1=c[3]then
      return [1,c];
    elif c[4]=1+c[2]  and c[1]+1=c[3]then
      return [-1,c];
    else 
      Print("Houston!\n");
      continue;
    fi;
end;

### This function returns the list of signs of the crossing in a pd of a KNOT
signs := function(pd)
  local c, s, max;
  s := [];
  max := Maximum(Flat(pd));

  for c in pd do
    if c[2] = max and c[1]<c[3] then 
      Add(s,-1);
      continue;
    elif c[1] = max and c[2]>c[4] then
      Add(s, 1);
      continue;
    elif c[1] = max and c[2]<c[4] then
      Add(s, -1);
      continue;
    elif c[4] = max and c[1]<c[3] then
      Add(s, 1);
      continue;
    elif c[4]<c[2]  and c[1]<c[3]then
      Add(s, 1);
      continue;
    elif c[4]>c[2]  and c[1]<c[3]then
      Add(s, -1);
      continue;
  else
      Print("Nothing?\n");
      return fail;
    fi;
  od;
  return s;
end;


#calculates n1 for a link of 2-componentes 
fulln1 := function(sspd)
  local c, s, max, N1;
  N1 := [];
  max := Maximum(Flat(sspd));

  for c in sspd do
    if c[3] =1 then 
      N1:=c[1]+1;
#      continue;
    elif c[4] = 1 and c[2]<>2 then
      N1:=c[2]+1;
#      continue;
    elif c[2] = 1 and 2<>c[4] then
      N1:=c[4]+1;
#      continue;
#    elif c[2]=2 and c[4]=1 and c[1]=3 then  
#      Print("hopf link?");
#     continue;
#    elif c[2]=2 and c[4]=1 and c[3]=3 then 
#      Print("hopf link?");
#     continue;
#    elif c[2]=1 and c[4]=2 and c[1]=3 then  
#      Print("hopf link?");
#     continue;
#    elif c[2]=1 and c[4]=2 and c[3]=3 then 
#      Print("hopf link?");
#     continue;
#     continue;
    else# Print("hofp?");
    fi;
  od;
#if N1=[] then 
#  Print("hopf??????");
#fi;
#Print("N1:");
  return N1;
end;