read("biquandles-maple.txt");

getrow:=proc(i,M)
#
# extracts row i from M
#
out:=[];
for j from 1 to coldim(M) do
 out:=[op(out),eval(M[i,j])];
od;
eval(out);
end;

getcol:=proc(i,M)
#
# extracts column i from M
#
out:=[];
for j from 1 to rowdim(M) do
 out:=[op(out),eval(M[j,i])];
od;
eval(out);
end;


zerow:=proc(v)
#
# check whether v is a zero vector
#
j:=true;
for i from 1 to coldim(matrix([v])) do
 if v[i] <> 0 then
  j:=false;
 fi;
od;
eval(j);
end;


getkernel:=proc(N)
#
# read off a basis for the kernel of N
#
M:=gaussjord(N);
out:=[];
lp:=[];
lead:={}; # set of column numbers of leading ones
for i from 1 to rowdim(matrix(M)) do
 j:=1;
 while M[i,j]=0 and j < coldim(M) do
  j:=j+1;
 od;
 if M[i,j]=1 then
  lead:=lead union {eval(j)};
  lp:=[op(lp),eval(j)];
 fi;
od;
for i from 1 to coldim(M) do
 if not eval(i) in lead then
  v:=[];
  for l from 1 to coldim(M) do
   v:=[op(v),0];
  od;
 for z from 1 to coldim(matrix([lp])) do
  v[lp[z]]:=eval(M[z,i]);
  od;
  v[i]:=-1;
  out:=[op(out),eval(v)];
 fi;
od;
eval(out);
end;

inspan:=proc(v,L)
#
# Checks whether v is in Span(L)
#
rowspace(matrix(L),'d1');
rowspace(matrix([op(L),v]),'d2');
if d1<>d2 then
  eval(false);
 else
  eval(true);
fi;
end;


ybcocy := proc(B)
#
# Computes the YB 2-cocycles of biquandle B
#
j:=1;
z1:=[];
for i from 1 to coldim(evalm(B[1])) do
 for j from 1 to coldim(evalm(B[1])) do
  z1:=[op(z1),0];
 od;
od;
n:=coldim(evalm(B[1]));
E:=[];
for i from 1 to n do
 for j from 1 to n do
  for k from 1 to n do
  p:=eval(z1);
  p[n*(i-1)+j]:= p[n*(i-1)+j]+1;
  p[n*(eval(B[2][i,j])-1)+k]:=p[n*(eval(B[2][i,j])-1)+k]+1;
  p[n*(eval(B[4][j,i])-1)+eval(B[4][k,B[2][i,j]])]:=
   p[n*(eval(B[4][j,i])-1)+eval(B[4][k,B[2][i,j]])]+1;
  p[n*(j-1)+k]:=p[n*(j-1)+k]-1;
  p[n*(i-1)+eval(B[4][k,j])]:=p[n*(i-1)+eval(B[4][k,j])]-1;
  p[n*(eval(B[2][i,B[4][k,j]])-1)+eval(B[2][j,k])]:=
   p[n*(eval(B[2][i,B[4][k,j]])-1)+eval(B[2][j,k])]-1;
  E:=[op(E),eval(p)];
  od;
 od;
od;
eval(getkernel(E));
end;

ybcob:=proc(B)
#
# Find YB 2-coboundaries
#
out:=[];
for i from 1 to coldim(matrix(B[1])) do
 cob:=[];
 for j from 1 to coldim(matrix(B[1])) do
  for k from 1 to coldim(matrix(B[1])) do
   temp:=0;
   if k=i then 
    temp:=temp+1;
   fi;
   if j=i then
    temp:=temp+1;
   fi;
   if B[2][j,k]=i then
    temp:=temp-1;
   fi;
   if B[4][k,j]=i then
    temp:=temp-1;
   fi;
   cob:=[op(cob),eval(temp)];
  od;
 od;
 out:=[op(out),eval(cob)];
