/*
The function M2(n,m) returns the Q-rank of M2(Cn+Cm). For the case of M2(Cn), please input M2(1,n);
The function M2m(n,m) returns the Q-rank of M2-(Cn+Cm). For the case of M2-(Cn), please input M2m(1,n);
*/


//M1+
function a(n);
if (n eq 1) or (n eq 2) then
return 1;
else
return EulerPhi(n) div 2;
end if;
end function;

//M1-
function b(n);
if n eq 2 then
return 0;
elif n eq 1 then
return 0;
else
return EulerPhi(n) div 2;
end if;
end function;

//M2-(Cn+Cnm) n ne 2
function c1(n,m);
dim:=1;
for i in [2..n*m] do
if IsPrime(i) and (n*m mod i eq 0)then
dim:=dim*(1-i^(-2));
end if;
end for;
if n eq 2 then
s:=1;
else
s:=EulerPhi(n) div 2;
end if;
return s*(1+m^2*n^3*dim/12);
end function;


//M2-(Cn+Cnm) n=2
function c2(m)
dim:=1;
for i in [2..2*m] do
if IsPrime(i) and (2*m mod i eq 0) then
dim:=dim*(1+1/i);
end if;
end for;
c:=(2*EulerPhi(m)+EulerPhi(2*m)) div 2;
return 1+m*EulerPhi(2*m)/6*dim-c;
end function;


//Ker(M2-->M2-(Cn+Cnm)
function f(n,m);
	  A:=FreeAbelianGroup(2);
	  A:=quo<A|n*A.1, (n*m)*A.2>;
	  record:=[];
	  dim:=0;
	  for x in A do
	  B:=sub<A|x>;
	  C:=quo<A|B>;
	  if Ngens(C) eq 1 then

	  if not (B in record) then
	  dim:=dim+a(Order(x))*b((n*n*m) div Order(x));
	  record:=record cat [B];
	  end if;

	  end if;
	  end for;

return dim;
end function;

//ker(M2-->M2-(Cn))
function g(n);
	  A:=FreeAbelianGroup(1);
	  A:=quo<A|n*A.1>;
	  record:=[];
	  dim:=0;
	  for x in A do
	  B:=sub<A|x>;
	  C:=quo<A|B>;
	  if Ngens(C) eq 1 then

	  if not (B in record) then
	  dim:=dim+a(Order(x))*b(n div Order(x));
	  record:=record cat [B];
	  end if;

	  end if;
	  end for;

return dim;
end function;

//M2-(Cn)
function h(n);
	  mu:=1;
	  for i in [2..n] do
	  if IsPrime(i) and (n mod i eq 0) then
	  mu:=mu*(1+1/i);
	  end if;
	  end for;
	  mu:=mu*n*EulerPhi(n)/2;

	  if n mod 2 eq 0 then
	  C2:=EulerPhi(n)+EulerPhi(n div 2);
	  else
	  C2:=EulerPhi(n);
	  end if;

	  return 1+mu/12-C2/2;
end function;


//M2(Cn)
function kn(n);
	  return h(n)+g(n);
end function;


//M2(Cn+Cnm)
function k(n,m);
	  if n eq 2 then
	  return c2(m)+f(n,m);
	  else
	  return c1(n,m)+f(n,m);
	  end if;

end function;


function M2(n,m);

	  if Gcd(n,m) eq 1 then
	  n:=n*m;
	  if n eq 2 then
	  return 0;
	  else
	  return kn(n);
	  end if;
	  else
	  prod:=n*m;
	  n:=Gcd(n,m);
	  m:=prod div n;
	  if n eq 2 and m eq 2 then
	  return 0;
	  else
	  return k(n, m div n);
	  end if;
	  end if;
end function;


function M2m(n,m);

	  if Gcd(n,m) eq 1 then
	  n:=n*m;
	  if n eq 2 then
	  return 0;
	  else
	  return h(n);
	  end if;
	  else
	  prod:=n*m;
	  n:=Gcd(n,m);
	  m:=prod div n;
	  if n eq 2 then
	  if m eq 2 then
	  return 0;
	  else 
	  return c2(m div 2);
	  end if;
	  else
	  return c1(n, m div n);
	  end if;
	  end if;
end function;