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, with labels, s, size
inversebiquandle:=function(biquandle)

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

n:=Size(biquandle.rperms);
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;

#The inverse solution of a given biquandle, withOUT labels, s, size
inversebiquandleSMALL:=function(biquandle)

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

n:=Size(biquandle.rperms);
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);
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;



check_virtual := function(c, set, bq,bi)
  if c[1] = 1 then
    if not r(bq, set[c[2][1]], set[c[2][4]]) = [set[c[2][2]], set[c[2][3]]] then
      return false;
    else
      return true;
    fi;
  fi;
  if c[1]=-1 then
    if not r(bq, set[c[2][3]], set[c[2][4]]) = [set[c[2][2]], set[c[2][1]]] then
      return false;
    else
      return true;
    fi;
  fi;
  if c[1]=0 then
    if not r(bi, set[c[2][1]], set[c[2][4]]) = [set[c[2][2]], set[c[2][3]]] 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) 
colorvirtual:= function(pd, bq, bi)
  local n, e, crossing, p, candidato, max, coloreos;

n:=Size(bq.rperms);
e:=[1..n];
max := Maximum(Flat(pd));
coloreos := [];

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

###################################
###################################


### this function computes all biquandle-colorings in a more efficient way
colormejorvirtual := function(pd, bq,bi)
  local  lista,N, 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:=[];
  N:= Size(bq.rperms);
  e:=[1..N];
 
  #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(bq, x,y)[1] ],