od;
M:=gaussjord(out);
out2:=[];
for i from 1 to rowdim(M) do
 j:=getrow(i,M);
 if not zerow(j) then
  out2:=[op(out2),eval(j)];
 fi;
od;
eval(out2);
end;


ybcohom:=proc(B)
#
# Computes the Yang-Baxter cohomology with Q coefficients
#
Z:=ybcocy(B);
Bn:=ybcob(B);
out:=[];
while coldim(matrix([Z]))<>0 do
 x:=eval(Z[1]);
 Z:=subsop(1=NULL,Z);
 if not inspan(eval(x),Bn) then
  out:=[op(out),eval(x)];
  Bn:=[op(Bn),eval(x)]
 fi;
od;
eval(out);
end;


typeI:=proc(v,B)
#
# Tests yb cocycle for type I condition
#
t1:=true;
for a from 1 to coldim(B[1]) do
 x:=0;
 y:=0;
 for j from 1 to coldim(B[1]) do
  if B[2][j,a]=a then 
   x:=eval(j);
  fi;
  if B[4][j,a]=a then 
   y:=eval(j);
  fi;
 od;
 if v[coldim(B[1])*(x-1)+a]<>0 then
  t1:=false;
 fi;
 if v[coldim(B[1])*(a-1)+y]<>0 then
  t1:=false;
 fi;
od;
eval(t1);
end;


redybcohom:=proc(B)
#
# list of YB cocycles which satisfy type I condition
#
L:=ybcohom(B);
out:=[];
for X in L do
 if typeI(X,B) then
  out:=[op(out),eval(X)];
 fi;
od;
eval(out);
end;


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


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


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

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


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


ybinv:=proc(v,phi,B)
#
#  Compute Yang-Baxter cocyle invariant for Gauss code v
#  where phi is a type I cocycle in C^2(B)
#
v2:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  v2:=[op(v2),eval(v[i])];
 fi;
od;
b:=[]; # position in v2 of corresponding crossing label
for i from 1 to coldim(matrix([v2])) do
 j:=1;
 while v2[i]+v2[j]<>0 do
  j:=j+1;
 od;
 b:=[op(b),eval(j)];
od;
nx:=[];
lt:=1;
c:=0;
for i from 1 to coldim(matrix([v]))-1 do
 if v[i+1]=0 then
  nx:=[op(nx),eval(lt)];
  lt:=eval(i+1-c);
  c:=c+1;
 elif v[i]<>0 then
  nx:=[op(nx),eval(i+1-c)];
 fi;
od;
#C:=bhomlist(gauss2biq(v),B);
C:=bhomlist2(v,B);
n:=coldim(B[1]);
f:=0;
for c in C do
 z:=0;
 for i from 1 to coldim(matrix([v2])) do
  if Re(v2[i])<0 then
   if Im(v2[i]) =0 then 
    z:=z+phi[n*(eval(c[i])-1)+eval(c[b[i]])];
   elif Im(v2[i])<>0 then
    z:=z-phi[n*(eval(c[nx[i]])-1)+eval(c[nx[b[i]]])];
   fi;
  fi;
 od;
 f:=f+t^eval(z);
od;
eval(f);
end;


ybinv2:=proc(v,B)
#
#  Compute Yang-Baxter cocyle invariant for Gauss code v
#  for all type I cocycles in C^2(B)
#
P:=redybcohom(B);
v2:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  v2:=[op(v2),eval(v[i])];
 fi;
od;
b:=[]; # position in v2 of corresponding crossing label
for i from 1 to coldim(matrix([v2])) do
 j:=1;
 while v2[i]+v2[j]<>0 do
  j:=j+1;
 od;
 b:=[op(b),eval(j)];
od;
nx:=[];
lt:=1;
c:=0;
for i from 1 to coldim(matrix([v]))-1 do
 if v[i+1]=0 then
  nx:=[op(nx),eval(lt)];
  lt:=eval(i+1-c);
  c:=c+1;
 elif v[i]<>0 then
  nx:=[op(nx),eval(i+1-c)];
 fi;
