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

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 M3(Zm+Zk+Zr) where m,k,r<=n.

3.The command Table3(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 3 i.e. M3(Zm+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 M3(Zm) where m<=n.

*/

//Create the function checking generation condition
function CheckGenerators(n,m,r,i1,j1,k1,i2,j2,k2,i3,j3,k3);

if (n eq 1) and (m eq 1) then//When G is cyclic 

if Gcd(Gcd(Gcd(k1,k2), k3), r) eq 1 then
return true;
else 
return false;
end if; 

elif (n eq 1) and (m ne 1) then//When G has two generators

Z:=RModule(IntegerRing(),2);
L:=sub<Z|[m,0],[0,r]>;
    Znm:=quo<Z|L>;
      gen:=sub<Znm|[j1,k1],[j2,k2], [j3,k3]>;

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

else//When G has three generators

Z:=RModule(IntegerRing(),3);
L:=sub<Z|[n,0,0],[0,m,0],[0,0,r]>;
  Znm:=quo<Z|L>;
    gen:=sub<Znm|[i1,j1,k1],[i2,j2,k2],[i3,j3,k3]>;
      
      if gen eq Znm then
      return true;
      else
      return false;
      end if;
      
      end if;

end function;

// Create the function sending ZN+ZM+ZK to standarad form compatible with quotient module in Magma
function StandardForm(n,m,r);
S:=Sort([n,m,r]);
n:=S[1];
m:=S[2];
r:=S[3];
prod:=n*m;
n:=Gcd(m,n);
m:=prod div n;
prod:=m*r;
m:=Gcd(m,r);
r:=prod div m;
return Sort([n,m,r]);
end function;

function FindRankFq(q,n,m,r);
//turn Zn+Zm+Zk to standard form 
Order:=StandardForm(n,m,r);
n:=Order[1];
m:=Order[2];
r:=Order[3];

A:=AssociativeArray();
B:=AssociativeArray();
C:=AssociativeArray();
M3NGens:=0; //number of generators 
Index:=AssociativeArray(RModule(IntegerRing(),9)); //record the order of generators

for i1 in [0..n-1] do
for j1 in [0..m-1] do
for k1 in [0..r-1] do
for i2 in [0..n-1] do
for j2 in [0..m-1] do
for k2 in [0..r-1] do
for i3 in [0..n-1] do
for j3 in [0..m-1] do
for k3 in [0..r-1] do
if CheckGenerators(n,m,r,i1,j1,k1,i2,j2,k2,i3,j3,k3) then
M3NGens:=M3NGens+1;
A[M3NGens]:=[i1,j1,k1];
B[M3NGens]:=[i2,j2,k2];
C[M3NGens]:=[i3,j3,k3];
Index[[i1,j1,k1,i2,j2,k2,i3,j3,k3]]:=M3NGens;
end if;
end for;
end for;
end for;
end for;
end for;
end for;
end for;
end for;
end for;

//Create the relation matrix
row:=((M3NGens*3) div 2)+2+M3NGens;
a:=SparseMatrix(RationalField(), row, M3NGens);
currentrow:=0;

for i in [1..M3NGens] do

//(M) relation
if A[i] le 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;
k1:=(A[i][3]-B[i][3]) mod r;
search:=[i1,j1,k1,B[i][1],B[i][2],B[i][3],C[i][1],C[i][2],C[i][3]];
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;
k2:=(B[i][3]-A[i][3]) mod r;
search:=[A[i][1],A[i][2],A[i][3],i2,j2,k2,C[i][1],C[i][2],C[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]+ 1;

end if;

//Symmetry relation
if (A[i] lt B[i]) and (B[i] lt C[i])  then//(a,b,c) with 6 symmetry relations

currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[A[i][1],A[i][2],A[i][3],C[i][1],C[i][2],C[i][3],B[i][1],B[i][2],B[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;

currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[B[i][1],B[i][2],B[i][3],A[i][1],A[i][2],A[i][3],C[i][1],C[i][2],C[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;

currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[B[i][1],B[i][2],B[i][3],C[i][1],C[i][2],C[i][3],A[i][1],A[i][2],A[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;


currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[C[i][1],C[i][2],C[i][3],A[i][1],A[i][2],A[i][3],B[i][1],B[i][2],B[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;

currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[C[i][1],C[i][2],C[i][3],B[i][1],B[i][2],B[i][3],A[i][1],A[i][2],A[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;

elif (A[i] eq B[i]) and (B[i] ne C[i]) then//(a,a,c)=(a,c,a)=(c,a,a)
currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[C[i][1],C[i][2],C[i][3],B[i][1],B[i][2],B[i][3],A[i][1],A[i][2],A[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;

currentrow:=currentrow+1;
a[currentrow, i]:=a[currentrow, i]+1;
search:=[B[i][1],B[i][2],B[i][3],C[i][1],C[i][2],C[i][3],A[i][1],A[i][2],A[i][3]];
target:=Index[search];
a[currentrow, target]:=a[currentrow, target]-1;
end if;
end for;

//Tabulate the result
Result:=AssociativeArray();
Result[0]:=[n,m];
Result[1]:=M3NGens-Rank(a);
primecount:=1;

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
b:=ChangeRing(a,FiniteField(p)); 
Result[primecount]:=M3NGens-Rank(b);
end if;
end for;

result:=Matrix(RationalField(), 1, primecount, [<1,i, Result[i]>: i in [1..primecount]]);
    
    return result;

end function;

procedure Table(q,n);

abcount:=0;
ab:=AssociativeArray();

for x in [1..n] do
for y in [x..n] do
for z in [y..n] do
if (StandardForm(x,y,z) eq [x,y,z]) or (StandardForm(x,y,z)[1] gt n) or (StandardForm(x,y,z)[2] gt n)or (StandardForm(x,y,z)[3] gt n) then 
abcount:=abcount+1;
ab[abcount]:=[x, y, z];
end if;
end for;
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+4;
table:=ZeroMatrix(RationalField(), row, column);
primecount:=0;

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

print "[M3 for Zn+Zm]";
print "[ m, n, r, Q-rank,Fq-rank...]";
print [table[1][j]:j in [1..primecount+4]];

for i in [1..abcount] do

finfld:=FindRankFq(q,ab[i][1],ab[i][2],ab[i][3]);
table[i+1,1]:=ab[i][1];
table[i+1,2]:=ab[i][2];
table[i+1,3]:=ab[i][3];

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

print [table[i+1][j]:j in [1..primecount+4]];

end for;
end procedure;

procedure Table3(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+4;
table:=ZeroMatrix(RationalField(), row, column);
primecount:=0;

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


print "[M3 for Groups of rank 3]";
print "[ ZN+ZN+ZN,  Q-rank,Fq-rank...]";
print [table[1][j]:j in [1..primecount+4]];

for i in [1..n] do

finfld:=FindRankFq(q,i,i,i);
table[i+1,1]:=i;
table[i+1,2]:=i;
table[i+1,3]:=i;

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

print [table[i+1][j]:j in [1..primecount+4]];

end for;

end procedure;

procedure 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+4;
table:=ZeroMatrix(RationalField(), row, column);
primecount:=0;

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


print "[M3 for cyclic Groups]";
print "[ ZN,  Q-rank,Fq-rank...]";
print [table[1][j]:j in [1..primecount+4]];

for i in [1..n] do

finfld:=FindRankFq(q,1,1,i);
table[i+1,1]:=1;
table[i+1,2]:=1;
table[i+1,3]:=i;

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

print [table[i+1][j]:j in [1..primecount+4]];

end for;

end procedure;

/*
Table3(5,3); TableCyclic(5,10); Table(5,3); 
*/