/*

1. The command FindM2(n,m) returns the group M2(Zn+Zm).

2. The command TorsionSubgroup(FindM2(n,m)) will return the torsion subgroup of M2(Zn+Zm).

3. The comman TorsionFreeRank(FindM2(n,m)) will return the Q-rank of M2(Zn+Zm)

*/

//Create the function checking whether (i1,j1) and (i2,j2) generates Zn+Zm (n<=m)
CheckGenerator1s:=function (n,m,i1,j1,i2,j2);

if n eq 1 then//Simplify the case when Z/n+Z/m=Z/mn
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;

 

//Create the function returning M2
FindM2:=function (n,m);
number:=0;
row:=0;
relation:=0;
A:= AssociativeArray();//where (i1,j1) is stored
B:= AssociativeArray();//where (i2,j2) is stored
Index:=AssociativeArray();//index the generator1 ((i1,j1),(i2,j2))
//Make sure Zn+Zm is reduced to the simpliefied form
if Gcd(n,m) eq 1 then
m:=m*n;
n:=1;
else
prod:=n*m;
n:=Gcd(n,m);
m:=prod div n;
end if;

//Find the generating symbols
for i1 in [0..n-1] do
for j1 in [0..m-1] do
for i2 in [0..n-1] do
for j2 in [0..m-1] do

if (CheckGenerator1s(n,m,i1,j1,i2,j2)) then
number:=number+1;
A[number] := [i1, j1];//Store(i1,j1)
B[number] := [i2, j2];//Store(i2,j2)
Index[i1*(10^9)+j1*(10^6)+i2*(10^3)+j2]:=number;//Store the index
end if;
end for;
end for;
end for;
end for;

//Create M2(Zn+Zm)
M:=FreeAbelianGroup(number);//Create the free groups on symbols
relation:={};

for i in [1..number] do

if A[i] eq B[i] then
//On symbols of the form (a,a)
search:=A[i][1]*(10^9)+A[i][2]*(10^6);
target:=Index[search];
relation:=Include(relation, 2*M.target-M.i);//Relation (M)
else
//On symbols of the form (a,b), a neq b
if A[i] le B[i] then//(a,b) and (b,a) will procude the same relation
search:=B[i][1]*(10^9)+B[i][2]*(10^6)+A[i][1]*(10^3)+A[i][2];
target:=Index[search];
relation:=Include(relation, M.target-M.i);//Symmetric Relation

i1:=(A[i][1]-B[i][1]) mod n;
j1:=(A[i][2]-B[i][2]) mod m;
i2:=(B[i][1]-A[i][1]) mod n;
j2:=(B[i][2]-A[i][2]) mod m;
Generator1:=[A[i], [i2, j2]];
Generator2:=[[i1, j1], B[i]];
search:=Generator2[1][1]*(10^9)+Generator2[1][2]*(10^6)+Generator2[2][1]*(10^3)+Generator2[2][2];
target1:=Index[search];
search:=Generator1[1][1]*(10^9)+Generator1[1][2]*(10^6)+Generator1[2][1]*(10^3)+Generator1[2][2];
target2:=Index[search];
relation:=Include(relation, M.target1+M.target2-M.i);//Relation(M)
end if;
end if;

end for;

M2:= quo<M|relation>;
torsion:=TorsionSubgroup(M2);
return torsion ;
end function;

/* Sample running:

print "This is a list of torsions of M_2";
for i in [1..20] do
print "";
print "";
print "";
print "Torsion of M_2(Z/" cat IntegerToString(i) cat ") is:";
FindM2(1,i);
print "";
print "";
print "";
end for;

for i in [1..5] do
print "";
print "";
print "";
print "Torsion of M_2(Z/" cat IntegerToString(i) cat "+Z/" cat IntegerToString(i) cat ") is:";
FindM2(i,i);
print "";
print "";
print "";
end for;


*/