od;
#C:=bhomlist(gauss2biq(v),B);
C:=bhomlist2(v,B);
n:=coldim(B[1]);
out:=[];
for phi in P do
 f:=0;
 for c in C do
  z:=0;
  for i from 1 to coldim(matrix([v2])) do
   if Re(v2[i])<0 then
    if Im(v2[i]) =0 then 
     z:=z+phi[n*(eval(c[i])-1)+eval(c[b[i]])];
    elif Im(v2[i])<>0 then
     z:=z-phi[n*(eval(c[nx[i]])-1)+eval(c[nx[b[i]]])];
    fi;
   fi;
  od;
  f:=f+t^eval(z);
 od;
 out:=[op(out),eval(f)];
od;
eval(out);
end;


getkernelmodn:=proc(N,n)
#
# read off a basis for the kernel of N
#
M:=Gaussjord(N) mod n;
out:=[];
lp:=[];
lead:={}; # set of column numbers of leading ones
for i from 1 to rowdim(matrix(M)) do
 j:=1;
 while M[i,j]=0 and j < coldim(M) do
  j:=j+1;
 od;
 if M[i,j]=1 then
  lead:=lead union {eval(j)};
  lp:=[op(lp),eval(j)];
 fi;
od;
for i from 1 to coldim(M) do
 if not eval(i) in lead then
  v:=[];
  for l from 1 to coldim(M) do
   v:=[op(v),0];
  od;
 for z from 1 to coldim(matrix([lp])) do
  v[lp[z]]:=eval(M[z,i]);
  od;
  v[i]:=-1 mod n;
  out:=[op(out),eval(v)];
 fi;
od;
eval(out);
end;


ybcocymodn := proc(B,m)
#
# Computes the YB 2-cocycles of biquandle B mod n
#
j:=1;
z1:=[];
for i from 1 to coldim(evalm(B[1])) do
 for j from 1 to coldim(evalm(B[1])) do
  z1:=[op(z1),0];
 od;
od;
n:=coldim(evalm(B[1]));
E:=[];
for i from 1 to n do
 for j from 1 to n do
  for k from 1 to n do
  p:=eval(z1);
  p[n*(i-1)+j]:= p[n*(i-1)+j]+1 mod m;
  p[n*(eval(B[2][i,j])-1)+k]:=p[n*(eval(B[2][i,j])-1)+k]+1 mod m;
  p[n*(eval(B[4][j,i])-1)+eval(B[4][k,B[2][i,j]])]:=
   p[n*(eval(B[4][j,i])-1)+eval(B[4][k,B[2][i,j]])]+1 mod m;
  p[n*(j-1)+k]:=p[n*(j-1)+k]-1 mod m;
  p[n*(i-1)+eval(B[4][k,j])]:=p[n*(i-1)+eval(B[4][k,j])]-1 mod m;
  p[n*(eval(B[2][i,B[4][k,j]])-1)+eval(B[2][j,k])]:=
   p[n*(eval(B[2][i,B[4][k,j]])-1)+eval(B[2][j,k])]-1 mod m;
  E:=[op(E),eval(p)];
  od;
 od;
od;
eval(getkernelmodn(E,m));
end;


ybcobmodn:=proc(B,n)
#
# Find YB 2-coboundaries mod n
#
out:=[];
for i from 1 to coldim(matrix(B[1])) do
 cob:=[];
 for j from 1 to coldim(matrix(B[1])) do
  for k from 1 to coldim(matrix(B[1])) do
   temp:=0;
   if k=i then 
    temp:=temp+1 mod n;
   fi;
   if j=i then
    temp:=temp+1 mod n;
   fi;
   if B[2][j,k]=i then
    temp:=temp-1 mod n;
   fi;
   if B[4][k,j]=i then
    temp:=temp-1 mod n;
   fi;
   cob:=[op(cob),eval(temp)];
  od;
 od;
 out:=[op(out),eval(cob)];
