/*

1. The command FindRankFq(p,n,m) returns a sequece consists of the Q-rank, Fq-rank where q<=p and q is a prime for the group B2(Zn+Zm).

2.The command Table(p,n) returns a table consists of the Q-rank, Fq-rank where q<=p and q is a prime for all groups B2(Zm+Zk) where m,k<=n.

3.The command Table2(p,n) returns a table consists of the Q-rank, Fq-rank where q<=p and q is a prime for all groups of rank 2 i.e. B2(Zm+Zm) where m<=n.

4.The command TableCyclic(p,n) returns a table consists of the Q-rank, Fq-rank where q<=p and q is a prime for all cyclic groups B2(Zm) where m<=n.

*/

// Create the function checking generation condition
function CheckGenerators(n,m,i1,j1,i2,j2);
// For cyclic case (a,b) generates Zp iff gcd(a,b,p)=1
if n eq 1 then

if Gcd(Gcd(j1,j2), m) eq 1 then
return true;
else 
return false;
end if; 

else

Z:=RModule(IntegerRing(),2);
L:=sub < Z|[n,0],[0,m]>;
Znm:=quo<Z|L>;
gen:=sub<Znm|[i1,j1],[i2,j2]>;

if gen eq Znm then
return true;
else
return false;
end if;

end if;

end function;

//Collect all (a,b) satisfying lexicographic order and generation condition
function B2Generators(n,m)
number:=0;
A:= AssociativeArray();
B:= AssociativeArray();
Index:=AssociativeArray(Universe([[1,2],[3,4]]));

for i1 in [0..n-1] do
for j1 in [0..m-1] do
for i2 in [i1..n-1] do
if i1 lt i2 then
sub:=0;
else
sub:=j1;
end if;
for j2 in [sub..m-1] do

if CheckGenerators(n,m,i1,j1,i2,j2) then
number:=number+1;
A[number] := [i1, j1];
B[number] := [i2, j2];
Index[[i1,j1,i2,j2]]:=number;
end if;

end for;
end for;
end for;
end for;

return A,B,Index,number;
end function;

function FindRankFq(q,n,m);
// Make sure the group is presented in the form compatible with quotient modules in Magma
prod:=n*m;
n:=Gcd(n,m);
m:=prod div n;

// Collect all symbols satisfying lexicographic order 
A,B,Index,number:=B2Generators(n,m);

//Build relation matrix
row:=number*3;
a:=SparseMatrix(RationalField(), row, number);
currentrow:=0;

for i in [1..number] do

//(B) relation
if A[i] ne B[i] then
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]-1;
i1:=(A[i][1]-B[i][1]) mod n;
j1:=(A[i][2]-B[i][2]) mod m;
search:=Sort([[i1,j1],[B[i][1],B[i][2]]]);//(a-b,b)
search:=[search[1][1],search[1][2],search[2][1],search[2][2]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]+ 1;

i2:=(B[i][1]-A[i][1]) mod n;
j2:=(B[i][2]-A[i][2]) mod m;
search:=Sort([[A[i][1],A[i][2]],[i2,j2]]);//(a,b-a)
search:=[search[1][1],search[1][2],search[2][1],search[2][2]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]+ 1;
else
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]-1;
search:=[0,0,B[i][1],B[i][2]];//(a,a)=(0,a)
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]+ 1;
end if;
end for;

//Tabulate the ranks in different fields
Result:=AssociativeArray();
Result[0]:=[n,m];
Result[1]:=number-Rank(a);
primecount:=1;
for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
b:=ChangeRing(a,FiniteField(p)); 
Result[primecount]:=number-Rank(b);
end if;
end for;
result:=Matrix(RationalField(), 1, primecount, [<1,i, Result[i]>: i in [1..primecount]]);
    return result;

end function;

function Table(q,n);
abcount:=0;
ab:=AssociativeArray();
for x in [1..n] do
for y in [x..n] do
if (x eq Gcd(x,y)) or ((x ne Gcd(x,y)) and (y)*(x div Gcd(x,y)) gt n)then 
abcount:=abcount+1;
ab[abcount]:=[x, y];
end if;
end for;
end for;

row:=abcount+1;
primecount:=0;
Prime:=AssociativeArray();

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
Prime[primecount]:=p;
end if;
end for;

column:=primecount+3;
table:=ZeroMatrix(RationalField(), row, column);
primecount:=0;

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
table[1,3+primecount]:=p;
end if;
end for;


for i in [1..abcount] do
finfld:=FindRankFq(q,ab[i][1],ab[i][2]);
table[i+1,1]:=ab[i][1];
table[i+1,2]:=ab[i][2];

for entry in [3..column] do
table[i+1, entry]:=finfld[1,entry-2];
end for; 

end for;

print "[B2 for Zn+Zm]";
print "[ m, n,Q-rank,Fq-rank...]";
return table;
end function;

function Table2(q,n);
ab:=AssociativeArray();
row:=n+1;
primecount:=0;
Prime:=AssociativeArray();

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
Prime[primecount]:=p;
end if;
end for;

column:=primecount+3;
table:=ZeroMatrix(RationalField(), row, column);
primecount:=0;

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
table[1,3+primecount]:=p;
end if;
end for;


for i in [1..n] do
finfld:=FindRankFq(q,i,i);
table[i+1,1]:=i;
table[i+1,2]:=i;
for entry in [3..column] do
table[i+1, entry]:=finfld[1,entry-2];
end for; 
end for;

print "[B2 for Groups of rank 2]";
print "[ n, n,Q-rank,Fq-rank...]";

return table;
end function;

function TableCyclic(q,n);
ab:=AssociativeArray();
row:=n+1;
primecount:=0;
Prime:=AssociativeArray();

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
Prime[primecount]:=p;
end if;
end for;

column:=primecount+3;
table:=ZeroMatrix(RationalField(), row, column);
primecount:=0;

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
table[1,3+primecount]:=p;
end if;
end for;


for i in [1..n] do
finfld:=FindRankFq(q,1,i);
table[i+1,1]:=1;
table[i+1,2]:=i;
for entry in [3..column] do
table[i+1, entry]:=finfld[1,entry-2];
end for; 
end for;

print "[B2 for Cyclic Groups]";
print "[ 1, n,Q-rank,Fq-rank...]";

return table;
end function;


/*
Table(5,10); TableCyclic(5,15); Table2(5,15);
*/