/*
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 M2^-(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 M2^-(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. M2^-(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 M2^-(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 M2Generators(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:=M2Generators(n,m);

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

for i in [1..number] do

//Minus sign relation (a,b)=(-a,b)
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
search:=[[(-A[i][1]) mod n ,(-A[i][2]) mod m], [B[i][1],B[i][2]]];
search:=Sort(search);
search:=[search[1][1],search[1][2],search[2][1],search[2][2]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]+ 1;

//Minus sign relation (a,b)=(a,-b)
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
search:=[[A[i][1] ,A[i][2]], [(-B[i][1]) mod n,(-B[i][2]) mod m]];
search:=Sort(search);
search:=[search[1][1],search[1][2],search[2][1],search[2][2]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]+ 1;


//(M) relation
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;

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 "[M2^- 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 "[M2^- 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 "[M2^- for Cyclic Groups]";
print "[ 1, n,Q-rank,Fq-rank...]";

return table;
end function;



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