// This MAGMA script tries to enumerate all the eta quotients
// of a given weight and level. This should work better for low
// weights then high weights. It searches for short vectors in lattices.

// This function takes a matrix of eta quotient
// q-expansions and a list of eta quotients
// and "reduces" the list and matrix. It returns
// a matrix whose rows are linearly independent
// and a reduced set of eta quotients

function matreduce(etaquomat,goodetaquolist)
  // Tricky rank checking
  rand_prime_list := [];
  goodrows := {};
  while #rand_prime_list lt 3 do
    rp := RandomPrime(25);
    if (Index(rand_prime_list,rp) eq 0) and (rp gt 100000) then
      Append(~rand_prime_list,rp);
    end if;
  end while;
  for rp in rand_prime_list do
    etaquomatfp := ChangeRing(etaquomat,GF(rp));
    ech := EchelonForm(Transpose(etaquomatfp));
    for i in [1..NumberOfRows(ech)] do
      pivot_flag := false;
      for j in [1..NumberOfColumns(ech)] do
        if (ech[i][j] ne 0) and (pivot_flag eq false) then
          pivot_flag := true;
          Include(~goodrows,j);
        end if;
      end for;
    end for;
  end for;
  etaquomat2 := Matrix([ etaquomat[i] : i in goodrows ]);
  goodetaquolist2 := [ goodetaquolist[i] : i in goodrows ];
  return etaquomat2, goodetaquolist2;
end function;

printf "Enter the weight k:";
readi k;
printf "Enter the 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

// Scale prodmatrix so that the columns are weighted by the actual order
// of vanishing
for m in [1..#divlist] do
  MultiplyColumn(~prodmatrix,EulerPhi(GCD(divlist[m],N div divlist[m])),m);
end for;
printf "prodmatrix = %o.\n",prodmatrix;
L := Lattice(prodmatrix);
ind := N*(&*[ 1 + 1/p : p in PrimeDivisors(N)]);
bound := Ceiling(((k/12)*ind)^2);
printf "Bound = %o.\n",bound;
A := ShortVectorsProcess(L,bound/#divlist,bound);
matinv := Transpose(inequmatrix)^(-1);
done := false;
count := 0;
etalist := [];
while done eq false do
  v, nm := NextVector(A);
  //printf "v = %o.\n",v;
  if (nm eq -1) then
    done := true;
  else
    if &and[ v[i] ge 0 : i in [1..#divlist]] then
      wt := 12*(&+[ v[i] : i in [1..#divlist]])/ind;
      if (wt eq k) then
        // We found one!
        newv := ZeroMatrix(Rationals(),1,#divlist);
        for m in [1..#divlist] do
          newv[1][m] := v[m]/EulerPhi(GCD(divlist[m],N div divlist[m]));
        end for;
        expovec := newv*matinv;
        printf "Exponent vector is %o.\n",expovec;
        Append(~etalist,expovec[1]);
        count := count + 1;
      end if;
    end if;
  end if; 
end while;

printf "We found %o eta quotients.\n",#etalist;
prec := Floor((k/12)*Index(Gamma0(N))+1);
R<q> := PowerSeriesRing(Rationals() : Precision := prec+1);
V2 := RMatrixSpace(Rationals(),1,prec+1);
// Get eta quotient q-expansion, sans power of q
indmax := Floor((1+Sqrt(24*prec+1))/2);
etaqexp := R!&+[(-1)^k*q^(Floor(k*(3*k-1)/2))  : k in [-indmax..indmax]];
etaqexp := etaqexp + BigO(q^(prec+1));
printf "Building q-expansions.";
etaqlist := [];
//printf "prec = %o.\n",prec;
for i in [1..#divlist] do
  ent := R!BigO(q^(prec+1));
  for j in [0..Min(prec,Floor((prec+1)/divlist[i]))] do
    ent := ent + Coefficient(etaqexp,j)*q^(divlist[i]*j);
  end for;
  Append(~etaqlist,ent);
end for;
printf "Done!\n";
printf "Dimension of M_%o(Gamma_0(%o)) = %o.\n",k,N,Dimension(ModularForms(N,k));
etaquomat := ZeroMatrix(Rationals(),0,prec+1);
printf "Building q-expansions...\n";
for m in [1..#etalist] do
  cureta := etalist[m];
  val1 := (1/24)*&+[ divlist[i]*cureta[i] : i in [1..#divlist]];
  val1 := Integers()!val1;
  newetaq := q^(val1)*&*[ etaqlist[i]^(Integers()!cureta[i]) : i in [1..#divlist]];
  vec1 := V2![[ Coefficient(newetaq,i) : i in [0..prec]]];
  etaquomat := VerticalJoin(etaquomat,vec1);
end for;
printf "Done.\n";
printf "Rank checking.\n";
etaquomatnew, newetalist := matreduce(etaquomat,etalist);
printf "Done.\n";
rk := Rank(etaquomatnew);
printf "The dimension of the space generated by eta quotients is %o.\n",rk;
printf "Codimension = %o.\n",Dimension(ModularForms(N,k))-rk;