// This MAGMA script attempts to prove that
// all modular forms of a given weight and level
// can be expressed via eta quotients.
// It does a "random walk" around
// the eta quotients in the space, and then checks to
// see if every modular function holomorphic except at infinity
// is a linear combination of eta quotients.
// This version attempts to be recursive, by dealing with divisors first.
// Also, when it finds an eta quotient, it computes its images under
// the Atkin-Lehner involutions.

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

datalist := [];

PrintFile("etamake9.out","L := [");

for bigN in [1..1000] do
if #EllipticPoints(Gamma0(bigN)) eq 0 then

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;
  if k*Index(Gamma0(N)) mod 12 ne 0 then
    done := true;
  end if;

  // Do the initial seeding by manually enumerating
  // the eta quotients in the space using ShortVectors.
  if (foundetaquos eq 0) and (done eq false) then
    done := true;
    printf "Finding all 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;
    ind := N*(&*[ 1 + 1/p : p in PrimeDivisors(N)]);
    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;
            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;

spanned, etalist, donelevs, bigetalist, etaquomat, foundetaquos := findetas(2,bigN,donelevs,bigetalist);
etark := NumberOfRows(etaquomat);

generated := false;

if (spanned eq true) then
  // Do weight 4
  wt4etalist := [ etalist[i] + etalist[j] : j in [i..#etalist], i in [1..#etalist]];
  k := 4;
  prec := Floor((k/12)*Index(Gamma0(bigN))+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));
  dim := Dimension(ModularForms(bigN,4));
  printf "Working on weight 4 eta quotients. We have dim M_4(Gamma_0(%o)) = %o.\n",bigN,dim;
  printf "Building q-expansions.";
  etaqlist := [];
  //printf "prec = %o.\n",prec;
  divlist := Divisors(bigN);
  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";
  wt4etaquomat := ZeroMatrix(Rationals(),0,prec+1);
  ind := 1;
  checktime := Min(#wt4etalist+1,dim);
  curetalist := [];
  done := false;
  while done eq false do
    // Compute Fourier expansion of next eta quotient
    rand := Random(#wt4etalist-1)+1;
    cureta := wt4etalist[rand];
    Remove(~wt4etalist,rand);
    Append(~curetalist,cureta);
    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]]];       
    wt4etaquomat := VerticalJoin(wt4etaquomat,vec1);
    ind := ind + 1;
    if (ind ge checktime) then
      printf "So far %o weight 4 eta quotients have been computed.\n",ind;
      newwt4etaquomat, newcuretalist := matreduce(wt4etaquomat,curetalist);
      wt4etaquomat := newwt4etaquomat;
      curetalist := newcuretalist;
      rk := #curetalist;
      printf "The rank is %o.\n",rk;
      if (rk eq dim) then
        done := true;
      else
        checktime := Min(#wt4etalist+ind,ind+500);
        if (#wt4etalist eq 0) then
          done := true;
        end if;
      end if;
    end if;
  end while;  
  printf "The space of weight 4 eta quotients obtained by multiplying weight 2
ones is %o.\n",rk;
  // Now do weight 6.  
  if (rk eq dim) then
    printf "The multiplication map M_2 x M_2 -> M_4 is surjective.\n";
    printf "Now testing weight 6.\n";
    wt6etalist := [ etalist[i] + curetalist[j] : i in [1..#etalist], j in [1..#curetalist]];
    k := 6;
    prec := Floor((k/12)*Index(Gamma0(bigN))+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));
    dim := Dimension(ModularForms(bigN,6));
    printf "We have dim M_6(Gamma_0(%o)) = %o.\n",bigN,dim;
    printf "Building q-expansions.";
    etaqlist := [];
    //printf "prec = %o.\n",prec;
    divlist := Divisors(bigN);
    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";
    wt6etaquomat := ZeroMatrix(Rationals(),0,prec+1);
    ind := 1;
    checktime := Min(#wt6etalist+1,dim);
    curetalist2 := [];
    done := false;
    while done eq false do
      // Compute Fourier expansion of next eta quotient
      rand := Random(#wt6etalist-1)+1;
      cureta := wt6etalist[rand];
      Remove(~wt6etalist,rand);
      Append(~curetalist2,cureta);
      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]]];       
      wt6etaquomat := VerticalJoin(wt6etaquomat,vec1);
      ind := ind + 1;
      if (ind ge checktime) then
        printf "So far %o weight 6 eta quotients have been computed.\n",ind;
        newwt6etaquomat, newcuretalist2 := matreduce(wt6etaquomat,curetalist2);
        wt6etaquomat := newwt6etaquomat;
        curetalist2 := newcuretalist2;
        rk := #curetalist2;
        printf "The rank is %o.\n",rk;
        if (rk eq dim) then
          done := true;
        else
          checktime := Min(#wt6etalist+ind,ind+800);
          if (#wt6etalist eq 0) then
            done := true;
          end if;
        end if;
      end if;
    end while;  
    printf "The space of weight 6 eta quotients obtained by multiplying weight 2
ones is %o.\n",rk;
    if (rk eq dim) then
      printf "The multiplication map M_2 x M_4 -> M_6 is surjective.\n";
      printf "Every even weight modular form for Gamma_0(%o) is a linear combination of products of weight 2 eta quotients.\n",bigN;
      printf "A basis for these eta quotients is given by the following:\n";
      printf "%o",etalist;
      generated := true;
    end if;
  end if;
end if;

if (spanned eq false) then
  k := 2;
  N := bigN;
  // Use Lattice methods to compute dimension of eta 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;
  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;
  foundetaquos := count;
  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);
  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;
  rk := Rank(etaquomat);
  printf "We found %o eta quotients.\n",#etalist;
  printf "The dimension of the space generated by eta quotients is %o.\n",rk;
  etark := rk;
  if (etark eq Dimension(ModularForms(bigN))) then
    spanned := true;
  end if;
end if;

PrintFile("etamake9.out",Sprintf("%o,",<bigN,generated,spanned,Dimension(ModularForms(bigN)),etark,foundetaquos>));

end if;
end for;