

/* general basis handling */


intrinsic Extract_Zp_basis (L : known_dim := 0) -> Any
{Input: L=list of matrices over a p-adic ring.
Option: known_dim=the dimension of the space spanned by L.
Output: Basis of the space spanned by L.
This function uses Nakayama`s lemma and works by computing an F_p-basis of the reductions.}
  MR := Parent(L[1]);
  d := Nrows(L[1]);
  R := CoefficientRing(MR);
  p := Prime(R);
  K := GF(p);
  MK := KMatrixSpace(K,d,d);
  LK := [MK!t : t in L];

  V := sub<MK | >;
  dim := 0;
  basis_indices := [];
  i := 1;
  while (i le #LK) and ((known_dim eq 0) or (dim lt known_dim)) do
    V1 := sub<MK | V, LK[i]>;
    dim1 := Dimension(V1);
    if dim1 gt dim then
      V := V1;
      dim := dim1;
      Append(~basis_indices,i);
    end if;
    i := i+1;
  end while;

  RS := RMatrixSpace(R,d,d);
  bas := [RS!(L[i]) : i in basis_indices];

  return basis_indices, bas;
end intrinsic;



/* compute the modular forms in a Z_p-orbit  as the duals of the Hecke algebra */

intrinsic Zp_forms (L : known_dim := 0) -> Any
{Input: L=list of Hecke operators over Z_p such that L[n] is the n-th Hecke operator.
Option: known_dim=dimension of the span of L (if known).
Output: out,bi, where
        out=list the elements of which form a basis of the space of modular forms described by the dual (with respect to the basis described by bi) of the algebra generated by the given Hecke operators. Any modular form is represented as a function N --> Z_p, sending n to a_n.
        bi=list the elements of which are indices i such that the L[i] form a Z_p-basis.} 

  bi,bas := Extract_Zp_basis(L : known_dim := known_dim); // gets indices such that the associated Hecke operators form a basis of the Hecke algebra
  RB := RMatrixSpaceWithBasis(bas); // this is the corresponding matrix space with basis, allowing us to compute coordinates

  out := [];
  for i := 1 to #bas do
    fn := function (n)
      return Coordinates(RB,RB!L[n])[i];
    end function;
    Append(~out,fn);
  end for;

  return out,bi;
end intrinsic;



/* functions testing higher congruences */


intrinsic function_val ( f :: Any : bound := 2000, prime_to := 1) -> RngIntElt
{Input: f=function N --> O, where O is a valuation ring.
Output: minimum of the valuations of f(n) for n starting from 1 up to the first where the evaluation of f(n) fails, or up to bound. Only n will be used that are prime to prime_to.}
  val := Valuation(Parent(f(1))!0);
  n := 1;
  while (val gt 0) and (n le bound)  and (Gcd(n,prime_to) eq 1) do
    try
      v := Valuation(f(n));
      if v lt val then val := v; end if;
    catch e
      return val;
    end try;
    n := n+1;
  end while;
  return val;
end intrinsic;



intrinsic belongs_to_mod_pm ( F :: SeqEnum, bi :: SeqEnum, g :: Any : bound := 2000, prime_to := 1) -> Any
{Input: F=list of functions N-->Z_p which are an echelonised basis, where the echelon part is given by bi[i].
        g=another function N-->O, where O is a valuation ring.
Output: <m,e> with maximum m such that modulo the uniformiser of O to the m-th power, g belongs to the space spanned by F; e is the ramification index of the coefficient ring of g.  Only coefficients will be compared the index of which is prime to prime_to. bound is an upper bound for the number of coefficients to be used.}
  h := function(n)
    a := g(n);
    R := Parent(a);
    for i := 1 to #bi do
      a := a - g(bi[i])*R!(F[i](n));
    end for;
    return a;
  end function;

  return <function_val(h : bound := bound, prime_to := prime_to), RamificationIndex(Parent(g(1)))>;
end intrinsic;


intrinsic belongs_to_mod_pm_lists ( F :: Tup, bi :: SeqEnum, g :: SeqEnum : bound := 2000, prime_to := 1) -> Any
{Input: F=list of lists representing functions N-->Z_p which are an echelonised basis, where the echelon part is given by bi.
        g=list representing a function N-->O, where O is a valuation ring.
Output: <m,e> with maximum m such that modulo the uniformiser of O to the m-th power, the function given by g belongs to the space spanned by by the functions in F; e is the ramification index of the coefficient ring of g. Only coefficients will be compared the index of which is prime to prime_to. bound is an upper bound for the number of coefficients to be used.}
  h := function(n)
    a := g[n];
    R := Parent(a);
    for i := 1 to #bi do
      a := a - g[bi[i]]*R!(F[i][n]);
    end for;
    return a;
  end function;

  return <function_val(h : bound := bound, prime_to := prime_to), RamificationIndex(Parent(g[1]))>;
end intrinsic;


intrinsic make_list ( F :: Any : bound := 0 ) -> SeqEnum
{Input: F=function with source N.
Output: a list of the values of F.
Option: If bound has a nonzero value, then only go up to that value. Otherwise go on until the evaluation fails.}
  L := [];

  n := 1;
  while (bound eq 0) or (n le bound) do
    try
      Append(~L,F(n));
    catch e
      return L;
    end try;

    n := n + 1;
  end while;

  return L;
end intrinsic;


intrinsic make_tup_of_lists ( F :: Any : bound := 0 ) -> Tup
{Input: F=list of functions with source N.
Output: a tuple the i-th entry of which is a list of the values of F[i].
Option: If bound has a nonzero value, then only go up to that value. Otherwise go on until the evaluation fails.}
  out := <>;
  for f in F do
    Append(~out, make_list(f : bound := bound));
  end for;

  return out;
end intrinsic;

// main function for the decomposition into Qpbar-eigenforms
intrinsic _deco ( alpha :: Any, L :: Any : use_char_poly := false) -> Any
{Input: alpha = a matrix restriction function; L= list of commuting Z_p-matrices.
Output: list of matrix restriction functions to the Q_p-eigenvectors.}
  // if no decomposition is requested:
  if IsEmpty (L) then return [alpha]; end if;

  // take out the first matrix M from the list; Lnew is the remaining list
  // in this function call, M is treated; Lnew is passed on for recursive treatment
  M := alpha(L[1]);
  Lnew := L[2..#L];

  d := Nrows(M);

  // if the matrix is already of size 1, then don't do anything, just continue
  if d eq 1 then
    return _deco (alpha, Lnew : use_char_poly := false);
  end if;

  F := FieldOfFractions(CoefficientRing(M)); // make sure coefficients are in a field, not a ring

  ch := CharacteristicPolynomial(M);

  if use_char_poly then
    mipo := [<Parent(ch)![-1,1],1>];
  else
    mipo := Factorisation(ch);  // factor the charpoly
  end if;

  out := []; // list for the output

  for fac in mipo do
    /* - for every irreducible factor of the charpoly, choose one of its roots;
       - compute the eigenspace with this eigenvalue
       - write a base change function restricting to the eigenspace
    */

    // create a field containing the roots of the factor of the charpoly
    Kloc := LocalField(F,fac[1]);
    // it seems better to work with the "classical" p-adic fields in Magma, rather than LocalFields
    K := RamifiedRepresentation(Kloc);

    repeat
      MK := MatrixAlgebra(K,d);
      Pol := PolynomialRing(K);

      r := Roots(Pol!fac[1]);

      c := r[1][1];
      newmat := MK!M - MK!c;

      W := Kernel(newmat);
      B := Basis(W);

      // if the precision is not high enough, it can happen that the kernel is zero
      // in that case, lower the precision by 10, and compute again until a nonzero kernel is found
      if #B eq 0 then 
        K := ChangePrecision(K,Precision(K)-10);
      end if;
    until #B gt 0;

    // create a function that, on input a matrix on the big space, returns its image for the action on W wrt the basis B
    // for some reason, the function Coordinates doesn't work; this is why I use Solution
    fn := function(s)
      K1 := K;
      prec := Precision(K1);
      while true do
        try
          MB := Matrix(B);
          T := alpha(s);

          /* for some reason, the Solution command sometimes crashes on elements in matrix algebras, but not in K-matrix spaces;
             this is why we stupidly change types here; no deep reason. */
          MB := KMatrixSpace(K1,Nrows(MB),Ncols(MB))!MB;
          T := KMatrixSpace(K1,Nrows(T),Ncols(T))!T;
          
          M := Matrix([Solution(MB,B[i]*T) : i in [1..#B]]);
          return M;
        catch e
          prec := prec-1;
          if prec lt 1 then
            print "Error occured!!!!!";
            return 0;
          end if;
          // Lowering the precision in order to find existing solution to system of linear equations
          K1 := ChangePrecision(K1,prec);
          // print "Lowering the precision when computing coefficient to", prec, use_char_poly;
        end try;
      end while;
    end function;

    // continue recursively with the remaining matrices, which will be restricted to the eigenspace just computed using fn
    A := _deco(fn, Lnew : use_char_poly := false);

    out := out cat A;
  end for;

  return out;
end intrinsic;

intrinsic Qpbar_eigenforms (L :: SeqEnum) -> Any
{Input: L=list of Hecke operators over Z_p such that L[n] is the n-th Hecke operator.
Output: list of common eigenvectors given as functions N --> O sending n to the corresponding eigenvalue of L[n].}

  // get basis of the algebra generated by the matrices in L
  _, bas := Extract_Zp_basis(L);

  // the identity function to start the recursion
  alpha := function(s) return s; end function;

  // list for the output
  out := [];

  // compute the common eigenspaces
  // this depends on the first operator in the basis to be the identity (this is always true if the commands from this package are used)
  dec := _deco(alpha,bas : use_char_poly := true);

  // for every common eigenspace, write down the eigenvalue function
  for d in dec do
    // the matrix d(L[n]) is scalar; write the function n --> scalar of d(L[n])
    fn := function(n) return d(L[n])[1,1]; end function;
    Append(~out,fn);
  end for;

  return out;
end intrinsic;



/* automated handling */

intrinsic WeakCong (L1 :: SeqEnum, L2 :: SeqEnum : prime_to := 1, bound := 2000) -> Any
{Input: L1,L2=lists of Hecke operators over Z_p such that Li[n] is the n-th Hecke operator.

Output: List of entries of the form <n,e>, where n means that there is a Qpbar-eigenform in the Hecke algebra generated by L2 that is congruent modulo pi^n to a modular form (not necessarily an eigenform) in the Zp-Hecke algebra generated by L1; n is maximal with this property; pi is a uniformiser of the p-adic coefficient field of the Qpbar-eigenform, which is of absolute ramification index e.
Option: prime_to indicates that in the comparison only Hecke operators with index prime to prime_to are used; bound is an upper bound for the number of coefficients to be compared.}
  print "computing Z_p-basis";
  F1,bi1 := Zp_forms(L1);
  print "computing Q_p^bar-eigenforms...";
  F2 := Qpbar_eigenforms(L2);
  print "starting comparisons...";

  orbit_list := [];  // initialise the output list

  // run through the Qp_orbits in F2
  for F in F2 do
    h := belongs_to_mod_pm(F1,bi1,F : bound := bound, prime_to := prime_to);

    print "getting ",h;

    Append(~orbit_list, h);
  end for;

  return orbit_list;
end intrinsic;

