twist.m

freeze;

////////////////////////////////////////////////////////////////
//                                                            //
//  Recoginising twists of newforms by Dirichlet characaters  //
//                                                            //
//  Steve Donnelly, October 2009                              //
//                                                            //
////////////////////////////////////////////////////////////////

debug := false;

intrinsic Twist(f::RngSerPowElt, chi::GrpDrchAElt) -> RngSerPowElt
{Twist the power series f by the character chi.}
    prec := AbsolutePrecision(f);
    // R := Parent(f);
    // K := Compositum(BaseRing(R), Codomain(chi));
    // _<q> := PowerSeriesRing(K);
    _<q> := Parent(f);
    coeffs := AbsEltseq(f);
    twisted_coeffs := [chi(n) * coeffs[n+1] : n in [0..#coeffs-1]];
    return &+[twisted_coeffs[n+1] * q^n : n in [0..#twisted_coeffs-1]] + O(q^prec);
end intrinsic;

intrinsic ApplyAut(sigma::Map, vec::ModTupFldElt) -> ModTupFldElt
{Apply the automorphism sigma to the vector vec.}
    return Vector([sigma(x) : x in Eltseq(vec)]);    
end intrinsic;

intrinsic ApplyAut(sigma::Map, f::RngSerPowElt) -> RngSerPowElt
{Apply the automorphism sigma to the power series f.}
    prec := AbsolutePrecision(f);
    coeffs := AbsEltseq(f);
    _<q> := Parent(f);
    return &+[sigma(coeffs[n+1])*q^n : n in [0..#coeffs-1]] + O(q^prec);    
end intrinsic;

// This one is to use with ComplexConjugate
intrinsic ApplyAut(sigma::Intrinsic, f::RngSerPowElt) -> RngSerPowElt
{Apply the automorphism sigma to the power series f.}
    prec := AbsolutePrecision(f);
    coeffs := AbsEltseq(f);
    _<q> := Parent(f);
    return &+[sigma(coeffs[n+1])*q^n : n in [0..#coeffs-1]] + O(q^prec);    
end intrinsic;

// Given f in a FldNum, express f as a poly in the gens

function relation(f, gens)
  F := Parent(f);
  assert F eq Universe(gens);
  assert IsAbsoluteField(F);
  Pol := PolynomialRing(Rationals(), #gens);
  degs := [Degree(MinimalPolynomial(gens[1]))];
  S := sub<F| gens[1] >;
  for i := 2 to #gens do 
    if f in S then
      Append(~degs, 1);
      continue; 
    end if;
    Append(~degs, Degree(MinimalPolynomial(gens[i], S)) );
    S := sub<F| gens[1..i] >;
  end for;
  assert f in S;
  powers := [[Pol.i^n : n in [0..degs[i]-1]] : i in [1..#gens]];
  mons := [Pol| &*tup : tup in CartesianProduct(powers)];
  M := ZeroMatrix(Rationals(), #mons, Degree(F));
  for i := 1 to #mons do
    InsertBlock(~M, Vector(Eltseq(F!Evaluate(mons[i],gens))), i, 1);
  end for;
  coeffs := Solution(M, Vector(Eltseq(F!f)));
  return &+ [coeffs[i]*mons[i] : i in [1..#mons]];
end function;

// Return iso F1 -> F2 of FldNums defined by sending gens1 to gens2
// where each Fi is generated by gensi

function isomorphism(gens1, gens2)
  F1 := Universe(gens1);
  F2 := Universe(gens2);
  assert IsAbsoluteField(F1) and IsAbsoluteField(F2);
  images := [Evaluate(relation(F1.i, gens1), gens2) : i in [1..Ngens(F1)]];
  return iso< F1 -> F2 | images >;
end function;

function is_isomorphism(gens1, gens2)
  F1 := Universe(gens1);
  F2 := Universe(gens2);
  if AbsoluteDegree(F1) ne AbsoluteDegree(F2) then
    return false, _;
  end if;
  phi := isomorphism(gens1, gens2);
  if forall{a : a in Generators(F1) | AbsoluteMinimalPolynomial(a) 
                                   eq AbsoluteMinimalPolynomial(a@phi)} then
    return true, phi;
  else
    return false, _;
  end if;
end function;

function eigenvalue(M, l)
  return Coefficient(Eigenform(M, l+1), l);
end function;


// For new irreducible M, get (non-multi) associated space
function associated_newspace(M)
  if not assigned M`associated_new_space then
    M := NewformDecomposition(M)[1];
  end if;
  if Type(M`associated_new_space) eq ModFrm then
    M := M`associated_new_space;
  end if;
  if assigned M`associated_newform_space then // M is multi
    M := M`associated_newform_space;
  end if;
  return M;
end function;


intrinsic IsTwist(M1::ModSymA, M2::ModSymA, p::RngIntElt : Bound:=0)
       -> BoolElt, GrpDrchAElt
{
Given two spaces of modular symbols representing (Galois orbits of) newforms,
this returns true iff some Galois conjugate (over Q) of M2 is the twist of M1 
by a nontrivial Dirichlet character of p-power conductor.  
If true, it also returns such a character.
}
// Note: when M1 = M2, this returns true iff the space is a self-twist 
//       by a NON-TRIVIAL character

  require IsPrime(p) : "The third argument should be prime";
  require IsNew(M1,p) : "The first argument is not a newspace at p";
  require IsNew(M2,p) : "The second argument is not a newspace at p";
  require #NewformDecomposition(M1) eq 1 and #NewformDecomposition(M2) eq 1: 
    "The arguments must correspond to a single Galois-conjugacy classes of newforms";

  M1 := associated_newspace(M1);
  M2 := associated_newspace(M2);

  N1 := Level(M1);
  N2 := Level(M2);
  N := LCM(N1, N2);
  r := Valuation(N,p);
  N0 := N div p^r;
  k := Weight(M1);

  // wlog M1 has the smaller p-level 
  if Valuation(N1,p) gt Valuation(N2,p) then
    bool, chi := IsTwist(M2, M1, p);
    if bool then
      return true, chi^-1;
    else
      return false, _;
    end if;
  end if;
  
  if k ne Weight(M2) or N ne N2 then
    return false, _;
  end if;

  eps1 := DirichletCharacter(M1);
  eps2 := DirichletCharacter(M2); 
  c1 := Valuation(Conductor(eps1), p);
  c2 := Valuation(Conductor(eps2), p);
  u := Min(r div 2, r - c2);

  // Loeffler-Weinstein Theorem 4.1 + erratum says:
  // If f1^chi = f2 then either
  //   (i) N(chi) divides p^u, or
  //  (ii) chi = eps2 * chi' where N(chi') divides both N(eps2) and p^u

  if c1 eq c2 then 
    // need to consider (ii)
    u := Max(u, c2);
  end if;
  // TO DO: stronger necessary condition for (ii): 
  //        N(eps1*eps2) < c2 (for some conjugates of eps1, eps2)
  // TO DO: for (ii), only try chi of the form eps2 * chi'

  // When trying to prove f1^chi = f2, we must compare eigenvalues 
  // up to the Sturm bound for some space known to contain both forms. 
  // We get the following upper bound for the level of f1^chi.
  v1 := Valuation(N1, p);
  if v1 eq 0 then
    twist_p_power := 2*u;
  else
    twist_p_power := 2*u + v1 - 1; // TO DO: why is this true?
  end if;
  Np := p^Max(r, twist_p_power);
  // Sturm bound for a subspace of Gamma1(N0*Np) with given character and 
  // Atkin-Lehner eigenvalues (see Buzzard & Stein "A mod 5 approach to 
  // modularity of icosahedral Galois representations", Lemma 1.4 and Cor 1.7)
  sturm_bound := Floor( k/12*Index(Gamma0(N0*Np))/2^#PrimeDivisors(N0) );
  if Bound gt 0 then
    Bound := Min(Bound, sturm_bound); 
  else
    Bound := sturm_bound;
    if IsVerbose("ModularSymbols") or IsVerbose("RepLoc") then
      print "[IsTwist] checking up to Sturm bound", Bound;
    end if;
  end if;

  test_primes_coprime := [l : l in PrimesUpTo(Bound) | GCD(l,N) eq 1];
  test_primes         := [l : l in PrimesUpTo(Bound) | l ne p];

  n := Exponent(UnitGroup(quo<Integers()|p^u>)); // order of chi divides n

  // if not twists, we should find out quickly
  b := Min(50, #test_primes_coprime);
  for l in test_primes_coprime[1..b] do
    pol1 := AbsoluteMinimalPolynomial(eigenvalue(M1,l)^n);
    if Evaluate(pol1, eigenvalue(M2,l)^n) ne 0 then
if debug then "quick false"; end if;
       return false, _;
    end if;
  end for;

  F1 := BaseRing(Parent(Eigenform(M1,1)));
  F2 := BaseRing(Parent(Eigenform(M2,1)));
  // Need to consider all composites F of F1 and F2.  
  if F1 eq Rationals() then
    composites := [<F2, Coercion(F2,F2)>];
  elif F2 eq Rationals() then
    composites := [<F1, Coercion(F2,F1)>];
  else
    F2a := AbsoluteField(F2);
    // Too expensive to just factor the defining poly of F2 over F1,
    // so first identify a large common subfield E1.
    b := Min(20, #test_primes_coprime);
    e1 := [eigenvalue(M1,l)^n : l in test_primes_coprime[1..b]];
    e2 := [eigenvalue(M2,l)^n : l in test_primes_coprime[1..b]];
    inds := [i : i in [1..#e1] | e1[i] notin Rationals()];
    assert inds eq [i : i in [1..#e1] | e2[i] notin Rationals()];
    if #inds eq 0 then
      E1 := Rationals();
      E1toF2 := Coercion(E1,F2);
      E1inF2a := sub< F2a | >;
    else
      if #inds gt 1 then 
        inds := inds[1..1]; 
      end if;
      e1 := [e1[i] : i in inds];
      e2 := [e2[i] : i in inds]; assert forall{a : a in e2 | a ne 0};
      F1a := AbsoluteField(F1);
      E1 := sub< F1a | ChangeUniverse(e1,F1a) >;
      E2 := sub< F2a | ChangeUniverse(e2,F2a) >;
      bool, E1toF2 := is_isomorphism( ChangeUniverse(e1,E1), ChangeUniverse(e2,E2) ); 
      if not bool then
if debug then "not an iso"; end if;
        return false, _;
      end if;
      E1inF2a := sub< F2a | E1.1@E1toF2 >;
    end if;
    // define embeddings of F2 via its absolute field
    pol2a := MinimalPolynomial(F2a.1, E1inF2a);
    pol2a_over_F1 := Polynomial([F1| (F2!c) @@ E1toF2 : c in Coefficients(pol2a)]);
    composites := [* *];
    for tup in Factorization(pol2a_over_F1) do 
      F := ext< F1 | tup[1] : DoLinearExtension >;
      F2toF := Coercion(F2,F2a) * hom< F2a->F | F.1 >;
      Append(~composites, <F, F2toF>);
    end for;
  end if;

  for tup in composites do 
    F, F2toF := Explode(tup);
    // Check if the values al/bl define a character chi with the right properties
    // TO DO: avoid lising all the elements?
    // TO D0: only consider c for which c^2*chi1 is a conjugate of chi2
    // TO DO: supply zeta to DirichletGroup
    chars := [* AssociatedPrimitiveCharacter(c) 
              : c in Elements(DirichletGroup(p^u, F)) 
              | not IsTrivial(c) *];
if debug then 
 chars0 := chars;
 printf "Considering %o characters: ", #chars; 
end if;

    for l in test_primes do

if debug and #chars gt 0 then
 e1l := eigenvalue(M1,l);
 e2l := eigenvalue(M2,l);
 c := chars[1];
end if;

      chars := [* c : c in chars | c(l)*F!eigenvalue(M1,l) eq F2toF(eigenvalue(M2,l)) *];

if debug then printf "%o ", <l, #chars>; end if;

      if IsEmpty(chars) then
        continue tup;
      end if;
    end for;
if debug then printf "\n"; end if;
    return true, chars[1];
 end for;
if debug then printf "\n"; end if;

 return false, _;
end intrinsic;


intrinsic IsMinimalTwist(M::ModSymA, p::RngIntElt : Bound:=0)
       -> BoolElt, ModSymA, GrpDrchAElt
{
Given a Hecke-irreducible space of modular symbols of level N,
returns true if the associated newform is not a twist of an
eigenform belonging to a space of lower level by a Dirichlet
character of p-power conducter.  If false, also returns the space
associated to the minimal newform, and the Dirichlet character.
}
  require IsPrime(p) : "The second argument should be prime";
  require IsNew(M,p): "The given space is not a newspace at p";
  require #NewformDecomposition(M) eq 1:
    "The given space must correspond to a single Galois-conjugacy class of newforms";

  M := associated_newspace(M);

  k := Weight(M);
  N := Level(M);
  r := Valuation(N,p);

  if r eq 0 then
    return true, M, DirichletGroup(1)!1;
  end if;

  N0 := N div p^r;
  eps := DirichletCharacter(M);
  c := Valuation(Conductor(eps), p);
  u := Min(r div 2, r - c);
if debug then "u =", u; end if;

  // Loeffler-Weinstein Theorem 4.1 + erratum says:
  // If we can reduce the p-power in the level, we can do so by 
  // twisting by a chi_p of conductor dividing p^u

  Chip := FullDirichletGroup(p^u); // possible twisting characters
  chips := Exclude(Elements(Chip), Chip!1); // (nontrivial)

  test_primes := [l : l in PrimesUpTo(100) | GCD(l,N) eq 1][1..10];

  for n := 0 to r-1 do
    // collect characters psi = eps*chip^2 that have p-level dividing p^n 
    Psi := FullDirichletGroup(N0*p^n);
    psis := {Psi| };
    for chip in chips do
      bool, psi := IsCoercible(Psi, eps*chip^2);
      if bool then
        Include(~psis, psi);
      end if;
    end for;
  //psis := GaloisConjugacyRepresentatives(Setseq(psis)); // TO DO: not sure about this
    psis := [* MinimalBaseRingCharacter(psi) : psi in psis *];
    for psi in psis do
      Mn := NewSubspace(CuspidalSubspace(ModularSymbols(psi, k, 1)), p);
      assert Level(Mn) eq N0*p^n;
if debug then "n =", n; Mn; end if;
 
      // if M does not twist to a subspace of Mn, we should find out here
      // order of chip divides e
      e := 2*Order(psi/eps);
      for l in test_primes do
        al := eigenvalue(M, l);
        Tl := HeckeOperator(Mn, l);
        if Rank(Evaluate(AbsoluteMinimalPolynomial(al^e), Tl^e)) eq Nrows(Tl) then
          continue psi;
if debug then "Rejected"; end if;
        end if;
      end for; 
 
      MMs := NewformDecomposition(Mn : Sort:=false);
      Sort(~MMs, func<x,y|Dimension(x)-Dimension(y)>);
      for MM in MMs do 
if debug then MM; end if;
        bool, tau := IsTwist(MM, M, p : Bound:=Bound);
        if bool then 
          return false, MM, tau;
        end if;
      end for;
    end for; // psi
  end for; // n

  return true, M, DirichletGroup(1)!1;
end intrinsic;

//////////////////////////////   TESTING    /////////////////////////////////

// Check that twists show up at the expected levels (or below),
// for levels dividing N0*p^E and weight k.
// For each newform space f of level N1 and chi of conductor p^c,
// the twist f^chi has level at most:
//   N1 * p^(2*c)   if N1 is prime to p,
//   N1 * p^(2*c-1) otherwise

function TestIsMinimalTwist(k, N0, p, E : Bound:=500, max_dim:=50)
  assert GCD(N0, p) eq 1;
  
  // list <MM, chi, ... >  for which we expect MM^chi to have p-level dividing p^E
  expected_twists := [* *];
  printf "Determining expected twists: "; time
  for e := 0 to E-1 do
    Me := CuspidalSubspace(ModularSymbols(Gamma1(N0*p^e), k, 1));
    c := (e eq 0) select E div 2
                  else  (E-e+1) div 2;
    Chic := FullDirichletGroup(p^c);
    chis := [chi : chi in Elements(Chic) | not IsTrivial(chi)];
    MMs := [MM : MM in NewformDecomposition(Me) | IsNew(MM,p) and Dimension(MM) le max_dim];
    expected_twists cat:=
      [* < assigned MM`associated_newform_space select MM`associated_newform_space else MM, 
           chi, Conductor(chi), Order(chi)> : MM in MMs, chi in chis *];
  end for;
  for tup in expected_twists do tup; end for;

  self_twists := [* *];
  non_minimal_twists := [* *];
  same_level_twists := [* *];
  twists := [* *];

  printf "Self twists: "; time
  for tup in expected_twists do
    bool, chi := IsTwist(tup[1], tup[1], p : Bound:=Bound);
    if bool then
      Append(~self_twists, <tup[1], chi, Conductor(chi), Order(chi)> );
    end if;
  end for;
  for tup in self_twists do tup; end for;

  printf "Non-minimal twists: "; time
  for tup in expected_twists do
    time bool := IsMinimalTwist(tup[1], p : Bound:=Bound);
    if not bool then
      Append(~non_minimal_twists, tup[1]);
    end if;
  end for;
  for MM in non_minimal_twists do MM; end for; 

  printf "Twists at same level: "; time
  for i in [1..#expected_twists], j in [i+1..#expected_twists] do
    Mi := expected_twists[i][1];
    Mj := expected_twists[j][1];
    if Level(Mi) eq Level(Mj) then
      time bool, chi := IsTwist(Mi, Mj, p : Bound:=Bound);
      if bool then
        Append(~same_level_twists, <Mi, chi, Conductor(chi), Order(chi)> );
        Append(~same_level_twists, <Mj, chi^-1, Conductor(chi), Order(chi)> );
      end if;
    end if;
  end for;
  for tup in same_level_twists do tup; end for;

  for e := 1 to E do
    printf "Looking for twists that show up at p-level %o\n", p^e;
    printf "Computing spaces: "; time
    MMs := NewformDecomposition(CuspidalSubspace(ModularSymbols(Gamma1(N0*p^e), k, 1)) : Sort:=false);
    MMs := [MM : MM in MMs | IsNew(MM, p)];
    for MM in MMs do
      time bool, Mchi, chi := IsMinimalTwist(MM, p : Bound:=Bound);
      if not bool then
        print <AssociatedNewSpace(Mchi), chi, Conductor(chi), Order(chi)>;
        Append(~twists, <AssociatedNewSpace(Mchi), chi, Conductor(chi), Order(chi)>);
      end if;
    end for;
  end for;

  return expected_twists, self_twists, non_minimal_twists, same_level_twists, twists;
end function;

function norms(M, n)
  return [AbsoluteNorm(eigenvalue(M,l)) : l in PrimesUpTo(n)];
end function;