/***********************************************************

   Computations related to Bas's paper cokerq.ps 
   William Stein, April 9, 2002.

 ***********************************************************/

function MakeLattice(B)
   Mat := RMatrixSpace(Rationals(),#B, #Eltseq(B[1]));
   X   := Mat!B;
   if Dimension(RowSpace(X)) eq 0 then
      return RowSpace(X);
   end if;
   return Lattice(X);
end function;

intrinsic S1F2IsGood(N::RngIntElt) -> BoolElt
{True if and only if S1F2(N) is defined.}
   if IsEven(N) then
      return false;
   end if;
   is_good := false;
   F := Factorization(N);
   if N mod 9 eq 0 then
      if #[p : p in F | p[1] mod 4 eq 3] gt 0 then
         is_good := true;
      end if;               
   end if;
   if not is_good and #[p : p in F | p[1] mod 3 eq 2] gt 0 then
      is_good := true;
   end if;            
   return is_good;   
end intrinsic;

intrinsic DimS1F2(N::RngIntElt : Verbose := true) -> RngIntElt, GrpAb
{The F_2 dimension of S_1(N, F_2) and Edixhoven's module M = S_2(2N,Z)/(S_2(2N,Z) meet W_2(S_2(2N,Z))).}
   require GCD(N,2) eq 1 : "Argument 1 must be odd.";
   require S1F2IsGood(N) : "Argument 1 must be divisible by a prime = 3 mod 4 and by 9, or by a prime = 2 mod 3.";   

   // Compute a basis of q-expansions for S_2(Gamma_0(2*N))^new.
   if Verbose then
      printf "Compute a basis of q-expansions for S_2(Gamma_0(2*%o))^new. ", N;
      t := Cputime();
   end if;

   S2N := CuspidalSubspace(ModularForms(2*N,2));
   if Dimension(S2N) eq 0 then
      return 0;
   end if;
   bnd := PrecisionBound(S2N) + 1;
   B2N := [PowerSeries(f,bnd) : f in Basis(NewSubspace(S2N,2))];

   q := PowerSeriesRing(IntegerRing()).1;
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;
 

   // Compute a basis of q-expansions for S_2(Gamma_0(N)).
   SN  := CuspidalSubspace(ModularForms(N,2));
   BN  := [PowerSeries(f,bnd) : f in Basis(SN)];
   BN2 := [Evaluate(f,q^2) + O(q^bnd): f in BN];
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;

   // Put them together.
   B := B2N cat BN cat BN2;

   // Compute matrix of W2 with respect to B
   if Verbose then
      printf "Compute matrix of W2 with respect to B. ";
   end if;
   W2 := MatrixAlgebra(RationalField(),Dimension(S2N))!0;
     // (1) the B2N part.
   U2 := HeckeOperator(NewSubspace(S2N,2),2);
   for i in [1..#B2N] do
      for j in [1..#B2N] do
         W2[i,j] := -U2[i,j];
      end for;
   end for;
     // (2) the BN cat BN2 part.
   for i in [1..#BN] do
      W2[i+#B2N,#BN+#B2N+i] := 2;
      W2[#BN+#B2N+i,i+#B2N] := 1/2;
   end for;
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;
   
   // Compute a saturated integral basis BZ for S2N
   if Verbose then
      printf "Compute a saturated integral basis BZ for S2N to precision ", bnd, " ";
   end if;
   BZ := [PowerSeries(f,bnd) : f in Basis(S2N)];
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;

   // Compute the change-of-basis matrix from B to BZ.
   if Verbose then
      printf "Compute the change-of-basis matrix from B to BZ. ";
   end if;
   S := MatrixAlgebra(RationalField(),Dimension(S2N))!
           (&cat[Eltseq(S2N!f) : f in B]);
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;

   // Compute the matrix of W2 with respect to BZ.
   if Verbose then
      printf "Compute the matrix of W2 with respect to BZ. ";
   end if;
   W2Z := S^(-1)*W2*S;   
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;
   // print "W2Z = ", W2Z;
   // print "W2 = ", W2;

   // Compute T^*/(T^* meet WT^*) 
   if Verbose then
      printf "Compute T^*/(T^* meet WT^*). ";
   end if;
   V := VectorSpace(RationalField(),Dimension(S2N));
   Tstar := MakeLattice(Basis(V));
   W2Tstar := MakeLattice([W2Z[i] : i in [1..Nrows(W2Z)]]); 
   t2 := Cputime();
   intersect := (Tstar meet W2Tstar);
   if Verbose then
      printf "(%o s to compute intersection) ",Cputime(t2); t2 := Cputime();
   end if;
   M := Tstar/intersect;
   if Verbose then
      printf "(%o s to compute quotient) ",Cputime(t2); 
   end if;
   if Verbose then
      printf "(%o s).\n", Cputime(t);   t := Cputime();
   end if;
   if Verbose then
      print "M = ", M;
   end if;
   I := Invariants(M);
   M4 := AbelianGroup([GCD(4,a) : a in I | a ne 0]);
   M4modM2 := AbelianGroup([2 : a in I | a ne 0 and a mod 4 eq 0]);
   return #Invariants(M4modM2), M;
end intrinsic;


