// Tpprog.m --- a program to compute conjectural p-adic valuations of
//              eigenvalues of T_p on spaces of cusp forms. Tested on
//              magma-2.7. Uses William Stein's modular symbols algorithms,
//              but only a very little (e.g. for computing dimensions
//              of spaces of cusp forms)
//
// By Kevin Buzzard (buzzard@ic.ac.uk)
//
// Please report any problems to buzzard@ic.ac.uk .
//
// Version 2.1, 10/8/2001.
//
// Version history:
//
// V2.1, minor tinker, rewriting some print statements for clarity, and
//       adding some more remarks.
//
// V2.0: complete rewrite, unobfuscating the code---V1 was very obscure.
//       Got rid of debugging info, and also disentangled various
//       routines.
//
// V1.2, 7/8/00, minor edit for p=2 because the naive notion
// of Gamma_0(N)-regular is definitely wrong.
//
// V1.1, 26/7/00, initial hack.
//
// program starts here.

// First a little routine constantsequence(n,r) that creates
// the constant sequence [r,r,...,r] of length n.

constantsequence := function(n,r)
  return [r : i in [1..n]];
end function;

// Next a routine which returns the number of cusps of X_0(N)

numberofcuspsGamma0 := function(N)
  return &+[EulerPhi(GreatestCommonDivisor(d,N div d)) : d in Divisors(N)];
end function;

// Now the main function, tpslopes(p,N,kmax), which returns
// a vector of length kmax whose k'th entry is my conjectural
// list of valuations of eigenvalues of T_p on S_k(Gamma_0(N)).

print"";
print "tpslopes(p,N,kmax) returns a vector of length kmax, whose k'th entry";
print "(1<=k<=kmax) is the conjectural sequence of valuations of eigenvalues";
print "of T_p on forms of level N, weight k, and trivial character.";
print "The conjecture is only valid if p doesn't divide N and if p is";
print "Gamma_0(N)-regular.";

print "Example of use:";
print"> c:=tpslopes(2,1,100);";
print"> c[100];";
print"[ 3, 9, 9, 17, 17, 26, 26, 32 ]";
print">                                ";
print"Hence I conjecture that the 2-adic valuations of the eigenvalues of";
print"T_2 on cusp forms of level 1 and weight 100 are [3,9,9,...,32],";
print"which indeed they are, as one can verify by an explicit computation.";