od;
M2:=matrix(out);
M:=Gaussjord(M2) mod n;
out2:=[];
for i from 1 to rowdim(M) do
 j:=getrow(i,M);
 if not zerow(j) then
  out2:=[op(out2),eval(j)];
 fi;
od;
eval(out2);
end;

ybcohommodn:=proc(B,n)
#
# Computes the Yang-Baxter cohomology with Z_n coefficients
#
Z:=ybcocymodn(B,n);
Bn:=ybcobmodn(B,n);
out:=[];
while coldim(matrix([Z]))<>0 do
 x:=eval(Z[1]);
 Z:=subsop(1=NULL,Z);
 if not inspanmodn(eval(x),Bn,n) then
  out:=[op(out),eval(x)];
  Bn:=[op(Bn),eval(x)]
 fi;
od;
eval(out);
end;

redybcohommodn:=proc(B,n)
#
# list of YB cocycles which satisfy type I condition
#
L:=ybcohommodn(B,n);
out:=[];
for X in L do
 if typeI(X,B) then
  out:=[op(out),eval(X)];
 fi;
od;
eval(out);
end;

ybinvmodn:=proc(v,phi,B,m)
#
#  Compute Yang-Baxter cocyle invariant for Gauss code v
#  where phi is a type I cocycle in C^2(B)
#
v2:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  v2:=[op(v2),eval(v[i])];
 fi;
od;
b:=[]; # position in v2 of corresponding crossing label
for i from 1 to coldim(matrix([v2])) do
 j:=1;
 while v2[i]+v2[j]<>0 do
  j:=j+1;
 od;
 b:=[op(b),eval(j)];
od;
nx:=[];
lt:=1;
c:=0;
for i from 1 to coldim(matrix([v]))-1 do
 if v[i+1]=0 then
  nx:=[op(nx),eval(lt)];
  lt:=eval(i+1-c);
  c:=c+1;
 elif v[i]<>0 then
  nx:=[op(nx),eval(i+1-c)];
 fi;
od;
#C:=bhomlist(gauss2biq(v),B);
C:=bhomlist2(v,B);
n:=coldim(B[1]);
f:=0;
for c in C do
 z:=0;
 for i from 1 to coldim(matrix([v2])) do
  if Re(v2[i])<0 then
   if Im(v2[i]) =0 then 
    z:=z+phi[n*(eval(c[i])-1)+eval(c[b[i]])] mod m;
   elif Im(v2[i])<>0 then
    z:=z-phi[n*(eval(c[nx[i]])-1)+eval(c[nx[b[i]]])] mod m;
   fi;
  fi;
 od;
 f:=f+t^eval(z);
od;
eval(f);
end;

ybinv2modn:=proc(v,B,m)
#
#  Compute Yang-Baxter cocyle invariant for Gauss code v
#  for all type I cocycles in C^2(B)
#
P:=redybcohommodn(B,m);
v2:=[];
for i from 1 to coldim(matrix([v])) do
 if v[i]<>0 then
  v2:=[op(v2),eval(v[i])];
 fi;
od;
b:=[]; # position in v2 of corresponding crossing label
for i from 1 to coldim(matrix([v2])) do
 j:=1;
 while v2[i]+v2[j]<>0 do
  j:=j+1;
 od;
 b:=[op(b),eval(j)];
od;
nx:=[];
lt:=1;
c:=0;
for i from 1 to coldim(matrix([v]))-1 do
 if v[i+1]=0 then
  nx:=[op(nx),eval(lt)];
  lt:=eval(i+1-c);
  c:=c+1;
 elif v[i]<>0 then
  nx:=[op(nx),eval(i+1-c)];
 fi;
