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

  The basis for sign = +1 is 

     Bplus = { [j] for j in F_p }  together with  { [j] + [j^p]  for j not in F_p }

  and for sign = -1 it is 

     Bminus= { [j] - [j^p]  for j not in F_p }.

  Suppose T is matrix of operator with respect to basis { [j] : j in ss }.
  For each j not in F_p, replace T by matrix got by either adding or 
  subtracting [j^p] row to [j] row.  Then:

  * To get matrix of T with respect to Bplus, choose one of [j] or [j^p] 
    for each j not in F_p.  Then delete the corresponding rows and columns
    for that choice. 

  * For Bminus, just delete all rows and columns correspodning to [j] 
    with j in F_p from the matrix wrt Bplus.


TODO:
 (1) Modify ModSS to return sparse matrices.
 (2) Optimize for the case when auxiliary level is trivial.
 (3) Basic arithmetic for sparse matrices (i.e., subtract off scalar)
 (4) Test and run for p <= 234431.

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

function RestrictMatrix(A, x)  // written by Allan Steel.
   F := BaseRing(Parent(A));
   if Type(x) eq SeqEnum then
      B := Matrix(F, x);
   else
      assert Type(x) in {ModTupRng, ModTupFld};
      if Dimension(x) eq 0 then
         return MatrixAlgebra(F,0)!0;
      end if;
      B := Matrix(F, BasisMatrix(x));
   end if; 
   R := MatrixAlgebra(F, Nrows(B)) ! Solution(B, B*A);
   return R;
end function;

function LinearCombinations(Comb, B) 
// Compute the linear combinations of the elements of B
// defined by the elements of Comb.
    n := #B;
    if n eq 0 then
      return B;
    end if;
    return [Parent(B[1])|&+[B[i]*v[i] : i in [1..n]] : v in Comb];
end function;

/*function KernelOn(A, B)
// Suppose B is a basis for an n-dimensional subspace
// of some ambient space and A is an nxn matrix.
// Then A defines a linear transformation of the space
// spanned by B.  This function returns the
// kernel of that transformation.
   if Type(B) eq ModTupFld then
      B := Basis(B);
   end if;
   n := #B;
   if n eq 0 then
      return sub<Universe(B)|0>;
   end if;
   assert Nrows(A) eq Ncols(A) and Nrows(A) eq n;
   K := Kernel(A);
   return VectorSpaceWithBasis(LinearCombinations(Basis(K),B));
end function;
*/

procedure SetSparseRow(~T, i, v)
   assert Type(T) eq MtrxSprs;
   assert Type(i) eq RngIntElt;

   for j in [1..Ncols(T)] do
      T[i,j] := v[j];
   end for;
end procedure;