tpslopes := function(p,N,kmax:Verbose:=false);

  // First some trivial checks.

  assert N gt 0 and IsPrime(p) and N in Integers() and (N mod p ne 0)
                and kmax in Integers() and kmax ge 0;
  
  // Next initialise "slopes" vector.

  slopes := [ [Integers()|] ]; // That is, a sequence of length one whose
                               // unique entry is a sequence of integers of
                               // length 0 !!! There's no slopes at weight
                               // one, of course.

  for k:=2 to p+2 do
    Append(~slopes,constantsequence(DimensionCuspFormsGamma0(N,k),0));
  end for;

  sl2:=slopes[2]; // remember this for later.

  // Now we must modify weight 2, in a way that makes me think that
  // perhaps we'd be better off predicting slopes at level Np instead.
  // We append a zero for each p-new form.

  slopes[2]:=constantsequence(
              DimensionCuspFormsGamma0(N*p,2)-DimensionCuspFormsGamma0(N,2),0);

  // We're done initializing now, if p>2. If p=2 then we have to modify
  // weight 4. Note "cat" is the magma 'concatenate sequences' operator.

  if p eq 2 then
    slopes[4] := slopes[2] cat 
                  constantsequence(DimensionCuspFormsGamma0(N,4)-#slopes[2],1);
    Append(~slopes,[]); // nothing at weight 5.
  end if;

  // That's the end of the initialization procedure. Now do we have
  // to do any hard work? Not if kmax is small.

  if #slopes ge kmax then
    // We can return already; best fix slopes[2] first though.
    slopes[2]:=sl2;
    return [slopes[i]:i in [1..kmax]];
  end if;

  // OK, if we get this far then we have to do some recursive defining!

  m:=numberofcuspsGamma0(N);
  if p gt 2 then G<teich>:=DirichletGroup(p,CyclotomicField(p-1)); end if;
  // (we don't need teich if p=2)
  trivialcharmodN:=DirichletGroup(N,CyclotomicField(p-1))!1;

  // main loop

  for k:=#slopes + 1 to kmax do

    // easy case, k odd.

    if k mod 2 eq 1 then Append(~slopes,[]); continue k; end if;

    // If we get here then k is even. We now associate three
    // integers a,b,c to k, and what algorithm we use for conjecturing
    // slopes at weight k will depend on these integers.

    a:=1; while p^a lt k-1 do a+:=1; end while; a-:=1;

    b:=1; while p^a*b lt k-1 do b+:=1; end while; b-:=1;

    c:=1+(((k-2-p^a*b)*p) div (p^a));

    // a can be arbitrarily large, it's about log_p(k). On the
    // other hand, 1<=b<=p-1 and 1<=c<=p.

    // The algorithm we use is one of three.

    if b+c lt p then alg:=1;
    elif b lt p-1 then alg:=2;
    else alg:=3;
    end if;

    // Now we define a vector v, which will be the conjectural
    // slopes at weight k. Its definition depends on which case
    // we're in.
    
    if alg eq 1 then

      // ********************* CASE 1 ****************************

      // (remark: never happens if p=2)
    
      k1:=k-b*(p-1)*p^(a-1);
      k2:=k-(b-1)*(p-1)*p^(a-1)-(b+c-1)*2*p^(a-1);
      v1:=slopes[k1];
      v2:=slopes[k2];
      base:=p^a*b+p^(a-1)*(c-1)+1;
      e:=k-base;
      epsilon:=teich^(base-1)*trivialcharmodN;
      s:=1+DimensionCuspForms(epsilon,1+e);

      assert #v1 + #v2 ge s-1;

      if #v1 ge s-1 then
         v:=[v1[i]:i in [1..s-1]];
      else
         v:=v1 cat [e-v2[s-i]:i in [1+#v1..s-1]];
      end if;

      assert #v eq s-1;

      v cat:= constantsequence(m,e);

      // ******************* END OF CASE 1 ***********************

    elif alg eq 2 then

      // ********************* CASE 2 ****************************

      // (remark: also never happens if p=2)

      k1:=k-((b+1)*p^(a-1)*(p-1));
      k2:=k-(p^(a-1)*(p-1));
      v1:=slopes[k1];
      v2:=slopes[k2];
      base:=(b+1)*p^(a-1)*(p-1)+1;
      e:=k-base;
      s:=1+DimensionCuspFormsGamma0(N*p,1+e);
      s2:=(s-1) div 2;

      assert #v2 ge s2;

      if #v1 ge s-1 then
        v:=[v1[i]:i in [1..s-1]];
      elif 2*#v1 ge s-1 then
        v:=v1 cat [e-v1[s-i]:i in [#v1+1..s-1]];
      else 
        w:=[v2[i]:i in [#v1+1..s2]];
        for i in w do assert i le (e div 2); end for;
        // this line can go later **DEBUG**
        wb:=[e-1-w[#w+1-i]:i in [1..#w]];
        v1b:=[e-v1[#v1+1-i]:i in [1..#v1]];

        v:=v1 cat w;
        if s mod 2 eq 0 then v cat:= [e div 2]; end if;
        v cat:= (wb cat v1b);
      end if;

      assert #v eq s-1;

      if e eq 1 then v cat:= constantsequence(m-1,e);
      else v cat:= constantsequence(m,e);
      end if;

      // ******************* END OF CASE 2 ***********************

    else assert alg eq 3;

      // ********************* CASE 3 ****************************
  
      // messiest case.

      k1:=k-p^a*(p-1);
      k2:=k-p^(a-1)*(p-1);
      v1:=slopes[k1];
      v2:=slopes[k2];
      base:=p^a*(p-1)+1;
      e:=k-base;
      s:=1+DimensionCuspFormsGamma0(N*p,1+e);
      s2:=(s-1) div 2;

      assert #v2 ge s2;

      if #v1 ge s-1 then
        v:=[v1[i]:i in [1..s-1]];
      elif 2*#v1 ge s-1 then
        v:=v1 cat [e-v1[s-i]:i in [#v1+1..s-1]];
      else
        w0:=[v2[i]:i in [#v1+1..s2]];
        w:=[Min(w0[i]+1,e div 2):i in [1..#w0]];
        wb:=[e-1-w[#w+1-i]:i in [1..#w]];
        v1b:=[e-v1[#v1+1-i]:i in [1..#v1]];

        v:=v1 cat w;
        if s mod 2 eq 0 then v cat:= [e div 2]; end if;
        v cat:= (wb cat v1b);
      end if;

      assert #v eq s-1;

      if e eq 1 then v cat:= constantsequence(m-1,e);
      else v cat:= constantsequence(m,e);
      end if;

      // ******************* END OF CASE 3 ***********************

    end if;

    // v currently has length about s. We now append to v the theta-like
    // slopes coming from weight k3, and this is the end of the construction.
    //  In fact, if v is already sufficiently large, it's already the end.

    if #v ge DimensionCuspFormsGamma0(N,k) then
      v:=[v[i]:i in [1..DimensionCuspFormsGamma0(N,k)]];
    else // need to add more
      k3:=2*base-k;
      v3:=slopes[k3];
      v cat:= [v3[i]+e:i in [1..DimensionCuspFormsGamma0(N,k)-#v]];
    end if;

    assert #v eq DimensionCuspFormsGamma0(N,k);

    // Finished!

    Append(~slopes,v);

  if Verbose and k mod 10 eq 0 then print "up to",k; end if;
  
  end for; // end of main loop.

  // Finally, fix weight 2 vector.
  slopes[2]:=sl2;

  return slopes;

end function;

// Finally, a little program to check that a prime p is Gamma_0(N)-regular.

isregular:=function(p,N);

  if N mod p eq 0 then return false; end if;

  ok:=true;

  for k:=2 to (p+3 div 2) by 2 do
    S:=CuspidalSubspace(ModularSymbols(N,k,-1));
    f:=CharacteristicPolynomial(HeckeOperator(S,p));
    l:=NewtonSlopes(f,p);
    if #l ge 1 and l[#l] cmpeq Infinity() then ok:=false;break k;end if;
    v:=[i:i in l];
    if v ne constantsequence(#v,0) then
      ok:=false;
      break k;
    end if;
  end for;
  if ok eq false then return false;end if;
  
  if p gt 2 then return true;end if;

  // if p=2 then we've only checked k=2 and we need to check k=4 too.

  S:=CuspidalSubspace(ModularSymbols(N,4,-1));
  f:=CharacteristicPolynomial(HeckeOperator(S,2));
  l:=NewtonSlopes(f,2);
  v:=[i:i in l];
  a:=DimensionCuspFormsGamma0(2*N,2)-DimensionCuspFormsGamma0(N,2);
  b:=DimensionCuspFormsGamma0(N,4)-a;
  if v ne constantsequence(a,0) cat constantsequence(b,1) then
    ok:=false;
  end if;

  return ok;
end function;
