############################################
#
#  June 15, 2007 -- corrected a bug in 
#  biqtest relating to type I move
#  condition.   ---Sam
#
############################################

read("quandles-maple.txt");

bhomtest := proc (A,B,v)
#
# # Tests whether v:A->B is a biquandle homomorphism
#
ret:=true;
if coldim(matrix([v]))<> coldim(A[1]) then ret:=false; fi;
for i from 1 to coldim(matrix([v])) do
 for j from 1 to coldim(matrix([v])) do
  if A[1][i,j] <> 0 and v[A[1][i,j]] <> B[1][v[i],v[j]] then
   ret:=false; break;
  fi;
 od;
od;
if ret then 
 for i from 1 to coldim(matrix([v])) do
  for j from 1 to coldim(matrix([v])) do
   if A[2][i,j] <> 0 and v[A[2][i,j]] <> B[2][v[i],v[j]] then
    ret:=false; break;
   fi;
  od;
 od;
fi;
if ret then 
 for i from 1 to coldim(matrix([v])) do
  for j from 1 to coldim(matrix([v])) do
   if A[3][i,j] <> 0 and v[A[3][i,j]] <> B[3][v[i],v[j]] then
    ret:=false; break;
   fi;
  od;
 od;
fi;
if ret then 
 for i from 1 to coldim(matrix([v])) do
  for j from 1 to coldim(matrix([v])) do
   if A[4][i,j] <> 0 and v[A[4][i,j]] <> B[4][v[i],v[j]] then
    ret:=false; break;
   fi;
  od;
 od;
fi;
eval(ret);
end;
 

opinv :=proc (A,i,j)
#
# # Finds the row x such that i=A[x,j]  
#
ret:=false;
c:=0;
z:=false;
for x from 1 to coldim(A) do
 if A[x,j]=i then
  ret:=x;
  c:=c+1;
 fi;
 if A[x,j]=0 then
  z:=true;
 fi;
od;
if c>1 then 
 ret:=false;
elif z then
 ret:=true;
fi;
eval(ret);
end;

