// This MAGMA script uses LattE to count the number of
// eta quotients in a given space of modular forms.

printf "Enter weight k:";
readi k;
printf "Enter level N:";
readi N;

// One constrait is sum r_delta = 2k

// Congruence constraints are
// sum delta r_delta = 0 (mod 24)
// sum (N/delta) r_delta = 0 (mod 24)
// Coefficient of log(p) in sum
// sum r_delta log(delta) is even.

F := Factorization(N);
divlist := Divisors(N);
G := FreeAbelianGroup(#divlist);
gens := [ G.i : i in [1..#divlist]];
grplist := [24,24];
for m in [1..#F] do
  Append(~grplist,2);
end for;
G2 := AbelianGroup(grplist);
homlist := [];
for n in [1..#divlist] do
  d := divlist[n];
  elt := d*G2.1+Floor(N/d)*G2.2;
  for m in [1..#F] do
    elt := elt + Valuation(d,F[m][1])*G2.(m+2);
  end for;
  Append(~homlist,elt);
end for;
phi := hom<G->G2|homlist>;
K := Kernel(phi);

basismatrix := ZeroMatrix(Rationals(),#divlist,#divlist);
for m in [1..#divlist] do
  vec := Eltseq(G!K.m);
  for n in [1..#divlist] do
    basismatrix[m][n] := vec[n];
  end for;
end for;

// Inequalities are for each d | N
// (N/24) sum_(delta | N) (gcd(d,delta)^2 r_delta)/(gcd(d,N/d) d delta) >= 0
// for all divisors d of N.

inequmatrix := ZeroMatrix(Rationals(),#divlist,#divlist);
for m in [1..#divlist] do
  d := divlist[m];
  for n in [1..#divlist] do
    delta := divlist[n];
    entry := (N/24)*(GCD(d,delta)^2)/(GCD(d,Floor(N/d))*d*delta);
    inequmatrix[m][n] := entry;
  end for;
end for;
prodmatrix := basismatrix*Transpose(inequmatrix);
// The matrix prodmatrix has rows that give
// the orders of vanishing of the basis elements at the cusps 1/d
// as d varies over divisors of N

// Now make it integral
denom := LCM({ Denominator(prodmatrix[m][n]) : m in [1..#divlist], n in [1..#divlist]});
prodmatrix := prodmatrix*denom;
lattematrix := Transpose(prodmatrix);
// Make a vector of weights
vec := ZeroMatrix(Rationals(),#divlist,1);
for m in [1..#divlist] do
  vec[m][1] := 1;
end for;
weightvec := basismatrix*vec;

// Write latte file
filename := Sprintf("wt%olev%o.latte",k,N);
System("rm " cat filename);
PrintFile(filename,Sprintf("%o %o",#divlist+1,#divlist+1));
for m in [1..#divlist] do
  str := "0";
  for n in [1..#divlist] do
    str := str cat Sprintf(" %o",lattematrix[m][n]);
  end for;
  PrintFile(filename,str);
end for;
str := Sprintf("%o",2*k);
for n in [1..#divlist] do
  str := str cat Sprintf(" %o",-weightvec[n][1]);
end for;
PrintFile(filename,str);
PrintFile(filename,Sprintf("linearity 1 %o",#divlist+1));

// Call latte, count eta quotients

syscall := Sprintf("./count %o 1> outfile1 2> outfile2\n",filename);
printf "Calling LattE. Output being written to outfile2...";
System(syscall);
printf "Done!\n";
numetaquos := StringToInteger(Read("numOfLatticePoints"));
printf "There are %o eta quotients in M_%o(Gamma_0(%o)).\n",numetaquos,k,N;
printf "Dimension of M_%o(Gamma_0(%o)) = %o.\n",k,N,Dimension(ModularForms(N,k));