// This MAGMA script attempts to prove that
// E_2 = R_2. To do this, it suffices 
// to check that M_r+2(Gamma_0(N)) is generated
// by eta quotients, where r is the Sturm bound weight.

// This script takes a level N and attempts to
// show that every holomorphic modular form of level N is
// a linear combination of eta quotients.

printf "Enter level N:";
readi bigN;

Ndivlist := Divisors(bigN);
donelevs := [ 0 : i in [1..#Ndivlist]];
bigetalist := <>;
for i in [1..#Ndivlist] do
  Append(~bigetalist,[ RSpace(Rationals(),#Divisors(Ndivlist[i])) | ]);
end for;

// V(d) operator
// It takes a power series f, and returns f|V(d)

function V(f,d)
  qq := Parent(f).1;
  start := Valuation(f);
  en := AbsolutePrecision(f)-1;
  ret := 0;  
  for n in [start..en] do
    ret := ret + Coefficient(f,n)*qq^(n*d);
  end for;
  ret := ret + BigO(qq^(AbsolutePrecision(f)*d));
  return ret;
end function;

// 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;

function findetas(k,N,donelevs,bigetalist)

  printf "Searching for eta quotients of weight k = %o and level N = %o.\n",k,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.
  // Also, find weight zero forms

  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

  // Do a search to find the eta quotients
  // Make a basis of the weight zero eta quotients in the space.

  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;
  Append(~grplist,0);
  G2 := AbelianGroup(grplist);
  homlist := [];
  for n in [1..#divlist] do
    d := divlist[n];
    elt := d*G2.1+Floor(N/d)*G2.2+G2.(#grplist);
    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);
  basismatrix2 := ZeroMatrix(Rationals(),#divlist-1,#divlist);
  for m in [1..#divlist-1] do
    vec := Eltseq(G!K.m);
    for n in [1..#divlist] do
      basismatrix2[m][n] := vec[n];
    end for;
  end for;
  // LLL the basis matrix
  basismatrix2 := RMatrixSpace(Rationals(),#divlist-1,#divlist)!BasisMatrix(Lattice(basismatrix2));

  // Start with one eta quotient
  VV := RSpace(Rationals(),#divlist);
  basismatrix := RMatrixSpace(Integers(),#divlist,#divlist)!basismatrix;

  etaquoset := {};
  queue := [];
  checked := {};
  queueset := {};
  goodetaquolist := [];

  prec := Floor((k/12)*Index(Gamma0(N))+1);
  R<q> := PowerSeriesRing(Rationals() : Precision := 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 "etaqlist is %o.\n",etaqlist;

  etaquomat := ZeroMatrix(Rationals(),0,prec+1);
  goodetaquolist := [];
  V2 := RMatrixSpace(Rationals(),1,prec+1);
  dim := Dimension(ModularForms(N,k));
  printf "The dimension is %o.\n",dim;

  foundetaquos := 0;

// Do recursive calls to start with
  atkinlist := [ Q : Q in Divisors(N) | GCD(Q,N div Q) eq 1 ];
  //printf "Atkin list is %o.\n",atkinlist;
  for d in Divisors(N) do
    if (donelevs[Index(Ndivlist,d)] eq 0) and (d lt N) then
      foundall, oldetaquolist, donelevs, bigetalist := findetas(k,d,donelevs,bigetalist);
    else
      oldetaquolist := bigetalist[Index(Ndivlist,d)];
    end if;
    olddivlist := Divisors(d);
    for e in Divisors(Floor(N/d)) do
      // Add old eta quotients and those hit with V(e) to etaquomat
      // and goodetaquolist
      for i in [1..#oldetaquolist] do
        newetaquo := VV!0;
        for mm in [1..#olddivlist] do
          newetaquo[Index(divlist,e*olddivlist[mm])] := oldetaquolist[i][mm];
        end for;
        for Q in atkinlist do
          newetaquo2 := VV!0;
          for m in [1..#divlist] do
            dd := GCD(N div Q, divlist[m]);
            newetaquo2[Index(divlist,(Q*dd^2) div divlist[m])] := newetaquo[m];
          end for;
          if not (newetaquo2 in etaquoset) then
            Include(~etaquoset,newetaquo2);
            foundetaquos := foundetaquos + 1;
            Append(~queue,<newetaquo2,0>);
            Include(~checked,newetaquo2);
            Include(~queueset,newetaquo2);
            val1 := (1/24)*&+[ divlist[i]*newetaquo2[i] : i in [1..#divlist]];
            val1 := Integers()!val1;
            //printf "The e = %o, i = %o Q = %o eta quotient is %o.\n",e,i,Q,newetaquo2;
            newetaq := q^(val1)*&*[ etaqlist[i]^(Integers()!newetaquo2[i]) : i in [1..#divlist]];
            vec1 := V2![[ Coefficient(newetaq,i) : i in [0..prec]]];       
            Append(~goodetaquolist,newetaquo2);
            etaquomat := VerticalJoin(etaquomat,vec1);
          end if;
        end for;
      end for;
    end for;
  end for;
  done := false;
  if foundetaquos gt 0 then
    etaquomat, goodetaquolist := matreduce(etaquomat,goodetaquolist);
    rk := NumberOfRows(etaquomat);
    rankcheck := foundetaquos + Max(250,dim-rk);
    if (rk eq dim) then
      done := true;
    end if;
  else
    rankcheck := dim;
  end if;
  e2 := 0;
  chi4 := KroneckerCharacter(-4);
  chi3 := KroneckerCharacter(-3);
  if (N gt 1) and (N mod 4 ne 0) then
    e2 := &+[ 1 + chi4(p) : p in PrimeDivisors(N)];
  end if;
  if (N gt 1) and (N mod 9 ne 0) then
    e3 := &+[ 1 + chi3(p) : p in PrimeDivisors(N)];
  end if;
  if (N eq 1) then
    e2 := 1;
    e3 := 1;
  end if;
  
  if (k mod 6 ne 0) and (e3 gt 0) then
    done := true;
  end if;
  if (k mod 4 ne 0) and (e2 gt 0) then
    done := true;
  end if;
  // If done is false, the space contains at least one eta quotient.

  // Do the initial seeding by manually enumerating
  // the eta quotients in the space using ShortVectors.
  if (done eq false) and (foundetaquos eq 0) then 
    done := true;
    foundetaquos := 0;
    printf "Finding *some* eta quotients by using short vectors.\n";
    prodmatrix3 := prodmatrix;
    for m in [1..#divlist] do
      MultiplyColumn(~prodmatrix3,EulerPhi(GCD(divlist[m],N div divlist[m])),m);
    end for;
    L := Lattice(prodmatrix3);
    //printf "prodmatrix = %o.\n",prodmatrix3;
    if (N gt 1) then
      ind := N*(&*[ 1 + 1/p : p in PrimeDivisors(N)]);
    else
      ind := 1;
    end if;
    bound := Ceiling(((k/12)*ind)^2);
    A := ShortVectorsProcess(L,bound/#divlist,bound);
    matinv := Transpose(inequmatrix)^(-1);
    donedone := false;
    count := 0;
    //printf "bound = %o.\n",bound;
    while donedone eq false do
      v, nm := NextVector(A);
      //printf "v = %o.\n",v;
      if (nm eq -1) then
        donedone := 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)[1];
            //printf "Exponent vector is %o.\n",expovec;
            foundetaquos := foundetaquos + 1;
            if (foundetaquos mod 100) eq 0 then
              printf "Found %o eta quotients so far.\n",foundetaquos;
              if (foundetaquos ge 1000) and (k ge 30) then
                donedone := true;
              end if;
            end if;
            Include(~etaquoset,expovec);
            Append(~queue,<expovec,0>);
            Include(~checked,expovec);
            Include(~queueset,expovec);
            Append(~goodetaquolist,expovec);
            val := (1/24)*&+[ divlist[i]*expovec[i] : i in [1..#divlist]];
            etaq := (q^val)*(&*[ etaqlist[i]^(Integers()!expovec[i]) : i in [1..#divlist]]);
            vec := V2![[ Coefficient(etaq,i) : i in [0..prec]]];
            etaquomat := VerticalJoin(etaquomat,vec);
          end if;
        end if;
      end if;
    end while;
  end if;

  if foundetaquos gt 0 then
    etaquomat, goodetaquolist := matreduce(etaquomat,goodetaquolist);
    rk := NumberOfRows(etaquomat);
    rankcheck := foundetaquos + dim-rk;
    if (rk eq dim) then
      done := true;
    end if;
  else
    rankcheck := dim;
  end if;

  passnum := 0;
  room := 3;
  lastcheck := 0;
  lastcheckcount := foundetaquos;

  if (done eq false) then
    printf "Starting the search for eta quotients in M_%o(Gamma_0(%o)).\n",k,N;
  end if;
  while (done eq false) do
    passnum := passnum + 1;
    if (passnum mod 500 eq 0) then
      printf "Pass #%o - found = %o, queue = %o, last rank = %o, next check = %o.\n",passnum,foundetaquos,#queue,rk,rankcheck;
    end if;
    rand := Random(#queue-1)+1;
    node := queue[rand][1];
    mult := queue[rand][2];
    Remove(~queue,rand);
    Exclude(~queueset,node);
    if #queue ge 50000 then
      room := 2;
    else
      room := 3;
    end if;
    testlist := [ node + VV!basismatrix2[i] : i in [1..#divlist-1]]
    cat [ node - VV!basismatrix2[i] : i in [1..#divlist-1]];
    for m in [1..#testlist] do
      v2 := testlist[m];
      if not (v2 in checked) then
        if not (v2 in etaquoset) then
          ordlist := v2*Transpose(inequmatrix);
          good := true;
          for m2 in [1..#divlist] do
            if ordlist[m2] lt 0 then
              good := false;
            end if; 
          end for;
          if good eq true then
            for Q in atkinlist do
              v3 := VV!0;
              for m in [1..#divlist] do
                dd := GCD(N div Q, divlist[m]);
                v3[Index(divlist,(Q*dd^2) div divlist[m])] := v2[m];
              end for;
              if not (v3 in etaquoset) then
                Include(~checked,v3);
                Include(~etaquoset,v3);
                Append(~queue,<v3,0>);
                Include(~queueset,v3);
                Append(~goodetaquolist,v3);
                foundetaquos := foundetaquos + 1;
                //if (foundetaquos ge 100) then
                //  room := 0;
                //end if;
                // Add to etaquomat
                val := (1/24)*&+[ divlist[i]*v3[i] : i in [1..#divlist]];
                val := Integers()!val;
                etaq := q^(val)*&*[ etaqlist[i]^(Integers()!v3[i]) : i in [1..#divlist]];
                vec := V2![[ Coefficient(etaq,i) : i in [0..prec]]];
                etaquomat := VerticalJoin(etaquomat,vec);
                numrows := NumberOfRows(etaquomat);
              end if;
            end for;
          else
            if mult lt room then
              if (v2 in queueset) then
                ind := Index([ m[1] : m in queue],v2);
                queue[ind][2] := Min(queue[ind][2],mult+1);
              else
                Append(~queue,<v2,mult+1>);
                Include(~queueset,v2);
              end if;
            end if;
          end if;
        end if;
      end if;
    end for;
    if (foundetaquos ge rankcheck) or ((passnum ge (lastcheck + 5000)) and (foundetaquos gt lastcheckcount)) then
      etaquomat, goodetaquolist := matreduce(etaquomat,goodetaquolist);
      rk := NumberOfRows(etaquomat);
      if rk eq dim then 
        done := true;
      else
        rankcheck := foundetaquos + Max(100,dim - rk);
        lastcheck := passnum;
        lastcheckcount := foundetaquos;
      end if;      
    end if;                    
    if (#queue eq 0) then
      done := true;
    end if;  
  end while;
  if foundetaquos gt 0 then
    etaquomat, goodetaquolist := matreduce(etaquomat,goodetaquolist);
  end if;
  rk := NumberOfRows(etaquomat);
  
  printf "We found %o eta quotients!\n",foundetaquos;
  printf "The dimension of M_%o(Gamma_0(%o)) = %o.\n",k,N,dim;
  printf "The dimension of the space spanned by the eta quotients is %o.\n",rk;
  foundall := false;
  if (rk eq dim) then
    foundall := true;
  end if;
  donelevs[Index(Ndivlist,N)] := 1;
  //printf "Parent of goodetaquolist is %o.\n",Parent(goodetaquolist);
  //This last line was modified. Before, it was goodetaquolist  
  bigetalist[Index(Ndivlist,N)] := SetToSequence(etaquoset);
  return foundall, goodetaquolist, donelevs, bigetalist, etaquomat, foundetaquos;
end function;

// Build magic Sturm bound eta quotient of level N
// and smallest possible weight

divlist := Divisors(bigN);
magiceta := [ MoebiusMu(Floor(bigN/d))*d : d in divlist ];
dv := GCD(magiceta);
magicetaold := magiceta;
magiceta := [ Floor(magiceta[i]/dv) : i in [1..#divlist]];
sum1 := &+[ divlist[i]*magiceta[i] : i in [1..#divlist]];
sum2 := &+[ Floor(bigN/divlist[i])*magiceta[i] : i in [1..#divlist]];
prod := &*[ divlist[i]^magiceta[i] : i in [1..#divlist]];
//printf "Initial magiceta = %o.\n",magiceta;
alpha := 1;
alpha := LCM(alpha,Floor(24/GCD(sum1,24)));
alpha := LCM(alpha,Floor(24/GCD(sum2,24)));
if not IsSquare(prod) then
  alpha := LCM(alpha,2);
end if;
wtchk := &+[ magiceta[i] : i in [1..#divlist]];
if (wtchk mod 2 eq 1) then
  alpha := LCM(alpha,4);
end if;
if (wtchk mod 4 eq 2) then
  alpha := LCM(alpha,2);
end if;

/*
if (bigN eq 529) then
  alpha := alpha*2;
end if;
*/
/*
if (bigN eq 121) then
  alpha := alpha*2;
end if;
*/

magicetaold := magiceta;
magiceta := [ alpha*magicetaold[i] : i in [1..#divlist]];
printf "An eta quotient achieving the Sturm bound is %o.\n",magiceta;
k := Floor((1/2)*&+[ magiceta[i] : i in [1..#divlist]]);
printf "The weight is %o.\n",k;

foundall, goodetaquolist, donelevs, bigetalist, etaquomat := findetas(k+2,bigN,donelevs,bigetalist);

dim := Dimension(ModularForms(bigN,k+2));
if NumberOfRows(etaquomat) eq dim then
  printf "Level %o is good!\n",bigN;
else
  printf "Things don't work for level %o.\n",bigN;
end if;