biqtest := proc(A)
#
# # Tests a vector of four square matrices for the
# # biquandle axioms
#
ret:=true;
n:=coldim(A[1]);
for i from 1 to n do if ret then
 for j from 1 to n do
  if i <> A[1][A[2][i,j],A[4][j,i]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do
  if j <> A[3][A[4][j,i],A[2][i,j]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do
  if i <> A[2][A[1][i,j],A[3][j,i]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do
  if j <> A[4][A[3][j,i],A[1][i,j]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[2][A[2][i,j],k] <>
       A[2][A[2][i,A[4][k,j]],A[2][j,k]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[4][A[4][k,j],i] <>
       A[4][A[4][k,A[2][i,j]],A[4][j,i]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[2][A[4][j,i],A[4][k,A[2][i,j]]] <>
       A[4][A[2][j,k],A[2][i,A[4][k,j]]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[1][A[1][i,j],k] <>
       A[1][A[1][i,A[3][k,j]],A[1][j,k]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[3][A[3][k,j],i] <>
       A[3][A[3][k,A[1][i,j]],A[3][j,i]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[1][A[3][j,i],A[3][k,A[1][i,j]]] <>
       A[3][A[1][j,k],A[1][i,A[3][k,j]]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do
 for j from 1 to n do
  x:=opinv(A[1],i,j);
  if x = false then 
   ret:=false;
  else
   if x<> A[2][i,A[3][j,x]] then 
    ret:=false;
   fi;
   if j <> A[4][A[3][j,x],i] then
    ret:=false;
   fi;
  fi;
 od;
od;
for i from 1 to n do
 for j from 1 to n do
  x:=opinv(A[2],i,j);
  if x = false then 
   ret:=false;
  else
   if x<> A[1][i,A[4][j,x]] then 
    ret:=false;
   fi;
   if j <> A[3][A[4][j,x],i] then
    ret:=false;
   fi;
  fi;
 od;
od;
for i from 1 to n do
 r1c:=0;
 for x from 1 to n do
  if x =A[4][i,x] and i=A[2][x,i] then
   r1c:=r1c+1;
  fi;
 od;
 if r1c<>1 then
  ret:=false;
 break;
 fi;
od;
for i from 1 to n do
 r1c:=0;
 for x from 1 to n do
  if x =A[3][i,x] and i=A[1][x,i] then
   r1c:=r1c+1;
  fi;
 od;
 if r1c <> 1 then
  ret:=false;
 break;
 fi;
od;
eval(ret);
end;


biqtest.old := proc(A)
#
# # Tests a vector of four square matrices for the
# # biquandle axioms
#
ret:=true;
n:=coldim(A[1]);
for i from 1 to n do if ret then
 for j from 1 to n do
  if i <> A[1][A[2][i,j],A[4][j,i]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do
  if j <> A[3][A[4][j,i],A[2][i,j]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do
  if i <> A[2][A[1][i,j],A[3][j,i]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do
  if j <> A[4][A[3][j,i],A[1][i,j]]  then
   ret:=false; 
   break;
  fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[2][A[2][i,j],k] <>
       A[2][A[2][i,A[4][k,j]],A[2][j,k]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[4][A[4][k,j],i] <>
       A[4][A[4][k,A[2][i,j]],A[4][j,i]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[2][A[4][j,i],A[4][k,A[2][i,j]]] <>
       A[4][A[2][j,k],A[2][i,A[4][k,j]]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[1][A[1][i,j],k] <>
       A[1][A[1][i,A[3][k,j]],A[1][j,k]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[3][A[3][k,j],i] <>
       A[3][A[3][k,A[1][i,j]],A[3][j,i]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do if ret then
 for j from 1 to n do if ret then
  for k from 1 to n do
    if A[1][A[3][j,i],A[3][k,A[1][i,j]]] <>
       A[3][A[1][j,k],A[1][i,A[3][k,j]]]
     then ret:=false; 
     break;
    fi;
  od; fi;
 od; fi;
od;
for i from 1 to n do
 for j from 1 to n do
  x:=opinv(A[1],i,j);
  if x<> A[2][i,A[3][j,x]] then 
   ret:=false;
  fi;
  if j <> A[4][A[3][j,x],i] then
   ret:=false;
  fi;
 od;
od;
for i from 1 to n do
 for j from 1 to n do
  x:=opinv(A[2],i,j);
  if x<> A[1][i,A[4][j,x]] then 
   ret:=false;
  fi;
  if j <> A[3][A[4][j,x],i] then
   ret:=false;
  fi;
 od;
od;
temp:=false;
for i from 1 to n do
 for x from 1 to n do
  if x =A[4][i,x] and i=A[2][x,i] then
   temp:=true; 
  fi;
 od;
 if not temp then 
 ret:=false;
 break;
 fi;
od;
temp:=false;
for i from 1 to n do
 for x from 1 to n do
  if x =A[1][i,x] and i=A[3][x,i] then
   temp:=true; 
  fi;
 od;
 if not temp then 
 ret:=false;
 break;
 fi;
od;
eval(ret);
end;

bfindzero:=proc(A)
#
# find the best zero entry in the biquandle
#
F:=false;
r:=0;
for k from 1 to 4 do 
 for i from 1 to coldim(A[1]) do 
  for j from 1 to coldim(A[1]) do
   if A[k][i,j]=0 then
    if not F then
     ret:=[k,i,j];
     F:=true;
    else
     c:=ratezero(A,k,i,j);
     if c > r then
      r:=eval(c);
      ret:=[k,i,j];
     fi;
    fi;
   fi;
  od; 
 od;
od;
if F then 
 eval(ret);
else
 eval(F)
fi;
end;


ratezero:=proc(A,m,i,j)
#
# rates a zero entry; higher number means more word completions
#
E:=0;
if m=1 then
 for k from 1 to coldim(A[1]) do
   if A[1][j,k]<>0 and A[3][k,j]<>0  then
    E:=E+1;
   fi;
   if A[3][i,k]<>0 and A[1][j,k]<>0  then
    E:=E+1;
   fi;
   if A[2][i,j]<>0 and A[4][j,i]<>0  then
    E:=E+1;
   fi;
   if A[2][i,k]<>0 and A[4][k,i]<>0  then
    E:=E+1;
   fi;
   if A[3][j,k]<>0 and A[4][k,i]<>0   then
    E:=E+1;
   fi;
   if A[2][k,j]<>0 and A[4][j,k]<>0  then
    E:=E+1;
   fi;
   if A[3][j,i]<>0  then
    E:=E+1;
   fi;
 od;
elif m=2 then 
 for k from 1 to coldim(A[1]) do
   if A[2][j,k]<>0 and A[4][k,j]<>0  then
    E:=E+1;
   fi;
   if A[1][k,j]<>0 and A[3][j,k]<>0  then
    E:=E+1;
   fi;
   if A[4][j,i]<>0 then
    E:=E+1;
   fi;
 od;
elif m=3 then 
 for k from 1 to coldim(A[1]) do
   if A[3][j,k]<>0 and A[1][k,j]<>0 then
    E:=E+1;
   fi;
   if A[4][k,j]<>0 and A[2][j,k]<>0 then
    E:=E+1;
   fi;
   if A[4][k,j]<>0 and A[2][j,k]<>0 then
    E:=E+1;
   fi;
   if A[1][j,i]<>0 then
    E:=E+1;
   fi;
 od;
elif m=4 then 
 for k from 1 to coldim(A[1]) do
   if A[2][j,i]<>0 then
    E:=E+1;
   fi;
   if A[3][i,k]<>0 and A[1][k,i]<>0  then
    E:=E+1;
   fi;
   if A[3][k,j]<>0 and A[1][j,k]<>0  then
    E:=E+1;
   fi;
 od;
fi;
rc:=0; 
cc:=0;
for k from 1 to coldim(A[1]) do
  if A[m][i,k]=0 then
   rc:=rc+1;
  fi;
  if A[m][k,j]=0 then
   cc:=cc+1;
  fi; 
od;
if rc=2 
 then E:=E+1; 
fi;
if cc=2 
 then E:=E+1; 
fi;
eval(E);
end;

avail := proc(M,r,c,m)
#
# Returns the smallest available entry for position (r,c) in M
# greater than or equal to m or "false" if none is availible
#
n:=coldim(M);
ret:=false;
A4:=false;
for i from 1 to n do
 if M[r,i]=i then 
  A4:=true;
  break;
 fi;
od;
for i from m to n do
 jei:=false;
 for j from 1 to n do
  if M[j,c]=i then
   jei:=true;
   break;
  fi;
 od;
 if not jei then
  if not (A4 and i=c) then
   ret:=i;
   break;
  fi;
 fi;
od;
eval(ret);
end;

permtest := proc (v)
#
# # tests whether a vector represents a permutation
#
ret:=true;
for i from 1 to coldim(matrix([v])) do
 for j from 1 to i-1 do
  if v[i]=v[j] and v[i]<>0 and v[j]<>0 then 
   ret:=false; 
  fi;
 od;
od;
#
ret;
end;


bpermtest:=proc (A)
#
# # tests whether all columns are permutations
#
ret:=true;
for k from 1 to 4 do
 for i from 1 to coldim(A[1]) do
  if not permtest(col(A[k],i)) then
   ret:=false;
  fi;
 od;
od;
eval(ret);
end;


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

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

bhomcount :=proc(A,B)
#
# Computes the biquandle counting invariant
#
eval(coldim(matrix([bhomlist(A,B)])));
end;


Obv:=proc(A)
#
# Obverse of a biquandle
#
eval([A[2],A[1],A[4],A[3]]);
end;


Flip:=proc(A)
#
# Flip of a biquandle
#
eval([A[3],A[4],A[1],A[2]]);
end;


biqlist := proc (A)
#
# find all biquandles which match starting pattern
#
n:=coldim(A[1]);
w:=[[evalm(A[1]),evalm(A[2]),evalm(A[3]),evalm(A[4])]];
out:=[];
cont:=true;
while cont do
 cont:=false;
 for i from 1 to coldim(matrix([w])) do
  B:=[evalm(w[1][1]),evalm(w[1][2]),evalm(w[1][3]),evalm(w[1][4])];
  f:=eval(bfindzero(B));
  w:=subsop(1=NULL,w);
  if not type(f,list) then 
   out:=[op(out),[evalm(B[1]),evalm(B[2]),evalm(B[3]),evalm(B[4])]];
  elif type(f,list) then
   j:=avail(B[f[1]],f[2],f[3],1);
   while type(j,integer) do
    tz:=j+1;
    X:=[evalm(B[1]),evalm(B[2]),evalm(B[3]),evalm(B[4])];
    X[f[1]][f[2],f[3]]:=eval(j);
    if bpermtest(X) then
     Z:=eval(biqfill(X));
     if type(Z[1],matrix) and type(Z[2],matrix) 
         and type(Z[3],matrix) and type(Z[4],matrix) then
      w:=[op(w),[evalm(Z[1]),evalm(Z[2]),evalm(Z[3]),evalm(Z[4])]];
      cont:=true;
     fi;
    fi;
    j:=avail(X[f[1]],f[2],f[3],tz);
   od;
  fi;
 od;
od;
eval(out);
end;


bisolist:= proc(A,B)
#
# lists isomorphisms from A to B
#
M:=bhomlist(A,B);
out:=[];
for i from 1 to coldim(matrix([M])) do
 if permtest(M[i]) then
  out:=[op(out),eval(M[i])];
 fi;
od;
eval(out);
end;

baut:=proc(A)
#
# Compute automorpishm group of A
#
eval(bisolist(A,A));
end;


breducelist:=proc(M)
#
# Removes isomorphic copies
#
X:=[eval(M[1])];
for i from 1 to coldim(matrix([M])) do
 isoc:=false;
 for j from 1 to coldim(matrix([X])) do
  if coldim(matrix([bisolist(M[i],X[j])]))<>0 
  or coldim(matrix([bisolist(M[i],Obv(X[j]))]))<>0 
  or coldim(matrix([bisolist(M[i],Flip(X[j]))]))<>0 
  or coldim(matrix([bisolist(M[i],Obv(Flip(X[j])))]))<>0 
 then
   isoc:=true;
  fi;
 od; 
 if not isoc then
  X:=[op(X),eval(M[i])];
 fi;
od;
eval(X);
end;

sfliptest:=proc(A)
#
# Checks whether a biquandle is self-flip
#
if coldim(matrix([bisolist(A,Flip(A))]))<>0 then
 eval(true);
else
 eval(false);
fi;
end;

sobvtest:=proc(A)
#
# Checks whether a biquandle is self-obverse
#
if coldim(matrix([bisolist(A,Obv(A))]))<>0 then
 eval(true);
else
 eval(false);
fi;
end;

biqfill := proc(B)
#
# Uses biquandle axioms to fill in zeroes
#
A:=[evalm(B[1]),evalm(B[2]),evalm(B[3]),evalm(B[4])];
n:=coldim(A[1]);
cont:=true;      # continue counter
C:=false;        # contradiction
while cont and not C do
 cont:=false;
 for i from 1 to n do
  for j from 1 to n do
   for k from 1 to n do
    if A[2][i,j]<>0 and A[2][j,k]<>0 and A[4][k,j]<>0 and 
      A[2][i,A[4][k,j]]<>0 and
      A[2][A[2][i,j],k]<>A[2][A[2][i,A[4][k,j]],A[2][j,k]] then 
     if A[2][A[2][i,j],k]<>0 and A[2][A[2][i,A[4][k,j]],A[2][j,k]]<>0 then
      C:=true;
     elif A[2][A[2][i,j],k]=0 and A[2][A[2][i,A[4][k,j]],A[2][j,k]]<>0 then
      A[2][A[2][i,j],k]:=eval(A[2][A[2][i,A[4][k,j]],A[2][j,k]]);
      cont:=true;
     elif A[2][A[2][i,j],k]<>0 and A[2][A[2][i,A[4][k,j]],A[2][j,k]]=0 then
      A[2][A[2][i,A[4][k,j]],A[2][j,k]]:=eval(A[2][A[2][i,j],k]);
      cont:=true;
     fi;
    fi;
 if not C then
    if A[4][k,j]<>0 and A[4][j,i]<>0 and A[2][i,j]<>0 and 
       A[4][k,A[2][i,j]]<>0 and
      A[4][A[4][k,j],i]<>A[4][A[4][k,A[2][i,j]],A[4][j,i]] then
     if A[4][A[4][k,j],i]<>0 and A[4][A[4][k,A[2][i,j]],A[4][j,i]]<>0 then
      C:=true;
     elif A[4][A[4][k,j],i]=0 and A[4][A[4][k,A[2][i,j]],A[4][j,i]]<>0 then
      A[4][A[4][k,j],i]:=eval(A[4][A[4][k,A[2][i,j]],A[4][j,i]]);
      cont:=true;
     elif A[4][A[4][k,j],i]<>0 and A[4][A[4][k,A[2][i,j]],A[4][j,i]]=0 then
      A[4][A[4][k,A[2][i,j]],A[4][j,i]]:=eval(A[4][A[4][k,j],i]);
      cont:=true;
     fi;
    fi; 
 fi;
 if not C then
    if A[4][j,i]<>0 and A[2][i,j]<>0 and A[4][k,A[2][i,j]]<>0 and 
       A[4][k,j]<>0 and
      A[2][j,k]<>0 and A[2][i,A[4][k,j]]<>0 and
      A[2][A[4][j,i],A[4][k,A[2][i,j]]]<>
      A[4][A[2][j,k],A[2][i,A[4][k,j]]] then
     if A[2][A[4][j,i],A[4][k,A[2][i,j]]]<>0 and
        A[4][A[2][j,k],A[2][i,A[4][k,j]]]<>0 then
      C:=true;
     elif A[2][A[4][j,i],A[4][k,A[2][i,j]]] =0 and
          A[4][A[2][j,k],A[2][i,A[4][k,j]]]<>0 then
      A[2][A[4][j,i],A[4][k,A[2][i,j]]]:=
      eval(A[4][A[2][j,k],A[2][i,A[4][k,j]]]);
      cont:=true;
     elif A[2][A[4][j,i],A[4][k,A[2][i,j]]]<>0 and
          A[4][A[2][j,k],A[2][i,A[4][k,j]]]=0 then
       A[4][A[2][j,k],A[2][i,A[4][k,j]]]:=
       eval(A[2][A[4][j,i],A[4][k,A[2][i,j]]]);
       cont:=true;
     fi;
    fi;
 fi;
 if not C then
    if A[1][i,j]<>0 and A[1][j,k]<>0 and A[3][k,j]<>0 and 
      A[1][i,A[3][k,j]]<>0 and
      A[1][A[1][i,j],k] <> A[1][A[1][i,A[3][k,j]],A[1][j,k]] then
     if A[1][A[1][i,j],k]<>0 and  A[1][A[1][i,A[3][k,j]],A[1][j,k]]<>0 then
      C:=true;
     elif A[1][A[1][i,j],k]=0 and A[1][A[1][i,A[3][k,j]],A[1][j,k]]<>0 then
      A[1][A[1][i,j],k]:=eval(A[1][A[1][i,A[3][k,j]],A[1][j,k]]);
      cont:=true;
     elif  A[1][A[1][i,j],k]<>0 and A[1][A[1][i,A[3][k,j]],A[1][j,k]]=0 then
      A[1][A[1][i,A[3][k,j]],A[1][j,k]]:=eval(A[1][A[1][i,j],k]);
      cont:=true;
     fi;
    fi;
 fi;
 if not C then
   if A[3][k,j]<>0 and A[3][j,i]<>0 and A[1][i,j]<>0 and 
      A[3][k,A[1][i,j]]<>0 and
      A[3][A[3][k,j],i] <> A[3][A[3][k,A[1][i,j]],A[3][j,i]] then
    if A[3][A[3][k,j],i]<>0 and A[3][A[3][k,A[1][i,j]],A[3][j,i]]<>0 then
     C:=true;
    elif A[3][A[3][k,j],i]=0 and A[3][A[3][k,A[1][i,j]],A[3][j,i]]<>0 then
     A[3][A[3][k,j],i]:=eval(A[3][A[3][k,A[1][i,j]],A[3][j,i]]);
     cont:=true;
    elif A[3][A[3][k,j],i]<>0 and A[3][A[3][k,A[1][i,j]],A[3][j,i]]=0 then
     A[3][A[3][k,A[1][i,j]],A[3][j,i]]:=eval(A[3][A[3][k,j],i]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[3][j,i]<>0 and A[1][i,j]<>0 and A[3][k,A[1][i,j]]<>0 and
      A[1][j,k]<>0 and A[3][k,j]<>0 and A[1][i,A[3][k,j]] <>0 and
      A[1][A[3][j,i],A[3][k,A[1][i,j]]] <> 
      A[3][A[1][j,k],A[1][i,A[3][k,j]]] then
    if A[1][A[3][j,i],A[3][k,A[1][i,j]]] <> 0 and
       A[3][A[1][j,k],A[1][i,A[3][k,j]]] <> 0 then
     C:=true;
    elif  A[1][A[3][j,i],A[3][k,A[1][i,j]]]=0 and
       A[3][A[1][j,k],A[1][i,A[3][k,j]]] <> 0 then
     A[1][A[3][j,i],A[3][k,A[1][i,j]]]:=
     eval(A[3][A[1][j,k],A[1][i,A[3][k,j]]]);
     cont:=true;
    elif  A[1][A[3][j,i],A[3][k,A[1][i,j]]]<>0 and
       A[3][A[1][j,k],A[1][i,A[3][k,j]]]=0 then
     A[3][A[1][j,k],A[1][i,A[3][k,j]]]:=
     eval(A[1][A[3][j,i],A[3][k,A[1][i,j]]]);
     cont:=true;
    fi;
   fi;
  fi;
# From RIII moves with other orientations
 if not C then
   if A[2][k,i]<>0 and A[1][j,i]<>0 and A[3][i,j]<>0 and
      A[2][k,A[1][j,i]] <>0 and 
      A[2][A[2][k,i],j]<> A[2][A[2][k,A[1][j,i]],A[3][i,j]]then
    if A[2][A[2][k,i],j] <> 0 and
      A[2][A[2][k,A[1][j,i]],A[3][i,j]] <> 0 then
     C:=true;
    elif A[2][A[2][k,i],j] =0 and
       A[2][A[2][k,A[1][j,i]],A[3][i,j]] <> 0 then
     A[2][A[2][k,i],j]:=
     eval(A[2][A[2][k,A[1][j,i]],A[3][i,j]]);
     cont:=true;
    elif  A[2][A[2][k,i],j]<>0 and
      A[2][A[2][k,A[1][j,i]],A[3][i,j]] =0 then
     A[2][A[2][k,A[1][j,i]],A[3][i,j]]:=
     eval(A[2][A[2][k,i],j]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[3][i,j]<>0 and A[1][j,i]<>0 and A[2][k,A[1][j,i]]<>0 and
     A[4][i,k]<>0 and A[2][k,i]<>0 and A[4][j,A[2][k,i]]<>0 and
     A[4][A[3][i,j],A[2][k,A[1][j,i]]]<> A[3][A[4][i,k],A[4][j,A[2][k,i]]] then
    if A[4][A[3][i,j],A[2][k,A[1][j,i]]]<>0 and
       A[3][A[4][i,k],A[4][j,A[2][k,i]]]<>0 then
     C:=true;
    elif A[4][A[3][i,j],A[2][k,A[1][j,i]]] =0 and
        A[3][A[4][i,k],A[4][j,A[2][k,i]]]<> 0 then
     A[4][A[3][i,j],A[2][k,A[1][j,i]]]:=
     eval(A[3][A[4][i,k],A[4][j,A[2][k,i]]]);
     cont:=true;
    elif A[4][A[3][i,j],A[2][k,A[1][j,i]]] <>0 and
       A[3][A[4][i,k],A[4][j,A[2][k,i]]]=0 then
     A[3][A[4][i,k],A[4][j,A[2][k,i]]]:=
     eval(A[4][A[3][i,j],A[2][k,A[1][j,i]]]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[1][j,i]<>0 and A[2][k,i]<>0 and A[4][i,k]<>0 and
      A[4][j,A[2][k,i]]<>0 and 
      A[4][A[1][j,i],k]<> A[1][A[4][j,A[2][k,i]],A[4][i,k]]  then
    if A[4][A[1][j,i],k] <> 0 and
       A[1][A[4][j,A[2][k,i]],A[4][i,k]] <> 0 then
     C:=true;
    elif A[4][A[1][j,i],k] =0 and
      A[1][A[4][j,A[2][k,i]],A[4][i,k]]  <> 0 then
     A[4][A[1][j,i],k]:=
     eval(A[1][A[4][j,A[2][k,i]],A[4][i,k]]);
     cont:=true;
    elif A[4][A[1][j,i],k] <>0 and
       A[1][A[4][j,A[2][k,i]],A[4][i,k]]=0 then
     A[1][A[4][j,A[2][k,i]],A[4][i,k]]:=
     eval(A[4][A[1][j,i],k]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[4][k,i]<>0 and A[3][j,i]<>0 and A[1][i,j]<>0 and
      A[4][k,A[3][j,i]]<>0 and
      A[4][A[4][k,i],j]<>A[4][A[4][k,A[3][j,i]],A[1][i,j]]  then
    if A[4][A[4][k,i],j] <> 0 and
       A[4][A[4][k,A[3][j,i]],A[1][i,j]] <> 0 then
     C:=true;
    elif A[4][A[4][k,i],j] =0 and
        A[4][A[4][k,A[3][j,i]],A[1][i,j]]<> 0 then
     A[4][A[4][k,i],j]:=
     eval(A[4][A[4][k,A[3][j,i]],A[1][i,j]]);
     cont:=true;
    elif A[4][A[4][k,i],j] <>0 and
      A[4][A[4][k,A[3][j,i]],A[1][i,j]] =0 then
     A[4][A[4][k,A[3][j,i]],A[1][i,j]]:=
     eval(A[4][A[4][k,i],j]);
     cont:=true;
   fi;
   fi;
 fi;
 if not C then
   if A[2][i,k]<>0 and A[4][k,i]<>0 and A[2][j,A[4][k,i]]<>0 and
      A[1][i,j]<>0 and A[3][j,i]<>0 and A[4][k,A[3][j,i]]<>0 and
      A[1][A[2][i,k],A[2][j,A[4][k,i]]]<>A[2][A[1][i,j],A[4][k,A[3][j,i]]] then
    if A[1][A[2][i,k],A[2][j,A[4][k,i]]] <> 0 and
       A[2][A[1][i,j],A[4][k,A[3][j,i]]] <> 0 then
     C:=true;
    elif A[1][A[2][i,k],A[2][j,A[4][k,i]]] =0 and
        A[2][A[1][i,j],A[4][k,A[3][j,i]]]<> 0 then
     A[1][A[2][i,k],A[2][j,A[4][k,i]]]:=
     eval(A[2][A[1][i,j],A[4][k,A[3][j,i]]]);
     cont:=true;
    elif A[1][A[2][i,k],A[2][j,A[4][k,i]]] <>0 and
      A[2][A[1][i,j],A[4][k,A[3][j,i]]] =0 then
     A[2][A[1][i,j],A[4][k,A[3][j,i]]]:=
     eval(A[1][A[2][i,k],A[2][j,A[4][k,i]]]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[4][k,i]<>0 and A[2][i,k]<>0 and A[2][j,A[4][k,i]]<>0 and
      A[3][j,i]<>0 and
      A[3][A[2][j,A[4][k,i]],A[2][i,k]]<>A[2][A[3][j,i],k]  then
    if A[3][A[2][j,A[4][k,i]],A[2][i,k]] <> 0 and
       A[2][A[3][j,i],k] <> 0 then
     C:=true;
    elif A[3][A[2][j,A[4][k,i]],A[2][i,k]] =0 and
      A[2][A[3][j,i],k]  <> 0 then
     A[3][A[2][j,A[4][k,i]],A[2][i,k]]:=
     eval(A[2][A[3][j,i],k]);
     cont:=true;
    elif A[3][A[2][j,A[4][k,i]],A[2][i,k]] <>0 and
      A[2][A[3][j,i],k] =0 then
     A[2][A[3][j,i],k]:=
     eval(A[3][A[2][j,A[4][k,i]],A[2][i,k]]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[3][k,j]<>0 and A[4][i,j]<>0 and A[2][j,i]<>0 and
      A[3][k,A[4][i,j]]<>0 and
      A[3][A[3][k,j],i]<>A[3][A[3][k,A[4][i,j]],A[2][j,i]]  then
    if A[3][A[3][k,j],i] <> 0 and
       A[3][A[3][k,A[4][i,j]],A[2][j,i]] <> 0 then
     C:=true;
    elif A[3][A[3][k,j],i] =0 and
        A[3][A[3][k,A[4][i,j]],A[2][j,i]]<> 0 then
     A[3][A[3][k,j],i]:=
     eval(A[3][A[3][k,A[4][i,j]],A[2][j,i]]);
     cont:=true;
    elif A[3][A[3][k,j],i] <>0 and
      A[3][A[3][k,A[4][i,j]],A[2][j,i]] =0 then
     A[3][A[3][k,A[4][i,j]],A[2][j,i]]:=
     eval(A[3][A[3][k,j],i]);
    cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[1][j,k]<>0 and A[3][k,j]<>0 and A[1][i,A[3][k,j]]<>0 and
      A[2][j,i]<>0 and A[4][i,j]<>0 and A[3][k,A[4][i,j]]<>0 and
      A[2][A[1][j,k],A[1][i,A[3][k,j]]]<>A[1][A[2][j,i],A[3][k,A[4][i,j]]] then
    if A[2][A[1][j,k],A[1][i,A[3][k,j]]] <> 0 and
       A[1][A[2][j,i],A[3][k,A[4][i,j]]] <> 0 then
     C:=true;
    elif A[2][A[1][j,k],A[1][i,A[3][k,j]]] =0 and
      A[1][A[2][j,i],A[3][k,A[4][i,j]]]  <> 0 then
     A[2][A[1][j,k],A[1][i,A[3][k,j]]]:=
     eval(A[1][A[2][j,i],A[3][k,A[4][i,j]]]);
     cont:=true;
    elif A[2][A[1][j,k],A[1][i,A[3][k,j]]] <>0 and
      A[1][A[2][j,i],A[3][k,A[4][i,j]]] =0 then
     A[1][A[2][j,i],A[3][k,A[4][i,j]]]:=
     eval(A[2][A[1][j,k],A[1][i,A[3][k,j]]]);
    cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[3][k,j]<>0 and A[1][j,k]<>0 and A[1][i,A[3][k,j]]<>0 and
     A[4][i,j]<>0 and
     A[4][A[1][i,A[3][k,j]],A[1][j,k]] <> A[1][A[4][i,j],k] then
    if A[4][A[1][i,A[3][k,j]],A[1][j,k]] <> 0 and
       A[1][A[4][i,j],k] <> 0 then
     C:=true;
    elif A[4][A[1][i,A[3][k,j]],A[1][j,k]] =0 and
      A[1][A[4][i,j],k]  <> 0 then
     A[4][A[1][i,A[3][k,j]],A[1][j,k]]:=
     eval(A[1][A[4][i,j],k]);
     cont:=true;
    elif A[4][A[1][i,A[3][k,j]],A[1][j,k]] <>0 and
      A[1][A[4][i,j],k] =0 then
      A[1][A[4][i,j],k]:=
     eval(A[4][A[1][i,A[3][k,j]],A[1][j,k]]);
    cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[1][i,j]<>0 and A[3][j,i]<>0 and A[3][k,A[1][i,j]]<>0 and
      A[2][k,j]<>0 and
      A[2][A[3][k,A[1][i,j]],A[3][j,i]]<>A[3][A[2][k,j],i]  then
    if A[2][A[3][k,A[1][i,j]],A[3][j,i]] <> 0 and
       A[3][A[2][k,j],i] <> 0 then
     C:=true;
    elif A[2][A[3][k,A[1][i,j]],A[3][j,i]] =0 and
       A[3][A[2][k,j],i] <> 0 then
     A[2][A[3][k,A[1][i,j]],A[3][j,i]]:=
     eval(A[3][A[2][k,j],i]);
     cont:=true;
    elif A[2][A[3][k,A[1][i,j]],A[3][j,i]] <>0 and
      A[3][A[2][k,j],i] =0 then
     A[3][A[2][k,j],i]:=
     eval(A[2][A[3][k,A[1][i,j]],A[3][j,i]]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[3][j,i]<>0 and A[1][i,j]<>0 and A[3][k,A[1][i,j]]<>0 and
      A[4][j,k]<>0 and A[2][k,j]<>0 and A[1][i,A[2][k,j]]<>0 and
      A[4][A[3][j,i],A[3][k,A[1][i,j]]]<>A[3][A[4][j,k],A[1][i,A[2][k,j]]] then
    if A[4][A[3][j,i],A[3][k,A[1][i,j]]] <> 0 and
       A[3][A[4][j,k],A[1][i,A[2][k,j]]] <> 0 then
     C:=true;
    elif A[4][A[3][j,i],A[3][k,A[1][i,j]]] =0 and
       A[3][A[4][j,k],A[1][i,A[2][k,j]]] <> 0 then
       A[4][A[3][j,i],A[3][k,A[1][i,j]]]:=
     eval(A[3][A[4][j,k],A[1][i,A[2][k,j]]]);
     cont:=true;
    elif A[4][A[3][j,i],A[3][k,A[1][i,j]]] <>0 and
      A[3][A[4][j,k],A[1][i,A[2][k,j]]] =0 then
      A[3][A[4][j,k],A[1][i,A[2][k,j]]]:=
     eval(A[4][A[3][j,i],A[3][k,A[1][i,j]]]);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[1][i,j]<>0 and A[2][k,j]<>0 and A[4][j,k]<>0 and 
      A[1][i,A[2][k,j]]<>0 and
      A[1][A[1][i,j],k]<>A[1][A[1][i,A[2][k,j]],A[4][j,k]]  then
    if A[1][A[1][i,j],k] <> 0 and
       A[1][A[1][i,A[2][k,j]],A[4][j,k]] <> 0 then
     C:=true;
    elif A[1][A[1][i,j],k] =0 and
      A[1][A[1][i,A[2][k,j]],A[4][j,k]]  <> 0 then
     A[1][A[1][i,j],k]:=
     eval(A[1][A[1][i,A[2][k,j]],A[4][j,k]]);
     cont:=true;
    elif A[1][A[1][i,j],k] <>0 and
     A[1][A[1][i,A[2][k,j]],A[4][j,k]]  =0 then
     A[1][A[1][i,A[2][k,j]],A[4][j,k]]:=
     eval(A[1][A[1][i,j],k]);
     cont:=true;
    fi;
   fi; fi;
   od;
  od;
 od; #fi;
# From axiom II
 for i from 1 to n do
  for j from 1 to n do
   if not C then
   if A[2][i,j]<>0 and A[4][j,i]<>0 then
    if i <> A[1][A[2][i,j],A[4][j,i]] and 
       A[1][A[2][i,j],A[4][j,i]]<> 0 then
     C:=true;
    elif A[1][A[2][i,j],A[4][j,i]]=0 then
     A[1][A[2][i,j],A[4][j,i]]:=eval(i);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[4][j,i]<>0 and A[2][i,j]<>0 then
    if j <> A[3][A[4][j,i],A[2][i,j]] and 
       A[3][A[4][j,i],A[2][i,j]]<>0 then
     C:=true;
    elif A[3][A[4][j,i],A[2][i,j]]=0 then
     A[3][A[4][j,i],A[2][i,j]]:=eval(j);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[1][i,j]<>0 and A[3][j,i]<>0 then
    if i<>A[2][A[1][i,j],A[3][j,i]] and 
       A[2][A[1][i,j],A[3][j,i]]<>0 then
     C:=true;
    elif A[2][A[1][i,j],A[3][j,i]]=0 then
     A[2][A[1][i,j],A[3][j,i]]:=eval(i);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if A[3][j,i]<>0 and A[1][i,j]<>0 then
    if j <> A[4][A[3][j,i],A[1][i,j]] and 
       A[4][A[3][j,i],A[1][i,j]]<>0 then
     C:=true;
    elif A[4][A[3][j,i],A[1][i,j]]=0 then
     A[4][A[3][j,i],A[1][i,j]]:=eval(j);
     cont:=true;
    fi;
   fi;
 fi;
 if not C then
   if type(opinv(A[1],i,j),integer) then
    x:=opinv(A[1],i,j);
    if A[3][j,x]<>0 then
     if x<>A[2][i,A[3][j,x]] and A[2][i,A[3][j,x]]<>0 then
      C:=true;
     elif A[2][i,A[3][j,x]]=0 then
      A[2][i,A[3][j,x]]:=eval(x);
      cont:=true;
     fi;
     if j<>A[4][A[3][j,x],i] and A[4][A[3][j,x],i]<>0 then 
      C:=true;
     elif A[4][A[3][j,x],i]=0 then
      A[4][A[3][j,x],i]:=eval(j);
      cont:=true;
     fi;
    fi;
   elif not opinv(A[1],i,j) then
    C:=true;
   fi;
 fi;
 if not C then
   if type(opinv(A[2],i,j),integer) then
    y:=opinv(A[2],i,j);
    if A[4][j,y]<>0 then
     if y<>A[1][i,A[4][j,y]] and A[1][i,A[4][j,y]]<>0 then
      C:=true;
     elif A[1][i,A[4][j,y]]=0 then
      A[1][i,A[4][j,y]]:=eval(y);
      cont:=true;
     fi;
     if j<>A[3][A[4][j,y],i] and A[3][A[4][j,y],i]<>0 then
      C:=true;
     elif A[3][A[4][j,y],i]=0 then
      A[3][A[4][j,y],i]:=eval(j);
      cont:=true;
     fi;
    fi;
   elif not opinv(A[2],i,j) then
    C:=true;
   fi;
   fi;
  od;
 od; #fi;
 if not C then
  for i from 1 to n do
   c4:=0;
   ch4:=false;
   for j from 1 to n do
    if A[1][i,j]=0 then
     ch4:=true;
    elif A[1][i,j]=j then 
     c4:=c4+1;
     if A[3][j,i]<>i then
      if A[3][j,i]=0 then
       A[3][j,i]:=i;
       cont:=true;
      else
       C:=true;
      fi;
     fi;
    fi;
   od;
   if c4<>1 and not ch4 then
    C:=true;
   fi;
  od;
 fi;
 if not C then
  for i from 1 to n do
   c4:=0;
   ch4:=false;
   for j from 1 to n do
    if A[2][i,j]=0 then
     ch4:=true;
    elif A[2][i,j]=j then 
     c4:=c4+1;
     if A[4][j,i]<>i then
      if A[4][j,i]=0 then
       A[4][j,i]:=i;
       cont:=true;
      else
       C:=true;
      fi;
     fi;
    fi;
   od;
   if c4<>1 and not ch4 then
    C:=true;
   fi;
  od;
 fi;
 if not C then
  for i from 1 to n do
   c4:=0;
   ch4:=false;
   for j from 1 to n do
    if A[3][i,j]=0 then
     ch4:=true;
    elif A[3][i,j]=j then 
     c4:=c4+1;
     if A[1][j,i]<>i then
      if A[1][j,i]=0 then
       A[1][j,i]:=i;
       cont:=true;
      else
       C:=true;
      fi;
     fi;
    fi;
   od;
   if c4<>1 and not ch4 then
    C:=true;
   fi;
  od;
 fi;
 if not C then
  for i from 1 to n do
   c4:=0;
   ch4:=false;
   for j from 1 to n do
    if A[4][i,j]=0 then
     ch4:=true;
    elif A[4][i,j]=j then 
     c4:=c4+1;
     if A[2][j,i]<>i then
      if A[2][j,i]=0 then
       A[2][j,i]:=i;
       cont:=true;
      else
       C:=true;
      fi;
     fi;
    fi;
   od;
   if c4<>1 and not ch4 then
    C:=true;
   fi;
  od;
 fi;
# fill in any obvious entries and detect subtle contradiction
 if not C then for i from 1 to 4 do
  if not C then for j from 1 to n do
   tct:=0;
   if not C then for k from 1 to n do
    if A[i][k,j]=0 then
     p:=k;
     tct:=tct+1;  
    fi;
   od; fi;
   if tct=1 and type(avail(A[i],p,j,1),integer) then
    A[i][p,j]:=avail(A[i],p,j,1);
    cont:=true;
   elif tct=1 and type(avail(A[i],p,j,1),boolean) then
    C:=true;
   fi;
  od; fi;
 od; fi;
od;
if C then 
 eval(not C);
else 
 eval(A); 
fi;
end;