od;
#C:=bhomlist(gauss2biq(v),B);
C:=bhomlist2(v,B);
n:=coldim(B[1]);
out:=[];
for phi in P do
 f:=0;
 for c in C do
  z:=0;
  for i from 1 to coldim(matrix([v2])) do
   if Re(v2[i])<0 then
    if Im(v2[i]) =0 then 
     z:=z+phi[n*(eval(c[i])-1)+eval(c[b[i]])];
    elif Im(v2[i])<>0 then
     z:=z-phi[n*(eval(c[nx[i]])-1)+eval(c[nx[b[i]]])];
    fi;
   fi;
  od;
  z:=z mod m;
  f:=f+t^eval(z);
 od;
 out:=[op(out),eval(f)];
od;
eval(out);
end;


inspanmodn:=proc(v,L,n)
#
# Checks whether v is in Span(L) mod n
#
W:=Gaussjord(matrix(L)) mod n;
L2:=[op(L),eval(v)];
W2:=Gaussjord(matrix(L2)) mod n;
d1:=0;
d2:=0;
for i from 1 to rowdim(W) do
 if not zerow(getrow(i,W)) then 
  d1:=d1+1;
 fi;
od;
for i from 1 to rowdim(W2) do
 if not zerow(getrow(i,W2)) then 
  d2:=d2+1;
 fi;
od;
if d1<>d2 then
  eval(false);
 else
  eval(true);
fi;
end;


Reverse:=proc(A)
#
# Reverse of a biquandle
#
n:=coldim(A[1]);
w:=[];
W:=[];
B:=[];
for i from 1 to n do
 w:=[op(w),0];
od;
for i from 1 to n do
 W:=[op(W),eval(w)];
od;
for i from 1 to 4 do
 B:=[op(B),matrix(W)];
od;
 for i from 1 to n do
  for j from 1 to n do
   for i1 from 1 to n do
    for j1 from 1 to n do
     if i1=A[1][i,j] and j1=A[3][j,i] then
      B[1][i1,j1]:=eval(i);
      B[3][j1,i1]:=eval(j);
     fi;
     if i1=A[3][i,j] and j1=A[1][j,i] then
      B[3][i1,j1]:=eval(i);
      B[1][j1,i1]:=eval(j);
     fi;
     if i1=A[2][i,j] and j1=A[4][j,i] then
      B[2][i1,j1]:=eval(i);
      B[4][j1,i1]:=eval(j);
     fi;
     if i1=A[4][i,j] and j1=A[2][j,i] then
      B[4][i1,j1]:=eval(i);
      B[2][j1,i1]:=eval(j);
     fi;
    od;
   od;
  od;
 od;
eval(B);
end;


abq:=proc(n,s,t)
#
# Computes the matrix of Alexander biquandle structure on Z/n 
# with s,t invertible
#
out:=true;
si:=0;
ti:=0;
for i from 1 to n do
 if s*i mod n = 1 then
  si:=i;
 fi;
 if t*i mod n = 1 then
  ti:=i;
 fi;
od;
if si=0 or ti=0 then
 out:=false;
fi;
w:=[];
W:=[];
B:=[];
for i from 1 to n do
 w:=[op(w),0];
od;
for i from 1 to n do
 W:=[op(W),eval(w)];
od;
for i from 1 to 4 do
 B:=[op(B),matrix(W)];
od;
for i from 1 to n do
 for j from 1 to n do
  if ti*i+(1-si*ti)*j mod n =0 then
   B[1][i,j]:=n;
  else
   B[1][i,j]:=ti*i+(1-si*ti)*j mod n;
  fi;
  if t*i+(1-s*t)*j mod n =0 then
   B[2][i,j]:=n;
  else
   B[2][i,j]:=t*i+(1-s*t)*j mod n;
  fi;
  if si*i mod n =0 then
   B[3][i,j]:=n;
  else
   B[3][i,j]:=si*i mod n;
  fi;
  if s*i mod n=0 then
   B[4][i,j]:=n;
  else
   B[4][i,j]:=s*i mod n;
  fi;
 od;
od;
if out then 
 eval(B);
else
 eval(false);
fi;
end;