intrinsic DeleteRows(A::., positions::SeqEnum) -> .
{}
   return RMatrixSpace(BaseRing(Parent(A)), Nrows(A) - #positions, Ncols(A))!
        &cat[Eltseq(A[i]) : i in [1..Nrows(A)] | not (i in positions)];
end intrinsic;

intrinsic DeleteRowsAndColumns(A::AlgMatElt, positions::SeqEnum) -> AlgMatElt
{}
   return MatrixAlgebra(BaseRing(Parent(A)), Nrows(A) - #positions)!
            DeleteRows(Transpose(DeleteRows(Transpose(A), positions)), positions);
end intrinsic;

intrinsic SparseDeleteRows(A::MtrxSprs, positions::SeqEnum) -> MtrxSprs
{}
   B := SparseMatrix(BaseRing(A), Nrows(A) - #positions, Ncols(A));   
   j := 1;
   for i in [1..Nrows(A)] do 
      if not (i in positions) then
         SetSparseRow(~B, j, A[i]);
         j := j + 1;
      end if;
   end for;
   return B;
end intrinsic;

intrinsic SparseDeleteRowsAndColumns(A::MtrxSprs, positions::SeqEnum) -> MtrxSprs
{}
   return SparseDeleteRows(Transpose(SparseDeleteRows(Transpose(A), positions)), positions);
end intrinsic;

intrinsic RationalBasisVecs(X::ModSS) -> SeqEnum
{}
   k := GF(Level(X));
   J := SupersingularPoints(X);
   return [i : i in [1..#J] | J[i][1] in k];
end intrinsic;

intrinsic NonrationalBasisVecs(X::ModSS) -> SeqEnum
{}
   p := Level(X);
   k := GF(p);
   J := SupersingularPoints(X);
   conjugates := [];
   answer := [];
   for i in [1..#J] do 
      j := J[i][1];
      if not (j in k) then
         if not (i in conjugates) then
            Append(~answer, i);
            jp := j^p;
            Append(~conjugates, Index(J, <jp, jp>));
         end if;
      end if;
   end for;
   return answer, conjugates;
end intrinsic;

intrinsic HeckeOperator(X::ModSS, sign::RngIntElt, n::RngIntElt, R::Rng) -> AlgMatElt
{Compute T_n on subspace of X on which Atkin-Lehner acts by sign.}
   printf "Computing Hecke operator T_%o on full space directly.\n",n;
   time T := SparseHeckeOperator(X,n);
   if sign eq 0 then
      return T;
   end if;
   rat := RationalBasisVecs(X);
   nonrat, conjugates := NonrationalBasisVecs(X);
   if sign eq +1 then
      for i in [1..#nonrat] do
         // T[conjugates[i]] := T[conjugates[i]] - T[nonrat[i]];
         SetSparseRow(~T, conjugates[i], T[conjugates[i]] - T[nonrat[i]]);
      end for;
      S := SparseDeleteRowsAndColumns(T,rat cat nonrat);
   else
      for i in [1..#nonrat] do
         // T[conjugates[i]] := T[conjugates[i]] + T[nonrat[i]];
         SetSparseRow(~T, conjugates[i], T[conjugates[i]] + T[nonrat[i]]);
      end for;
      S := SparseDeleteRowsAndColumns(T,nonrat);   
   end if;
   return ChangeRing(S, R);
end intrinsic;

intrinsic HeckeOperator(X::ModSS, sign::RngIntElt, W::., n::RngIntElt, R::Rng) -> AlgMatElt
{Compute T_n on subspace W of subspace of X on which Atkin-Lehner acts by sign.}
   Tn := Matrix(HeckeOperator(X,sign,n,R));
   require Nrows(Tn) eq Degree(W) : "Argument 3 must have degree=", Nrows(Tn);
   return SparseMatrix(RestrictMatrix(Tn,W));
end intrinsic;

intrinsic HasseInterval(prime::RngIntElt) -> SeqEnum
{}
   require IsPrime(prime) : "Argument 1 must be prime.";
   return [Ceiling(-2*Sqrt(prime))..Floor(2*Sqrt(prime))];
end intrinsic;

function MyKernel(A)
   USE_SPARSE := true;
   if not USE_SPARSE then
      return Kernel(A);  
   end if;
   return Kernel(SparseMatrix(A));
end function;

function MyEigenspace(A, a)
   assert Type(A) eq MtrxSprs;
   assert Nrows(A) eq Ncols(A);
   for i in [1..Nrows(A)] do 
      A[i,i] := A[i,i]-a;
   end for;
   return Kernel(A);
end function;

function MyRank(A, a)
   for i in [1..Nrows(A)] do 
      A[i,i] := A[i,i]-a;
   end for;
   return Rank(A);
end function;

intrinsic EllipticFactors(X::ModSS, sign::RngIntElt, K::Fld,
                          W::., start_prime::RngIntElt, stop_prime::RngIntElt) -> SeqEnum
{}
   require IsPrime(start_prime) : "Argument 4 must be prime.";

   if Dimension(W) eq 0 then
      return [];
   end if;
   if Dimension(W) eq 1 then
      return [W];
   end if;

   T := HeckeOperator(X, sign, W, start_prime, K);
   ans := [];
   for a in HasseInterval(start_prime) do 
      printf "Computing kernel of T_%o - (%o)\n", start_prime, a;
      time r := MyRank(T,a);
      print "rank = ", r;
      if r eq Nrows(T) then
         continue;
      end if;
      time ker := MyEigenspace(T,a);
      V := RSpaceWithBasis(LinearCombinations(Basis(ker),Basis(W)));
      for A in EllipticFactors(X, sign, K, V, NextPrime(start_prime), stop_prime) do
         Append(~ans, A);
      end for;
   end for;
   return ans;
end intrinsic;

intrinsic EllipticFactors(X::ModSS, sign::RngIntElt, K::Fld, stop_prime::RngIntElt) -> SeqEnum
{}
   T := HeckeOperator(X, sign, 2, K);
   ans := [];
   for a in HasseInterval(2) do
      time r := MyRank(T,a);
      print "rank = ", r;
      if r eq Nrows(T) then
         continue;
      end if;
      printf "Computing kernel of T_2 - (%o)\n", a;
      time V := MyEigenspace(T,a);
      for A in EllipticFactors(X, sign, K, V, 3, stop_prime) do
         Append(~ans, A);
      end for;
   end for;
   return ans;
end intrinsic;