[pd[1][2][3],r(bq, x,y)[2] ],
[pd[1][2][4],y]
]; 
      fi;
      if pd[1][1]=-1 then
        precandidatos:=[
[pd[1][2][1],r(bq, x,y)[2] ],
[pd[1][2][2],r(bq, x,y)[1] ], 
[pd[1][2][3],x],
[pd[1][2][4],y]
]; 
      fi;
      if pd[1][1]=0 then
        precandidatos:=[
[pd[1][2][1],x],
[pd[1][2][2],r(bi, x,y)[1] ],
[pd[1][2][3],r(bi, x,y)[2] ],
[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(bq, x,y)[1] ],
[c[2][3],r(bq, x,y)[2] ],
[c[2][4],y]
];     
      fi;
        if c[1]=-1 then
        precandidatos:=[
[c[2][1],r(bq, x,y)[2] ],
[c[2][2],r(bq, x,y)[1] ], 
[c[2][3],x],
[c[2][4],y]
]; 
        fi;
        if c[1]=0 then
        precandidatos:=[
[c[2][1],x],
[c[2][2],r(bi, x,y)[1] ],
[c[2][3],r(bi, x,y)[2] ],
[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;


### this function computes all biquandle-colorings in a more efficient way
colormejorsingular := function(pd, bq,tau)
  local  lista,N, 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:=[];
  N:= Size(bq.rperms);
  e:=[1..N];
 
  #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(bq, x,y)[1] ],
[pd[1][2][3],r(bq, x,y)[2] ],
[pd[1][2][4],y]
]; 
      fi;
      if pd[1][1]=-1 then
        precandidatos:=[
[pd[1][2][1],r(bq, x,y)[2] ],
[pd[1][2][2],r(bq, x,y)[1] ], 
[pd[1][2][3],x],
[pd[1][2][4],y]
]; 
      fi;
      if pd[1][1]=2 then
        precandidatos:=[
[pd[1][2][1],x],
[pd[1][2][2],r(tau, x,y)[1] ],
[pd[1][2][3],r(tau, x,y)[2] ],
[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(bq, x,y)[1] ],
[c[2][3],r(bq, x,y)[2] ],
[c[2][4],y]
];     
      fi;
        if c[1]=-1 then
        precandidatos:=[
[c[2][1],r(bq, x,y)[2] ],
[c[2][2],r(bq, x,y)[1] ], 
[c[2][3],x],
[c[2][4],y]
]; 
        fi;
        if c[1]=2 then
        precandidatos:=[
[c[2][1],x],
[c[2][2],r(tau, x,y)[1] ],
[c[2][3],r(tau, x,y)[2] ],
[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;







#connected components of a 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, re, igualdades;
igualdades:=[];

for x in [1..biquandle.size] do
 Add( igualdades,[x,biquandle.s[x]] );
od;
re:=relationgenerated(igualdades);
return re;
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;





#
#Uncfg
#
#
#
#sistema de representantes de biQ/s
clasesS:=function(biquandle)
local c,eq,x,n, igualdades;
 igualdades:=[];
n:=biquandle.size;

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


reprs:=function(S,u)
local c;
for c in S do
if u in c then return c[1];
fi;od;
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 the reduced Unc)
gen_eq_f:=function(q)
local A,d1,d2,n,x,E,y,z,u,reps,c,clases ,s,Equations,w,S;

n:=Size(q.rperms);
S:=[];
for x in [1..n] do
for y in [1..n] do
Add(S,["f",[x,y]]);
od;od;

s:=oficialbiquandle(q).s;
E:=[];
for x in [1..n] do

#[f4)
u:=[
[ ["f" ,[x,s[x]]] , [ ]]
, [ [] , [] ] 
];
 Add(E,ShallowCopy(u)); 
od;

E:=Set(E);


clases:=[];
#las igualdades de la condicion 2
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do
 w:=Set([ ["f",[y,z] ] , ["f",[r(q,x,y)[1],r(q,r(q,x,y)[2],z)[1] ]] ]);
   Add(clases,w);
od;od;od;

clases:=relationgenerated(clases);

reps:=[];
for c in clases do
Add(reps, c[1]);
od;


Equations:=E;
A:=[1,2,3,4];
#ecuaciones de cociclo
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do
  A[1]:=reprs(clases, ["f",[x,y]]);
  A[2]:=reprs(clases, ["f",[r(q,x,y)[2],z]]);
  A[3]:=reprs(clases, ["f",[x,r(q,y,z)[1]]]);
  d1:=r(q,x,r(q,y,z)[1])[2];
  d2:=r(q,y,z)[2];
  A[4]:=reprs(clases, ["f",[d1,d2]]);
  
  w:=Set([ [A[1],A[2]],[A[3],A[4]] ]);
  if Size(w)=1 then
   continue;
  else
   Add(Equations,w);
  fi;
od;
od;
od;

Equations:=Set(Equations);


for u in reps do
Equations:=Difference(Equations, [ [  [u, u],[u ,u] ] ]);
Equations:=Difference(Equations, [ [  [u,[]],[[],u] ] ]);
Equations:=Difference(Equations, [ [  [u,[]],[u,[]] ] ]);
Equations:=Difference(Equations, [ [  [[],u],[[],u] ] ]);
Equations:=Difference(Equations, [ [  [[],u],[u,[]] ] ]);
od;

return rec(generators:=clases,equations:=Equations);
end;





rep:=function(S,p)
local c;
if p=[] then return[]; fi;

 for c in S do
  if p in c then 
   return c[1];
  fi;
 od;
 return fail;
end;


##
#find using inverses
LookWithInverses:=function(equations)
local eq,a,b,c,d,I,C,A,g;
I:=[];


for eq in equations do
a:=eq[1][1];  b:=eq[1][2];  c:=eq[2][1];  d:=eq[2][2];
if a=b and a=[] then Add(I,[c,[d,-1]]);
fi;
od;

C:=relationgenerated(I);
I:=[];
for c in C do 
A:=[];
 for g in c do
 if g=[] then continue;
 else
 if g[2]=-1 then continue; else Add(A,g);  fi;
 fi;
 Add(I,A);
 od;
od;

return I;
end;




##
recycle_eq:=function(S,Eq)
local eq,a,aa,SS,c,reps,u,EE,i,b,d,x,p1,p2;
SS:=S;
for eq in Eq do
a:=eq[1][1];b:=eq[1][2];c:=eq[2][1];d:=eq[2][2];
if a=c then Add(SS,[b,d]);fi;
if b=d then Add(SS,[a,c]);fi;

if a=[] and b=c then Add(SS,[[],d]);fi;
if b=[] and a=d then Add(SS,[[],c]);fi;
if c=[] and a=d then Add(SS,[[],b]);fi;
if d=[] and b=c then Add(SS,[[],a]);fi;

if a=[] and d=[] then Add(SS,[b,c]);fi;
if b=[] and c=[] then Add(SS,[a,b]);fi;
od;
SS:=relationgenerated(SS);

EE:=[];
for eq in Eq do
a:=eq[1][1];  b:=eq[1][2];  c:=eq[2][1];  d:=eq[2][2];
a:=rep(SS,a); b:=rep(SS,b); c:=rep(SS,c); d:=rep(SS,d);
if b=[] then p1:=[b,a]; else p1:=[a,b];fi;
if d=[] then p2	:=[d,c]; else p2:=[c,d];fi;
x:=Set([p1,p2]);
if Size(x)=2 then Add(EE,x); fi;
od;
EE:=Set(EE);

for c in SS do
u:=c[1];
EE:=Difference(EE,[   [[u,u ],[u,u]]  ]);
EE:=Difference(EE,[   [[u,[]],[[],u]] ]);
EE:=Difference(EE,[   [[u,[]],[u,[]]] ]);
EE:=Difference(EE,[   [[[],u],[[],u]] ]);
EE:=Difference(EE,[   [[[],u],[u,[]]] ]);
od;


SS:=relationgenerated(UnionSet(SS,LookWithInverses(EE)));


return  rec(generators:=SS,equations:=EE);
end;





uncf:=function(q)
local geneq,geneq2,k,i,E,S;
geneq:=gen_eq_f(q);
k:=Size(geneq.generators);

for i in [1..k] do
#Print(i,"*");
S:=geneq.generators;
E:=geneq.equations;
geneq2:=recycle_eq(S,E);
if geneq2=geneq then return geneq;fi;
geneq:=geneq2;
od;
return geneq;
end;




checkPair:=function(Q,qinv)
local V,x,y,z,A,B,C,D,n;
n:=Size(Q.rperms);
V:=true;
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do
A:=r(qinv,x,y)[1];
B:=r(qinv,x,y)[2];
C:=r(qinv,y,z)[1];
D:=r(qinv,y,z)[2];
if 
[
r(qinv, r(qinv,x,y)[1],r(Q,r(qinv,x,y)[2]   ,z)[1]  )[1],
r(qinv, r(qinv,x,y)[1],r(Q,r(qinv,x,y)[2]   ,z)[1]  )[2],
 r(Q,r(qinv,x,y)[2],z)[2]
]
<>
[
r(Q,x,r(qinv,y,z)[1])[1],
r(qinv, r(Q,x,r(qinv,y,z)[1] )[2]  , r(qinv,y,z)[2])[1],
r(qinv, r(Q,x,r(qinv,y,z)[1] )[2]  , r(qinv,y,z)[2])[2]
]

then return false; fi;
od;od;od;
return true;
end;
#
checkcycle:=function(sigma,Q,q,n)
local i;

for i in [1..n] do 
if
 CycleStructurePerm(PermList(q.rperms[i]))
<>
 CycleStructurePerm(PermList(Q.rperms[i^sigma]))
or
 CycleStructurePerm(PermList(q.lperms[i]))
<>
 CycleStructurePerm(PermList(Q.lperms[i^sigma]))
then
 return false;
fi;
od;
return true;
end;


#
qhomcycle:=function(q,Q)
local n,Sn,G,sigma;
n:=Size(q.rperms);
if Size(Q.rperms)<>n then return false; fi;

Sn:=SymmetricGroup(n);
G:=[];
for sigma in Sn do
if checkcycle(sigma,Q,q,n)=false then continue;
else
Add(G,sigma);
fi;
od;
return G;
end;

#
checkmorphism:=function(p,bq)
local x,y,n;
n:=Size(bq.rperms);
 for x in [1..n] do for y in [1..n] do
  if r(bq,x,y)[1]^p <> r(bq,x^p,y^p)[1] or r(bq,x,y)[2]^p <> r(bq,x^p,y^p)[2]
  then return false;
  fi;
 od;od;
return true;
end;



#compute the list of isomorphisms of a biquandle
qmorphismCycle:=function(bq)
local p,G,iso,n;
n:=Size(bq.rperms);
iso:=[];
G:=qhomcycle(bq,bq);

for p in G do
if checkmorphism(p,bq) then
Add(iso,p);
fi;
od;

return iso;
end;




#compute the list of isomorphisms of a biquandle, without using the cycle-preserving condition first
qmorphism:=function(bq)
local p,G,iso,n;
n:=Size(bq.rperms);
iso:=[];
#G:=qhomcycle(bq,bq);
G:=SymmetricGroup(n);

for p in G do
if checkmorphism(p,bq) then
Add(iso,p);
fi;
od;

return iso;
end;





#
listaisoCycle:=function(biquandle)
local x,y,i,f,ff,lista,L,R,iso,n,u,U,U1,orb;

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

U:=qmorphismCycle(biquandle);
if Size(U)=1 then 
 iso:=AsList(SymmetricGroup(n));
else
 orb:=OrbitsDomain(AsGroup(U), AsList(SymmetricGroup(n)), OnRight);
 iso:=List(orb,c->Representative(c));
fi;

for f in iso do

 L := NullMat(n,n);
 R := NullMat(n,n);

 for x in [1..n] do
 for y in [1..n] do
 R[x][y]:=biquandle.rperms[x^f][y^f]^Inverse(f);
 L[x][y]:=biquandle.lperms[x^f][y^f]^Inverse(f);
 od;
 od;
 Add(lista, rec(rperms:=ShallowCopy(R),lperms:=ShallowCopy(L)) );
od;

return lista;
end;


#
pares:=function(lb,li)
local qb,qi,qii,A,B,c;
B:=[];
for qb in lb do
A:=[];
Print("*");
for qi in li do
Print("|");
for qii in  listaisoCycle(qi) do
if checkPair(qb,qii)=true then Add(A, qii);Print(".");
fi;
od;
od;
Add(B,[qb,A]);
od;
return B;
end;



listaisoG:=function(G,qi)
local x,y,i,f,ff,lista,L,R,iso,n,u,U,U1,orb;

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

for f in G do

 L := NullMat(n,n);
 R := NullMat(n,n);

 for x in [1..n] do
 for y in [1..n] do
 R[x][y]:=qi.rperms[x^f][y^f]^Inverse(f);
 L[x][y]:=qi.lperms[x^f][y^f]^Inverse(f);
 od;
 od;
 Add(lista, rec(rperms:=ShallowCopy(R),lperms:=ShallowCopy(L)) );
od;

return lista;
end;


##


##
#

#
isopar:=function(q,li)
local G,l1,l2,qq;

G:=qmorphismCycle(q);
l1:=[];
l2:=li;

while l2<>[] do
qq:=l2[1];
Add(l1,qq);
l2:=Difference(l2,[qq]);
l2:=Difference(l2,listaisoG(G,qq));
od;
return l1;
end;

paresiso:=function(lb,li)
local P,p,B;
B:=[];
Print("biquadnles:",Size(lb)," invo:",Size(li),"\n");
P:=pares(lb,li); 
for p in P do 
Add(B,[p[1],isopar(p[1],p[2])]);
od;
return B;
end;





simple_equalities_g:=function(pp) #pp=[q,invo]
local n,x,y,z,bqi,S,SS,q,qi,u;
q:=pp[1];
qi:=pp[2];
n:=Size(qi.rperms);
S:=[];
bqi:=oficialbiquandle(qi);

for x in [1..n] do
Add(S,[ [] , ["g", [x , bqi.s[x] ] ] ]);#g1
od;

S:=relationgenerated(S);

for x in [1..n] do
for y in [1..n] do
for z in [1..n] do
u:=Set( [["g",[y,z]],["g", [r(q,x,y)[1],r(qi,r(q,x,y)[2],z)[1]]] ]);#m1
if Size(u)=1 then 
 continue;
else
 Add(S,ShallowCopy(u));
fi;
u:=Set( [["g",[y,z]],["g",[r(qi,x,y)[1],r(q,r(qi,x,y)[2],z)[1]]] ]);#m2
if Size(u)=1 then continue;
else
  Add(S,ShallowCopy(u));
fi;
od;od;od;

for x in [1..n] do
for y in [1..n] do
Add(S,[["g",[x,y]],["g",[x,y ]]]);
od;od;
#Print(S,"\n","\n");
SS:=relationgenerated(S);

return SS;
end;



double_equalities_g:=function(pair)
local n,x,y,z,binvo,S,q,invo,E,eq, reps,u;
q:=pair[1];
invo:=pair[2];
n:=Size(invo.rperms);
S:=simple_equalities_g(pair);

reps:=[];
for u in S do
Add(reps,u[1]);
od;

E:=[];
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do
Add(E, [rep(S,[x,y]),rep(S,r(invo,x,y)),[],[]]); #g2
Add(E,[ rep(S,[x,y]),rep(S,[r(invo,x,y)[2],z]),rep(S,[x,r(invo,y,z)[1]]),rep(S,[r(invo,x,r(invo,y,z)[1])[2],r(invo,y ,z )[2]])]);#g3
Add(E,[ rep(S,[y,z]),rep(S,[r(invo,x,r(invo,y,z)[1])[2],r(invo,y,z)[2]]),rep(S,[x,y]),rep(S,[r(invo,x,y)[1],r(invo,r(invo,x,y)[2] ,z )[1]])]);#g4
Add(E,[ rep(S,[y,z]),rep(S,[x,r(invo,y,z)[1]]),rep(S,[r(invo,x,y)[2],z]),rep(S,[r(invo,x,y)[1],r(invo,r(invo,x,y)[2],z)[1]])]);#g5
Add(E,[ rep(S,[x,y]),rep(S,[r(invo,x,y)[2],z]),rep(S,[x,r(q,y,z)[1]]),rep(S,[r(invo,x,r(q,y,z)[1])[2],r(q,y,z)[2]])]);#m3
Add(E,[ rep(S,[y,z]),rep(S,[r(q,x,r(invo,y,z)[1])[2] ,r(invo,y,z)[2]]),rep(S,[x,y]),rep(S,[r(invo,x,y)[1],r(q,r(invo,x,y)[2],z)[1]])]);#m4
od;od;od;

for u in reps do
E:=Difference(E, [  [u,u,u,u]  ]);
E:=Difference(E, [  [u,[],[],u]  ]);
E:=Difference(E, [  [u,[],u,[]]  ]);
E:=Difference(E, [  [[],u,[],u]  ]);
E:=Difference(E, [  [[],u,u,[]]  ]);
od;

for eq in E do
E:=Difference(E,[[eq[3],eq[4],eq[1],eq[2]]]);
E:=UnionSet(E,[eq]);
od;

return rec(generators:=S,equations:=E);
end;
#

gen_eq_f_g:=function(pair)
local n,x,y,z,binvo,S,q,invo,E,eq, reps,u;
q:=pair[1];
invo:=pair[2];
n:=Size(invo.rperms);
S:=simple_equalities_g(pair);

reps:=[];
for u in S do
Add(reps,u[1]);
od;

E:=[];
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do

#[m5) r(invo,
u:=Set([
[ ["g" ,[x,r(invo,y,z)[1]]] , ["f" ,[r(invo,x,r(invo,y,z)[1])[2],r(invo,y,z)[2]] ] ],
[ ["f",[x,y]] , ["g",[ r(q,x,y)[2],z]] ]
]);
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
#[m6)
u:=Set([
[
["g",[x,y]],["f",[r(invo,x,y)[1],r(invo,r(invo,x,y)[2],z)[1] ]] 
],
[
["f",[y,z]],["g",[r(invo,x,r(q,y,z)[1])[2],r(q,y,z)[2]	 ]]
]
]);
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
od;od;od;


E:=Set(E);

for u in reps do
E:=Difference(E,[   [[u,u ],[u,u]]  ]);
E:=Difference(E,[   [[u,[]],[[],u]] ]);
E:=Difference(E,[   [[u,[]],[u,[]]] ]);
E:=Difference(E,[   [[[],u],[[],u]] ]);
E:=Difference(E,[   [[[],u],[u,[]]] ]);
od;

return rec(generators:=S,equations:=E);
end;
#


#
gen_eq_f4:=function(pair)
local n,x,y,z,tau,S,q,E,eq, reps,u;
q:=pair[1];
tau:=pair[2];
n:=Size(q.rperms);
S:=[];
for x in [1..n] do
for y in [1..n] do
Add(S,[["f",[x,y]]]);
od;od;

E:=[];
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do

#[f4)
u:=[
[  ["f" ,[x,y]] , ["f",[r(q,x,y)[2],z]] ]
,
[  ["f",[x,r(tau,y,z)[1]] ] ,
 ["f", [r(q,x,r(tau,y,z)[1] )[2],r(tau,y,z)[2] ] ]]
];
Add(E,ShallowCopy(u)); 
od;od;od;
E:=Set(E);

return rec(generators:=S,equations:=E);
end;
#






#
gen_eq_h:=function(q)
local n,x,y,z,S,E,eq, reps,u;
n:=Size(q.rperms);
S:=[];
for x in [1..n] do
for y in [1..n] do
Add(S,[["h",[x,y]]]);
od;od;

E:=[];
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do

#[h1)
u:=[ 
[  ["h" ,[y,z]] , []]
,
[ ["h",[ r(q,x,y)[1], r(q,r(q,x,y)[2] ,z)[1]   ]] , [] ]
];
Add(E,ShallowCopy(u)); 
od;od;od;

E:=Set(E);

return rec(generators:=S,equations:=E);
end;
#



gen_eq_c:=function(pair)

local n,x,y,z,tau,S,q,E,eq, reps,u;
q:=pair[1];
tau:=pair[2];
n:=Size(q.rperms);
S:=[];
for x in [1..n] do
for y in [1..n] do
Add(S,[["f",[x,y]]]);
Add(S,[["h",[x,y]]]);
od;od;

E:=[];
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do

#c1)
u:=[
[ ["f",[x, r(q,y,z)[1]]] , ["h" , [r(q,x,r(q,y,z)[1])[2] ,r(q,y,z)[2]] ]],
[ ["h",[x,y]] , ["f",[r(tau,x,y)[2],z ] ]]];
Add(E,ShallowCopy(u)); 

#c2)
u:=[
[ ["f",[y,z]] , ["h" , [r(q,x,r(q,y,z)[1])[2] ,r(q,y,z)[2]] ]],
[ ["h",[x,y]] , ["f",[r(tau,x,y)[1],r(q,r(tau,x,y)[2],z)[1] ] ]
]];
Add(E,ShallowCopy(u)); 

#c3)
u:=[
[ ["h",[x,y]] , [] ],
[ ["f",[x,y]] , ["h",r(q,x,y) ] ]
];
Add(E,ShallowCopy(u)); 

#c4)
u:=[
[ ["h",r(q,x,y)] , [] ],
[ ["h",[x,y]] , ["f",r(tau,x,y) ] ]
];
Add(E,ShallowCopy(u)); 
od;od;od;

E:=Set(E);


return rec(generators:=S,equations:=E);
end;
#





gen_eq_g:=function(pair)
local n,x,y,z,binvo,S,q,invo,E,eq, reps,u;
q:=pair[1];
invo:=pair[2];
n:=Size(invo.rperms);
S:=simple_equalities_g(pair);

reps:=[];
for u in S do
Add(reps,u[1]);
od;

E:=[];
for x in [1..n] do
for y in [1..n] do
for z in [1..n] do

u:= Set([ [rep(S,["g",[x,y]]),rep(S,["g",r(invo,x,y)])]   ,   [[],[]] ]); #g2
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;


u:=Set([
[rep(S,["g",[x,y]]),rep(S,["g",[r(invo,x,y)[2],z]])],
[rep(S,["g",[x,r(invo,y,z)[1]]]),rep(S,["g",[r(invo,x,r(invo,y,z)[1])[2],r(invo,y ,z )[2]]])]]);#g3
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
u:=Set([[rep(S,["g",[y,z]]),rep(S,["g",[ r(invo,x,r(invo,y,z)[1])[2], r(invo,y,z)[2] ] ]  )],
[rep(S,["g",[x,y]]),rep(S,["g",[r(invo,x,y)[1],r(invo,r(invo,x,y)[2] ,z )[1] ] ])]]);#g4
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
u:=Set([[rep(S,["g",[y,z]]),rep(S,["g",[x,r(invo,y,z)[1]]])],
[rep(S,["g",[r(invo,x,y)[2],z]]),rep(S,["g",[r(invo,x,y)[1],r(invo,r(invo,x,y)[2],z)[1]]])]]);#g5
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
u:=Set([[rep(S,["g",[x,y]]),rep(S,["g",[r(invo,x,y)[2],z]])],
[rep(S,["g",[x,r(q,y,z)[1]]]),rep(S,["g",[r(invo,x,r(q,y,z)[1])[2],r(q,y,z)[2]]])]]);#m3
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
u:=Set([[rep(S,["g",[y,z]]),rep(S,["g",[r(q,x,r(invo,y,z)[1])[2] ,r(invo,y,z)[2]]])],
[rep(S,["g",[x,y]]),rep(S,["g",[r(invo,x,y)[1],r(q,r(invo,x,y)[2],z)[1]]])]]);#m4
if Size(u)=1 then continue; else Add(E,ShallowCopy(u)); fi;
od;od;od;


E:=Set(E);

for u in reps do
E:=Difference(E,[   [[u,u ],[u,u]]  ]);
E:=Difference(E,[   [[u,[]],[[],u]] ]);
E:=Difference(E,[   [[u,[]],[u,[]]] ]);
E:=Difference(E,[   [[[],u],[[],u]] ]);
E:=Difference(E,[   [[[],u],[u,[]]] ]);
od;

return rec(generators:=S,equations:=E);
end;
#





#REVISAR QUE DA FAIL!!!!!!!
uncfg:=function(pair)
local gg,gf,gfg,k,i,Eq,S,Un,Un2;

gg:=gen_eq_g(pair);
gfg:=gen_eq_f_g(pair);
gf:=gen_eq_f(pair[1]);

S:=relationgenerated(UnionSet(gfg.generators, UnionSet(gg.generators,gf.generators)));

k:=Size(S);
Eq:=UnionSet(gfg.equations,UnionSet(gf.equations, gg.equations));

Un:=rec(generators:=S,equations:=Eq);

for i in [1..k] do
#Print(i,"\n");
Un2:=recycle_eq(S,Eq);
if Un2=Un then return Un;fi;
Un:=Un2;
S:=Un.generators;
Eq:=Un.equations;
od;
return Un;
end;


uncfh:=function(pair)#TERMINAR LOS GEN_EQ_F/FSTAU/H!!!!!!!!!!!!!!!!!!!!!
local gf123,gf4,gh,gc,k,i,Eq,S,Un,Un2,x,y,n;
n:=Size(pair[1].rperms);

gf123:=gen_eq_f(pair[1]);
gf4:=gen_eq_f4(pair);
gh:=gen_eq_h(pair[1]);
gc:=gen_eq_c(pair);

S:=[[[]]];
for x in [1..n] do
for y in [1..n] do
Add(S,[["f",[x,y]]]);
Add(S,[["h",[x,y]]]);
od;od;

k:=Size(S);
Eq:=UnionSet(gf123.equations, gf4.equations);
Eq:=UnionSet(gh.equations, Eq);
Eq:=UnionSet(gc.equations, Eq);

Un:=rec(generators:=S,equations:=Eq);

for i in [1..k] do
#Print(i,"\n Un=",Un,"\n\n\n");
Un2:=recycle_eq(S,Eq);
if Un2=Un then return Un;fi;
Un:=Un2;
S:=Un.generators;
Eq:=Un.equations;
od;
return Un;
end;





let:=function(i)
if i=0 then  return "1" ;fi;
if i=1 then  return "a" ;fi;
if i=2 then  return "b" ;fi;
if i=3 then  return "c" ;fi;
if i=4 then  return "d" ;fi;
if i=5 then  return "e" ;fi;
if i=6 then  return "f" ;fi;
if i=7 then  return "g" ;fi;
if i=8 then  return "h" ;fi;
if i=9 then  return "i" ;fi;
if i=10 then  return "j" ;fi;
if i=11 then  return "k" ;fi;
if i=12 then  return "l" ;fi;
if i=13 then  return "m" ;fi;
if i=14 then  return "n" ;fi;
if i=15 then  return "o" ;fi;
if i=16 then  return "p" ;fi;
if i=17 then  return "q" ;fi;
if i=18 then  return "r" ;fi;
if i=19 then  return "s" ;fi;
if i=20 then  return "t" ;fi;
if (i in [1..20])=false then return i; fi;
end;


ff:=function(S,u)
local i;

if u in S[1] then Print();fi;
for i in [2..Size(S)] do
if u in S[i] then Print(let(i-1)); fi;
od;
end;

#
print_universal_cocyclepair:=function(pair)
local Un,S,Eq,c,k,i,eq,j;
Un:=uncfg(pair);
S:=Un.generators;
Eq:=Un.equations;

Print(Size(S)-1," generadores, ",Size(Eq)," ecuaciones.","\n","\n", "clases:","\n");

j:=0;
for c in S do 
k:=Size(c)-1;
Print(let(j),"=");
j:=j+1;
for i in [1..k] do
Print(c[i],"=");
od;
Print(c[k+1],"\n","\n");
od;
Print("\n","ecuaciones:","\n");
for eq in Eq do
ff(S,eq[1][1]); ff(S,eq[1][2]); Print("="); ff(S,eq[2][1]); ff(S,eq[2][2]); Print("\n");
od;
end;


#
print_gen_eq:=function(Un)
local S,Eq,c,k,i,eq,j;
S:=Un.generators;
Eq:=Un.equations;

Print(Size(S)-1," generadores, ",Size(Eq)," ecuaciones.","\n","\n", "clases:","\n");

j:=0;
for c in S do 
k:=Size(c)-1;
Print(let(j),"=");
j:=j+1;
for i in [1..k] do
Print(c[i],"=");
od;
Print(c[k+1],"\n");
od;
Print("ecuaciones:","\n");
for eq in Eq do
ff(S,eq[1][1]); ff(S,eq[1][2]); Print("="); ff(S,eq[2][1]); ff(S,eq[2][2]); Print("\n");
od;
end;

#
print_Ab:=function(Ab)
local S,Eq,c,k,i,eq,j,t;
S:=Ab.generators;
Eq:=Ab.equations;

Print(Size(S)-1," generadores, ",Size(Eq)," ecuaciones.","\n","\n", "clases:","\n");

j:=0;
for c in S do 
k:=Size(c)-1;
Print("f",j,"=");
j:=j+1;
for i in [1..k] do
Print(c[i],"=");
od;
Print(c[k+1],"\n");
od;
Print("ecuaciones:","\n");
for eq in Eq do

for t in eq[1] do 
if t=[] then Print(" ");
else
Print(t[1],t[2]);
fi;
od;
Print("=");
for t in eq[2] do 
Print(t[1],t[2]);
od;
Print("\n");

od;
end;



#
print_universal_singular_cocyclepair:=function(pair)
local Un,S,Eq,c,k,i,eq,j;
Un:=uncfh(pair);
S:=Un.generators;
Eq:=Un.equations;

Print(Size(S)-1," generadores, ",Size(Eq)," ecuaciones.","\n","\n", "clases:","\n");

j:=0;
for c in S do 
k:=Size(c)-1;
Print(let(j),"=");
j:=j+1;
for i in [1..k] do
Print(c[i],"=");
od;
Print(c[k+1],"\n","\n");
od;
Print("\n","ecuaciones:","\n");
for eq in Eq do
ff(S,eq[1][1]); ff(S,eq[1][2]); Print("="); ff(S,eq[2][1]); ff(S,eq[2][2]); Print("\n");
od;
end;

#
isgtrivial:=function(p)
local n,gen,x,y ;
n:=Size(p[1].rperms);
gen:=uncfg(p).generators[1];
for x in [1..n] do
for y in [1..n] do
if (["g",[x,y]] in gen)=false then 
Print("g(",x,",",y,") is maybe non trivial ");
return false;
fi;
od;od;
return true;
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;


invariantVirtual:=function(pd,biquandle,involutive)
local generators,
inv,max, W, W1,W2, W3, BWg, coloreo,  Un, cruce, Fnc,
n, rk, N1, N2,i, j,k,l;

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

#calculates U_nc equations and generators
Un:=uncfg([biquandle,involutive]);
generators:=Un.generators;
generators:=Difference(generators,[generators[1]]); #le quitamos a [] como generador
rk:=Size(generators);
Fnc := FreeGroup(rk);

#and given a coloring calculates the invariant

n:=componentesconexas(pd);
#n:=1;

if n=1 then
  for coloreo in colormejorvirtual(pd, biquandle,involutive) do
    W := One(Fnc);
    for i in [1..max] do
      for cruce in pd do
        if (cruce[2][1]=i) or (cruce[1]=0 and cruce[2][4]=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)-1;
  for coloreo in colormejorvirtual(pd, biquandle,involutive) do
    W1 := One(Fnc);
    for i in [1..N1] do
      for cruce in pd do
        if (cruce[2][1]=i) or (cruce[1]=0 and cruce[2][4]=i) then
          W1:=W1*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    W2 := One(Fnc);
    for j in [N1+1..max] do
      for cruce in pd do
        if (cruce[2][1]=j) or (cruce[1]=0 and cruce[2][4]=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 colormejorvirtual(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;



linkingVirtual:=function(pd)
local generators, inv,max, W, W1,W2, W3, BWg, coloreo,  Un, cruce, Fnc,
n, rk, N1, N2,i, j,k,l, biquandle,involutive;

biquandle:=bialexander(2,1,1);
involutive:=bialexander(2,1,1);

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

#calculates U_nc equations and generators
Un:=uncfg([biquandle,involutive]);
generators:=Un.generators;
generators:=Difference(generators,[generators[1]]); #le quitamos a [] como generador
rk:=Size(generators); # rk=4
Fnc := FreeGroup("a","b","h","H");

#and given a coloring calculates the invariant

n:=componentesconexas(pd);
n:=2;
#OJO EL CASO N=3!!!!)

if n=1 then
  for coloreo in colormejorvirtual(pd, biquandle,involutive) do
    W := One(Fnc);
    for i in [1..max] do
      for cruce in pd do
        if (cruce[2][1]=i) or (cruce[1]=0 and cruce[2][4]=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)-1;
  for coloreo in colormejorvirtual(pd, biquandle,involutive) do
    W1 := One(Fnc);
    W2 := One(Fnc);
    for i in [1..N1] do
      for cruce in pd do
        if (cruce[2][1]=i) or (cruce[1]=0 and cruce[2][4]=i) then
          W1:=W1*BWgen(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    for j in [N1+1..max] do
      for cruce in pd do
        if (cruce[2][1]=j) or (cruce[1]=0 and cruce[2][4]=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 colormejorvirtual(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;
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 
return false;
fi;
od;od;od;
return true;
end;

s_biquandle:=function(q)
local x,y, n,s;
n:=Size(q.rperms);
s:=[];
for x in [1..n] do
Add(s,0);
od;

for x in [1..n] do
 for y in [1..n] do
  if r(q,x,y)=[x,y] then
   s[x]:=y;
  fi;
 od;
od;
if Size(Difference(s,[[0]]))<n then return false; fi;
return s;
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 sign 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;


##


##
listaisoCycleRel:=function(binv,biquandle)#da la lista de iso con identidad en biquandle.
local x,y,i,f,ff,n,lista,L,H,R,orb,Uq,Ui,iso;
n:=Size(biquandle.rperms);
lista:=[];

Uq:=qmorphismCycle(biquandle);
Ui:=qmorphismCycle(binv);

H:=IntersectionSet(Uq,Ui);

if Size(Uq)=1 then 
 iso:=Uq;
else
 if Size(H)=1 then
 iso:=AsList(Uq);
 else
 orb:=OrbitsDomain( AsGroup(IntersectionSet(Uq,Ui)), AsList(Uq), OnRight);
 iso:=List(orb,c->Representative(c));
 fi;
fi;


for f in iso do

 L := NullMat(n,n);
 R := NullMat(n,n);

 for x in [1..n] do
 for y in [1..n] do
 R[x][y]:=binv.rperms[x^f][y^f]^Inverse(f);
 L[x][y]:=binv.lperms[x^f][y^f]^Inverse(f);
 od;
 od;
 Add(lista, rec(rperms:=ShallowCopy(R),lperms:=ShallowCopy(L)) );
od;

return [biquandle,lista];
end;

repspair:=function(Q,listagrande)
#de una lista de reps. de pares Q, biquandles involutivos
local l1,l2,q,G;
#G:=qmorphismCycle(Q);
G:=qmorphism(Q);

l1:=[];
l2:=listagrande;

while l2<>[] do
 Print(" ",Size(l1),"-");
 q:=l2[1];
 Add(l1,q);
 l2:=Difference(l2,[q]);
 Print(Size(l2),"|");
 #l2:=Difference(l2,listaisoCycleRel(q,Q)[2]);
 l2:=Difference(l2,listaisoG(G,q));
od;
Print(Size(l1));
return l1;
end;


pairs:=function(lbiq,linv)
local lista,l, listaq,q,qi,Qi;

lista:=[];
for q in lbiq do
listaq:=[];
 for qi in linv do
# l:=listaisoCycleRel(qi,q )[2];
 l:=listaisoCycle(qi);
  for Qi in l do
  if checkPair(q,Qi) then 
  Add(listaq,Qi);
  fi;
  od;
 od;
#Add(lista, [q,listaq]);
Add(lista, [q,repspair(q,listaq)]);
od;
return lista;
end;

pairsfull:=function(lbiq,linv)
local lista,l, listaq,q,qi,Qi,c;

c:=Size(lbiq);
lista:=[];
for q in lbiq do
Print(c," ");
c:=c-1;
 for qi in linv do
 l:=listaisoCycle(qi);
  for Qi in l do
  if checkPair(q,Qi) then 
  Add(lista,[q,Qi]);
  fi;
  od;
 od;
od;

return lista;
end;



listaisoPairs:=function(p)
local x,y,i,f,ff,lista,Lb,Li,Rb,Ri,iso,n,u,U,Ub,Ui,orb;

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


#Ub:=qmorphismCycle(p[1]);
#Ui:=qmorphismCycle(p[2]);
#U:=IntersectionSet(Ub,Ui);

#if Size(U)=1 then 
 iso:=AsList(SymmetricGroup(n));
#else
# orb:=OrbitsDomain(AsGroup(U), AsList(SymmetricGroup(n)), OnRight);
# iso:=List(orb,c->Representative(c));
#fi;

for f in iso do
 Lb := NullMat(n,n);
 Rb := NullMat(n,n);
 Li := NullMat(n,n);
 Ri := NullMat(n,n);

 for x in [1..n] do
 for y in [1..n] do
 Rb[x][y]:=p[1].rperms[x^f][y^f]^Inverse(f);
 Lb[x][y]:=p[1].lperms[x^f][y^f]^Inverse(f);
 Ri[x][y]:=p[2].rperms[x^f][y^f]^Inverse(f);
 Li[x][y]:=p[2].lperms[x^f][y^f]^Inverse(f);
 od;
 od;
 Add(lista, [rec(rperms:=ShallowCopy(Rb),lperms:=ShallowCopy(Lb)),rec(rperms:=ShallowCopy(Ri),lperms:=ShallowCopy(Li))] );
od;

return lista;
end;

repsPairs:=function(listagrande)#de una lista de pares
local l1,l2,p;
l1:=[];
l2:=listagrande;

while l2<>[] do
Print(" ",Size(l1),"-");
p:=l2[1];
Add(l1,p);
l2:=Difference(l2,[p]);
Print(Size(l2),"|");

l2:=Difference(l2,listaisoPairs(p));
od;
Print(Size(l1));
return l1;
end;





#give the connected components of a biquandle
Tvirtual:=function(Q,bi)

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]] ]);
clases:=UnionSet(clases,[  [x,r(bi,x,y)[2]],[y,r(bi,x,y)[1]] ]);
od;
od;

clases:=relationgenerated(clases);

return clases;

end;


connectedpairs:=function(lbq,linv)
local A,q,qinv,qi,k;

A:=[];
for q in lbq do 
for qinv in linv do 
for qi in listaisoCycle(qinv) do 
k:=Size(Tvirtual(q,qi));Print(".");
 if k=1 and checkPair(q,qi) then Print("*");
  Add(A,[q,qi]);
 fi;
od;od;od;
return A;
end;
#


dominio:=function(q)
local i,j,n,SL,SR,D;
n:=Size(q.lperms);
SL:=[];
SR:=[];

for i in [1..n] do
if q.lperms[i]<>[] then Add(SL,i);fi;
od;
for i in [1..n] do
if q.rperms[i]<>[] then Add(SR,i);fi;
od;
D:=[];
for i in SL do
for j in SR do
Add(D,[i,j]);
od;od;
return D;
end;



#
listaiso:=function(biquandle)
local x,y,i,f,ff,lista,L,R,iso,n, Dom;

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

iso:=SymmetricGroup(n);

for f in iso do
 L := NullMat(n,n);
 R := NullMat(n,n);
 for x in [1..n] do
 for y in [1..n] do
 R[x][y]:=biquandle.rperms[x^f][y^f]^Inverse(f);
 L[x][y]:=biquandle.lperms[x^f][y^f]^Inverse(f);
 od;
 od;
 lista:=UnionSet(lista, ShallowCopy([rec(rperms:=ShallowCopy(R),lperms:=ShallowCopy(L))]) );
od;

return lista;
end;
#


#
listaisopartial:=function(uplas)
local x,y,f,lista,iso,R,N,n;

N:=Size(uplas);
n:=Size(uplas[1]);
lista:=[];

iso:=SymmetricGroup(N);

for f in iso do
 R := NullMat(N,n);
 for x in [1..N] do
 for y in [1..n] do
 R[x][y]:=uplas[x^f][y^f]^Inverse(f);
 od;
 od;
#Print(R);
 Add(lista, R);#ShallowCopy(R));
od;

return Set(lista);
end;
#







#de una lista de biquandles da una lista de representantes de clases de isomorfismo
repsCycle:=function(listagrande)
local l1,l2,q;
l1:=[];
l2:=listagrande;

while l2<>[] do
Print(" ",Size(l1),"-");
q:=l2[1];
Add(l1,q);
l2:=Difference(l2,[q]);
Print(Size(l2),"|");

l2:=Difference(l2,listaisoCycle(q));
od;
Print(Size(l1));
return l1;
end;
#

#
comparacion:=function(A,B)
local a,C,l;
C:=[];
for a in A do
l:=IntersectionSet(listaiso(a),B);
if Size(l)>1 then
    Add(C,l);
fi;
od;

return C;
end;




#
allpairs:=function(lbq,linv)
local A,B,q,qinv,qi,k,n,lq,li,liall;
n:=Size(lbq[1].rperms);

lq:=UnionSet(lbq,linv);
lq:=UnionSet(lbq,[bialexander(n,1,1)]);
lq:=repsCycle(lq);
li:=UnionSet(linv,[bialexander(n,1,1)]);

liall:=[];
for q in li do 
liall:=UnionSet(liall, listaisoCycle(q));
od; 
Print("\n","all involutives:",Size(liall),"\n");


A:=[];
for q in lq do 
B:=[];
for qi in liall do 
# k:=Size(Tvirtual(q,qi));Print(".");
# if k=1 and 
  if checkPair(q,qi) then Print("*");
  Add(B,qi);
 fi;
od;
Add(A,[q,isopar(q,B)]);
od;
return A;
end;





ret:=function(qi)
local A,x,y,n;
n:=Size(qi.rperms);

A:=[];
for x in [1..n] do
for y in [x+1..n] do
if qi.rperms[x]= qi.rperms[y] then Add(A,[x,y]);
fi;
od;
od;
for x in [1..n] do
Add(A,[x,x]);
od;

return relationgenerated(A);
end;
#

ListPermn:=function(p,n)
local l,i,s;
 l:=ListPerm(p);
 s:=Size(l);
 for i in [s+1..n] do
  l[i]:=i;
 od;
 return l;
end;


#q is involutive, f\in Aut : f^2=1, r(x,y):=r(fx,fy)
iS:=function(f)
local x,y,rperms,q,lperms,n;
  n:=Size(f);
  q:=bialexander(n,1,1);
  lperms := NullMat(n,n);
  rperms := NullMat(n,n);
 for x in [1..n] do
 for y in [1..n] do
 lperms[x][y]:=f[q.lperms[x][y]];
 rperms[x][y]:=ListPermn(Inverse(PermList(f)),n)[q.lperms[x][y]];
# rperms[x][y]:=f[q.lperms[x][y]];
 od;od;
 return rec( lperms := lperms, rperms := rperms);
end;


#
virtualis:=function(biquandle)
local q,U,f,n,A,M,CC;
q:=biquandle;
M:=qmorphismCycle(biquandle);
n:=Size(q.rperms);

U:=[];
for f in M do
Add(U,PermList(f));
od;

if Size(U)=1 then 
return [[q,iS(M[1])]];
else
U:=AsGroup(U);
A:=[];
CC:=List(ConjugacyClasses(U),c->Representative(c));

for f in CC do
Add(A,[q,iS(ListPermn(f,n))]);
od;
return A;
fi;
end;

listvirtualis:=function(listbiq)
local A,q;
A:=[];
for q in listbiq do
A:=UnionSet(A,virtualis(q));
od;
return A;
end;



#devuelve true si todas las igualdades de la compatibilidad son ciertas, devuelve false si no.
checkSingular:=function(bi,tau)
local x,y,z,v,w,r12,r23,rt12,rt23,N;
N:=Size(bi.rperms);

# chequeo de S tau=tau S
for x in [1..N] do;
for y in [1..N] do;
 v:=r(bi,x,y);
 w:=r(tau,x,y);
 if
  r(bi,w[1],w[2])<> r(tau,v[1],v[2]) then 
  return false;
 fi;
od;
od;

#hay que chequear (1xtau)(Sx1)(1xS)=(Sx1)(1xS)(1xtau)
#                (tau x1)(1xS)(Sx1)=(1xS)(Sx1)(tau x1)
for x in [1..N] do;
for y in [1..N] do;
for z in [1..N] do;
r12:=r(bi,x,y);
r23:=r(bi,y,z);
rt12:=r(tau,x,y);
rt23:=r(tau,y,z);

if (
[r(tau,r12[1],r(bi,r12[2],z)[1])[1],
r(tau,r12[1],r(bi,r12[2],z)[1])[2],
r(bi,r12[2],z)[2]]
<>
[r(bi,x,rt23[1])[1],
r(bi,r(bi,x,rt23[1])[2],rt23[2])[1],
r(bi,r(bi,x,rt23[1])[2],rt23[2])[2]
 ])
or
##las 1eras tres son RIV(a) las 2das 3 son RIV(b)
([r(bi,x,r23[1])[1],
r(tau,r(bi,x,r23[1])[2],r23[2])[1],
r(tau,r(bi,x,r23[1])[2],r23[2])[2]
]
<>
[r(bi,rt12[1],r(bi,rt12[2],z)[1])[1],
r(bi,rt12[1],r(bi,rt12[2],z)[1])[2],
r(bi,rt12[2],z)[2] ])
then 
return false;
fi;
od;od;od;
return true;
end;



#de una lista de biquandles da una lista de representantes de clases de isomorfismo
repspartial:=function(listagrande)
local l1,l2,q;
l1:=[];
l2:=listagrande;

while l2<>[] do
Print(" ",Size(l1),"-",Size(l2)," ");
q:=l2[1];
Add(l1,q);
l2:=Difference(l2,[q]);
l2:=Difference(l2,listaisopartial(q));
od;

return l1;
end;


checkBijective:=function(tau)
local n,x,y,A;
n:=Size(tau.rperms);
A:=[];
for x in [1..n] do
for y in [1..n] do
Add(A, r(tau,x,y));
od;od;

return Size(Set(A))=n^2;
end;




singularFlipList := function(n)
 local Sn, sigma,set, singlist, l1,lperms,tau ,cclasses,CC,all;

Sn:=SymmetricGroup(n);
set:=[]; 
for sigma in Sn do
Add(set, ListPermn(sigma,n));
od;

cclasses:=List(ConjugacyClasses(Sn),c->Representative(c));
CC:=[];
for sigma in cclasses do
Add(CC, ListPermn(sigma,n));
od;


singlist:=[];

CC:=set;

all:=[];
for l1 in CC do
for lperms in Tuples(set,n-1) do;
Add(lperms,l1);
Add(all, lperms);
od;
od;
#Print(all);

for lperms in all do
tau:=rec(lperms:=lperms, rperms:=lperms ) ;
if checkBijective(tau) then Add(singlist, tau); fi;
od;

return singlist;
end;

checkpb:=function(upla,n)
local i,j,A;
A:=[];
for i in [1..n] do
for j in [1..n] do
Add(A,[upla[i][j],upla[j][i]]);
od;
od;
return Size(Set(A))=n^2;
end;

checkpb2:=function(pair)
local s,t,A;
s:=pair[1];
t:=pair[2];
A:=[[s[1],s[1]],[s[2],t[1]],[t[1],s[2]],[t[2],t[2]]];
return Size(Set(A))=4;
end;

checkpb3:=function(terna)
local s,t,A;
s:=terna[1];
t:=terna[2];
t:=terna[3];
A:=[[s[1],s[1]],[s[2],t[1]],[t[1],s[2]],[t[2],t[2]]];
#hacer el chequeo biy en tres!
return Size(Set(A))=4;
end;




FlipList := function(n)
local Sn, sigma,set, singlist, l1,l2,lperms,tau ,cclasses,CC,all,all3,all4,all5,pair,terna,quator;

Sn:=SymmetricGroup(n);
set:=[]; 
for sigma in Sn do
Add(set, ListPermn(sigma,n));
od;

cclasses:=List(ConjugacyClasses(Sn),c->Representative(c));
CC:=[];
for sigma in cclasses do
Add(CC, ListPermn(sigma,n));
od;

#set: todas las permutaciones, CC=reps de clases de conj

singlist:=[];

all:=[];
for l1 in set do
for l2 in set do
if checkpb2([l1,l2]) then Add(all, [l1,l2]);Print(".");
else Print("*"); fi;
od;
od;

Print("\n size all2 =",Size(all),"\n");
all:=repspartial(all);
Print("\n size all2 =",Size(all),"\n");

if n=2 then return all;fi;

if n>2 then
all3:=[];
for pair in all do
for tau in set do
if
checkpb([pair[1],pair[2],tau],3) then
Add(all3,[pair[1],pair[2],tau]);Print(".");
else Print("*");
fi;
od;
od;

Print("\n size all3 =",Size(all3),"\n");
all3:=repspartial(all3);
Print("\n size all3 =",Size(all3),"\n");


if n=3 then return all3; 
else

all4:=[];
for terna in all3 do
for tau in set do
if
checkpb([terna[1],terna[2],terna[3],tau],4) then
Add(all4,[terna[1],terna[2],terna[3],tau]);#Print(".");
#else Print("*");
fi;
od;
od;

Print("\n size all4 =",Size(all4),"\n");
all4:=repspartial(all4);
Print("\n size all4 =",Size(all4),"\n");


if n=4 then return all4; 
else

all5:=[];
for quator in all4 do
for tau in set do
if
checkpb([quator[1],quator[2],quator[3],quator[4],tau],5) then
Add(all5,[quator[1],quator[2],quator[3],quator[4],tau]);
else Print("*");Print(".");
fi;
od;
od;

#Print("\n size all5 =",Size(all5),"\n");
#all5:=repspartial(all5);
#Print("\n size all5 =",Size(all5),"\n");

#HACER LISTA CLASE DE ISO HASTA UN PEQUEÑO n=2,3,4.. 
#

#for lperms in all do
#tau:=rec(lperms:=lperms, rperms:=lperms ) ;
#if checkBijective(tau) then Add(singlist, tau); fi;
#od;

return all5;

fi;
fi;
fi;
end;














imagen:=function(q)
local n,x,y,A;
 n:=Size(q.lperms);
A:=[];
for x in [1..n] do 
for y in [1..n] do 
Add(A,r(q,x,y));
od;od;
return(Set(A));
end;







#de una lista de biquandles da una lista de representantes de clases de isomorfismo
reps:=function(listagrande)
local l1,l2,q;
l1:=[];
l2:=listagrande;

while l2<>[] do
Print(" ",Size(l1),"-",Size(l2),"--");
q:=l2[1];
Add(l1,q);
l2:=Difference(l2,[q]);
l2:=Difference(l2,listaiso(q));
od;

return l1;
end;




checkPartialInj:=function(q)
local D,u,i,j,I;
D:=dominio(q);
I:=[];
for u in D do
i:=u[1];
j:=u[2];
Add(I,r(q,i,j));
od;
if Size(Set(I))=Size(D) then return true;
else return false;
fi;
end;


generateBijective:=function(n)
local q,q0,perms0,i,j, initialList, Sn,sigma,set,s1,s2,s3,s4;

perms0:=NullMat(n,n);
for i in [1..n] do
#for j in [1..n] do
#perms0[i][j]:=[];
perms0[i]:=[];
od;#od;

Sn:=SymmetricGroup(n);
set:=[]; 
for sigma in Sn do
Add(set, ListPermn(sigma,n));
od;

initialList:=[];
q0:=rec(lperms:=StructuralCopy(perms0),rperms:=StructuralCopy(perms0));

for s1 in set do
for s2 in set do
for s3 in set do
#for s4 in set do
q:=StructuralCopy(q0);
q.lperms[1]:=ShallowCopy(s1);
q.lperms[2]:=ShallowCopy(s2);
q.rperms[1]:=ShallowCopy(s3);
#q.rperms[2]:=ShallowCopy(s4);
if checkPartialInj(q)
then 
Add(initialList,q);
fi;
q:=StructuralCopy(q0);
od;od;od;#od;

return initialList;
end;


generateBijectiveFlip:=function(n)
local q,q0,perms0,i,j, initialList, Sn,sigma,set,s1,s2,s3,s4;

perms0:=NullMat(n,n);
for i in [1..n] do
#for j in [1..n] do
#perms0[i][j]:=[];
perms0[i]:=[];
od;#od;

Sn:=SymmetricGroup(n);
set:=[]; 
for sigma in Sn do
Add(set, ListPermn(sigma,n));
od;

initialList:=[];
q0:=rec(lperms:=StructuralCopy(perms0),rperms:=StructuralCopy(perms0));

for s1 in set do
for s2 in set do
q:=StructuralCopy(q0);
q.lperms[1]:=ShallowCopy(s1);
q.lperms[2]:=ShallowCopy(s2);
q.rperms[1]:=ShallowCopy(s1);
q.rperms[2]:=ShallowCopy(s2);
if checkPartialInj(q)
then 
Add(initialList,q);
fi;
q:=StructuralCopy(q0);
od;od;

return initialList;
end;




#devuelve true si todas las igualdades de la compatibilidad son ciertas, devuelve false si no.
checkSingularParcial:=function(bi,tau)
local u,x,y,z,v,w,r12,r23,rt12,rt23,N, Dom;
N:=Size(bi.rperms);
Dom:=dominio(tau);


# chequeo de S tau=tau S
for u in Dom do
x:=u[1];
y :=u[2];
 v:=r(bi,x,y);
 w:=r(tau,x,y);
 if
  ((r(bi,w[1],w[2]) in Dom) and (v in Dom)) then
    if r(bi,w[1],w[2]) <> r(tau,v[1],v[2]) then    return false; fi;
 fi;
od;

#hay que chequear (1xtau)(Sx1)(1xS)=(Sx1)(1xS)(1xtau)
#                (tau x1)(1xS)(Sx1)=(1xS)(Sx1)(tau x1)
for u in Dom do
x:=u[1];
y:=u[2];

for z in [1..N] do;
if [y,z] in Dom then
r12:=r(bi,x,y);
r23:=r(bi,y,z);

 rt12:=r(tau,x,y);
 rt23:=r(tau,y,z);

 w:=[r12[1],r(bi,r12[2],z)[1]];
 if w in Dom then
  if (
   [r(tau,w[1],w[2])[1],
   r(tau,w[1],w[2])[2],
   r(bi,r12[2],z)[2]]
   <>
   [r(bi,x,rt23[1])[1],
   r(bi,r(bi,x,rt23[1])[2],rt23[2])[1],
   r(bi,r(bi,x,rt23[1])[2],rt23[2])[2]
    ])
   then 
   return false;
  fi;
 fi;


 w:=[r(bi,x,r23[1])[2],r23[2]];
 if w in Dom then
   ##las 1eras tres son RIV(a) las 2das 3 son RIV(b)
   if
   ([r(bi,x,r23[1])[1],
   r(tau,w[1],w[2])[1],
   r(tau,w[1],w[2])[2]
   ]
   <>
   [r(bi,rt12[1],r(bi,rt12[2],z)[1])[1],
   r(bi,rt12[1],r(bi,rt12[2],z)[1])[2],
   r(bi,rt12[2],z)[2] ])
   then 
   return false;
  fi;
 
 fi;

fi;
od;
od;
return true;
end;


#
generateSingularPair:=function(biq)
local A,B,tau,n,q,q0,q1,q2,q3,perms0,i,j, initialList, Sn,sigma,set,s1,s2,s3,s4;

n:=Size(biq.lperms);
if n=2 then Print("OJO n=2!");else

initialList:=[];
for tau in generateBijective(n) do
if checkSingularParcial(biq,tau) then
Add(initialList,tau);
fi;
od;
Print("Size(initial List)= ", Size(initialList),"\n");

Sn:=SymmetricGroup(n);
set:=[]; 
for sigma in Sn do
Add(set, ListPermn(sigma,n));
od;

A:=[];
for q0 in initialList do
 for s1 in set do
 q1:=StructuralCopy(q0);
 q1.rperms[2]:=ShallowCopy(s1);
 if checkPartialInj(q1) and checkSingularParcial(biq,q1)  then 
 for s2 in set do
 q2:=StructuralCopy(q1);
 q2.lperms[3]:=ShallowCopy(s2);
 if checkPartialInj(q2) and checkSingularParcial(biq,q2)  then 
 for s3 in set do
# q:=StructuralCopy(q0);
# q.rperms[2]:=ShallowCopy(s1);
# q.lperms[3]:=ShallowCopy(s2);
 q3:=StructuralCopy(q2);
 q3.rperms[3]:=ShallowCopy(s3);
 if checkPartialInj(q3) and checkSingularParcial(biq,q3)  then 
 Add(A,q3);
 fi;
# q:=StructuralCopy(q0);
 od;
 fi;od;
 fi;od;
od;
Print("Size(second List)= ", Size(A),"\n");

if n=3 then return A;
else

B:=[];
for q0 in A do
 for s1 in set do
 for s2 in set do
 q:=StructuralCopy(q0);
 q.lperms[4]:=ShallowCopy(s1);
 q.rperms[4]:=ShallowCopy(s2);
 if checkPartialInj(q) and checkSingularParcial(biq,q) 
 then 
 Add(B,q);
 fi;
 q:=StructuralCopy(q0);
 od;od;
od;

Print("Size(third List)= ", Size(B),"\n");
return B;
fi;
fi;end;


#
upla2rec:=function(upla)
return rec(lperms:=upla,rperms:=upla);
end;

#
uplas2list:=function(uplas)
local upla,A;
A:=[];
for upla in uplas do
Add(A,upla2rec(upla));
od;
return A;
end;




#gives the permutation (as a list) of the multiplication by lambda in Z/nZ
multlambda:=function(n,lambda)
local e,i,ms,x;
e := Enumerator(ZmodnZ(n));
ms:=ListPermn((),n);

for x in e do
ms[Position(e,x)]:=Position(e,lambda*x);
od;
return ms;
end;


tauphi:=function(phi,n,s,t)
local x,y,u,v,e,lperms,rperms;

e := Enumerator(ZmodnZ(n));
lperms := NullMat(n,n);
rperms := NullMat(n,n);

 for x in e do
    for y in e do
u:=s* e[phi[Position(e,y-x/s)]]+x;
v:=t*e[phi[Position(e,x-s*y)]]+y;
      lperms[Position(e, x)][Position(e, y)] := Position(e,u);
      rperms[Position(e, y)][Position(e, x)] := Position(e,v);
    od;
  od;
return rec( lperms := lperms, rperms := rperms);
end;


#los tau singulares compatibles con bialexander cuando (1-st) es una unidad
singularbialexander := function(n, s, t)
local ms,mt,m1,Sn,G,listphi,g,phi,e,x,y,u,v,listsing,lperms,rperms,tau;
  
Print("\n n=",n,", 1-st=",1-s*t,"\n");

ms:=PermList(multlambda(n,s));
mt:=PermList(multlambda(n,t));
m1:=PermList(multlambda(n,-1));
Print("ms=",ms,", mt=",mt,", m(-1)=",m1,"\n");

Sn:=SymmetricGroup(n);
G:=Centralizer(Sn,m1);
G:=Centralizer(G,ms);
G:=Centralizer(G,mt);
Print("Size(Centralizer)=",Size(G),"\n");


listphi:=[];
for g in G do
Add(listphi, ListPermn(g,n));
od;

listsing:=[];
for phi in listphi do
Print("*");
tau:=tauphi(phi,n,s,t);
if checkBijective(tau) then
Add(listsing, tauphi(phi,n,s,t));Print(".");
fi;
od;


return listsing;
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
BWgenSV := function(cruce, coloreo, clases,Fnc)
local bwg,i;

bwg:=One(Fnc);
  if cruce[1]=1 then 
    for i in [1..Size(clases)] do
      if ["f",[coloreo[cruce[2][1]],coloreo[cruce[2][4]]]] in clases[i] then
        bwg := GeneratorsOfGroup(Fnc)[i];
      fi;
    od;
  fi;
# de alguna forma hay que poner ["f",[cruce coloreado]] in classes[i] (caso 1), o ["h",... (caso 2) o ["g",... (caso 0)

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


  if cruce[1]=2 then 
    for i in [1..Size(clases)] do
#Print("\n\n");
#Print(clases[i],"\n");
#Print(["h",[coloreo[cruce[2][1]],coloreo[cruce[2][4]]]],"\n");
      if ["h",[coloreo[cruce[2][1]],coloreo[cruce[2][4]]]] in clases[i] then
        bwg := GeneratorsOfGroup(Fnc)[i];
      fi;
    od;
  fi;

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

return bwg;
end;




invariantSingular:=function(pd,biquandle,tau)
local generators,
inv,max, W, W1,W2, W3, BWg, coloreo,  Un, cruce, Fnc,
n, rk, N1, N2,i, j,k,l;

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

#calculates U_nc equations and generators
Un:=uncfh([biquandle,tau]);
generators:=Un.generators;
generators:=Difference(generators,[generators[1]]); #le quitamos a [] como generador
rk:=Size(generators);
Fnc := FreeGroup(rk);

#and given a coloring calculates the invariant

n:=componentesconexas(pd);
#n:=1;
 if n>3 then Display("ATENTION: More than 3 connected components");fi;

if n=1 then
  for coloreo in colormejorsingular(pd, biquandle,tau) do
    W := One(Fnc);
    for i in [1..max] do
      for cruce in pd do
        if (cruce[2][1]=i) or (cruce[1]=2 and cruce[2][4]=i) then
          W:=W*BWgenSV(cruce, coloreo, generators,Fnc);
        fi;
      od;
    od;
    Add(inv,[coloreo,W]);
  od;
  return(inv);
elif n=2 then
  N1:=n1(pd)-1;
  for coloreo in colormejorsingular(pd, biquandle,tau) do
    W1 := One(Fnc);
    for i in [1..N1] do
      for cruce in pd do
        if (cruce[2][1]=i) or (cruce[1]=2 and cruce[2][4]=i) then
          W1:=W1*BWgenSV(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    W2 := One(Fnc);
    for j in [N1+1..max] do
      for cruce in pd do
        if (cruce[2][1]=j) or (cruce[1]=2 and cruce[2][4]=j) then
          W2:=W2*BWgenSV(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 colormejorsingular(pd, biquandle,tau) 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*BWgenSV(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*BWgenSV(cruce, coloreo, generators, Fnc);
        fi;
      od;
    od;
    Add(inv,[coloreo,[W1,W2,W3]]);
  od;
return inv;
fi;

end;


ordenar:=function(list)
local S,A,a,i,c;
A:=[];
S:=Set(list);

for a in S do
 c:=0;
 for i in [1..Size(list)] do
  if a=list[i] then c:=c+1; fi;
 od;
 for i in [1..c] do
  Add(A,a);
 od;
od;
return A;
end;

#quita un elemento x de una lista. (si x esta n veces, entonces queda n-1 veces)
#y la devuelve ordenada
quitar:=function(list,x)
local S,a,i,c,A;

S:=Set(list);
A:=[];

for a in S do
 c:=0;
 for i in [1..Size(list)] do
  if a=list[i] then c:=c+1; fi;
 od;
if a=x then
 for i in [1..c-1] do
  Add(A,a);
 od; 
else
 for i in [1..c] do
  Add(A,a);
 od;
fi;

od;

return A;
end;


recycle_eq_ab:=function(S,Eq)
local eq,a,aa,SS,c,reps,u,EE,EEE,i,b,d,x,p1,p2,l,Int;
SS:=S;

for eq in Eq do
if Size(eq[1])=1 then
Add(SS, [eq[1][1],eq[2][1]]);
fi;
od;

SS:=relationgenerated(SS);

EE:=[];
for eq in Eq do
l:=Size(eq[1]);
for i in [1..l] do
eq[1][i]:=ShallowCopy(rep(SS,eq[1][i]));
eq[2][i]:=ShallowCopy(rep(SS,eq[2][i]));
od;
Add(EE,eq);
od;



for eq in EE do
Int:=ShallowCopy(IntersectionSet(eq[1],eq[2]));
for x in Int do
eq[1]:=quitar(eq[1],x);
eq[2]:=quitar(eq[2],x);
od;
od;

EE:=Set(EE);

for eq in EE do
eq[1]:=ordenar(eq[1]);
eq[2]:=ordenar(eq[2]);
od;

EE:=Set(EE);

EEE:=[];
for eq in EE do
Add(EEE,ordenar(eq));
od;

EE:=Difference(EEE,[[[],[]]]);
#SS:=relationgenerated(UnionSet(SS,LookWithInverses(EE)));


return  rec(generators:=SS,equations:=EE);
end;





gen_eq_ab:=function(q,tau)
local n,x,y,z,S,E,eq,s,t, reps,u,e1,e2;
n:=Size(q.rperms);

s:=s_biquandle(q);


E:=[];
#f2')
for x in [1..n] do
Add(E,[[["f",[x,s[x]]] ],[[]]]);
for y in [1..n] do
#c3')
u:=[
[  ["f" ,[x,y]] , ["h",r(q,x,y)] ]
,
[  ["f",r(tau,x,y)],["h",[x,y]] ]
];
Add(E,ShallowCopy(u)); 
for z in [1..n] do
#[f1')
e1:=
[  ["f" ,[x,y]] , ["f",[r(q,x,y)[2],z]] ,["f",[r(q,x,y)[1] ,r(q,r(q,x,y)[2],z)[1] ]] ];
e2:=
[  ["f",[x,r(q,y,z)[1]] ],["f", [r(q,x,r(q,y,z)[1] )[2],r(q,y,z)[2] ] ], ["f",[y,z]] ];
u:=[ordenar(e1),ordenar(e2)];
#!
#Print(Size(e1)=Size(e2));
Add(E,ShallowCopy(u)); 
#c1')
e1:=[  ["h" ,[y,z]] , ["f",[x,r(tau,y,z)[1] ]] ,["f",[r(q,x,r(tau,y,z)[1])[2],r(tau,y,z)[2] ]] ];
e2:=[ ["f",[x,y]], ["f",[r(q,x,y)[2],z] ],["h", [r(q,x,y)[1],r(q,r(q,x,y)[2],z)[1] ] ]  ];
u:=[ordenar(e1),ordenar(e2)];
#Print(Size(e1)=Size(e2));
Add(E,ShallowCopy(u)); 

#c2')
e1:=[ ["f",[y,z]], ["f",[x,r(q,y,z)[1]] ],["h", [r(q,x,r(q,y,z)[1])[2],r(q,y,z)[2] ] ]  ];
e2:=[ ["h" ,[x,y]] , ["f",[r(tau,x,y)[2],z ]] ,["f",[r(tau,x,y)[1],r(q,r(tau,x,y)[2],z)[1] ]] ];
u:=[ordenar(e1),ordenar(e2)];
#!
#Print(Size(e1)=Size(e2));
Add(E,ShallowCopy(u)); 

#conmutativity
#for t in [1..n] do
#Add(E, [ [ ["f",[x,y]],["f",[z,t]] ] , [["f",[z,t]],["f",[x,y]]] ]);
#Add(E, [ [ ["h",[x,y]],["f",[z,t]] ] , [["f",[z,t]],["h",[x,y]]] ]);
#Add(E, [ [ ["h",[x,y]],["h",[z,t]] ] , [["h",[z,t]],["h",[x,y]]] ]);
#od;

od;od;od;


E:=Set(E);

for eq in E do
if Size(Set(eq))=1 then E:=Difference(E,[eq]);fi;
if Size(Set(eq))=2 then 
E:=Difference(E,[eq]);
Add(E,ShallowCopy(Set(eq)));
fi;
od;


return Set(E);
end;
#


ge:=function(n)
local S,x,y;

S:=[[[]]];
for x in [1..n] do
for y in [1..n] do
Add(S,[["f",[x,y]]]);
Add(S,[["h",[x,y]]]);
od;od;
return S;
end;

Abfh:=function(q,tau)#TERMINAR LOS GEN_EQ_F/FSTAU/H!!!!!!!!!!!!!!!!!!!!!
local k,i,Eq,S,Ab,Ab2,x,y,n;
n:=Size(q.rperms);


S:=[[[]]];
for x in [1..n] do
for y in [1..n] do
Add(S,[["f",[x,y]]]);
Add(S,[["h",[x,y]]]);
od;od;

k:=Size(S);

Eq:=gen_eq_ab(q,tau);

Ab:=rec(generators:=S,equations:=Eq);

for i in [1..k] do
#Print(i,"\n Un=",Un,"\n\n\n");
Ab2:=recycle_eq_ab(S,Eq);
if Ab2=Ab then return Ab;fi;
Ab:=Ab2;
S:=Ab.generators;
Eq:=Ab.equations;
od;
return Ab;
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
BWgenSV := function(cruce, coloreo, clases,Fnc)
local bwg,i;

bwg:=One(Fnc);
  if cruce[1]=1 then 
    for i in [1..Size(clases)] do
      if ["f",[coloreo[cruce[2][1]],coloreo[cruce[2][4]]]] in clases[i] then
        bwg := GeneratorsOfGroup(Fnc)[i];
      fi;
    od;
  fi;
# de alguna forma hay que poner ["f",[cruce coloreado]] in classes[i] (caso 1), o ["h",... (caso 2) o ["g",... (caso 0)

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


  if cruce[1]=2 then 
    for i in [1..Size(clases)] do
#Print("\n\n");
#Print(clases[i],"\n");
#Print(["h",[coloreo[cruce[2][1]],coloreo[cruce[2][4]]]],"\n");
      if ["h",[coloreo[cruce[2][1]],coloreo[cruce[2][4]]]] in clases[i] then
        bwg := GeneratorsOfGroup(Fnc)[i];
      fi;
    od;
  fi;

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

return bwg;
end;



BWcruces:= function(link, coloreo, clases,Fnc)
local bwg,cruce;

bwg:=One(Fnc);
for cruce in link do
bwg := bwg*BWgenSV(cruce, coloreo, clases,Fnc);
od;
return bwg;
end;

#hay que definir Abfh:=function(q,tau) que deberia devolver un grupo abeliano libre,
# donde se habran hecho alguinas simplificaciones, y una lista de relaciones

statesum:=function(link, q,tau)
local Ab,generators,rk,Fnc,col,clases,state,state0,state1,c,bw0,bw1;

#calculates Ab^fh equations and generators
Ab:=Abfh(q,tau);
generators:=Ab.generators;
generators:=Difference(generators,[generators[1]]); #le quitamos a [] como generador
rk:=Size(generators);
clases:=generators;

Fnc := FreeAbelianGroup(rk); #abelianizado..?

state0:=[];
for col in colormejorsingular(link,q,tau) do
Add(state0,BWcruces(link, col, clases,Fnc));
od;
state1:=Set(state0);


state:=[];
for bw0 in state1 do
c:=0;
for bw1 in state0 do
if bw0=bw1 then c:=c+1; fi;
od;
Add(state,[c,bw0]);
od;
return state;
end;