// heckeop: a way to compute char polys of Hecke operators on huge spaces.
// 
// Version 1.1, 9/1/2001, original version.
//
// The function charpolyheckeop takes a level N, a prime p prime to N,
// and a weight k, and computes the char poly of the Hecke operator T_p at
// level Gamma_0(N) and weight k. The idea is that it computes mod q
// for lots of primes q, and then glues the results together via the
// CRT command.
// 
// Of course it would be easy to modify things so as to slightly relax
// the assumptions on the character and p not dividing N, but currently
// I'm just interested in p=59 and N=1.

Qx<x>:=PolynomialRing(Rationals());  

charpolyheckeop:=function(N,p,k:Timings:=false)

  maxeig:=2*p^((k-1)/2);

  if Timings then tt:=Cputime();end if;

  d:=DimensionCuspFormsGamma0(N,k);
  if d eq 0 then return PolynomialRing(Rationals())!1; end if;

  maxvec:=[Binomial(d,i)*maxeig^i:i in [1..d]];
  // i'th entry of maxvec is bound for coefft of x^(d-i).

  max:=Maximum(maxvec);   // this tells us just how far we have to go.
                          // Perhaps a slightly smarter program would
                          // check that we are in the 99.9% of cases
                          // (for high weight) where this max is 
                          // just the bound for the constant coefft,
                          // and compute the char poly for fewer primes,
                          // and just the det for all primes after this.
                         
  l:=2;pr:=1;v:=[];

  // l is the prime after the one we're up to. pr is the product of the
  // good primes we've had so far. v is the list of char polys mod p for
  // all good primes so far.

  // I am *mildly* concerned about precision. I want to make sure
  // that the product of the first n good primes really is at least 2*max.
  // But can I make pr a real variable known to a few DPs? This would
  // probably make things a lot quicker...

  // Now build up a table of char polys mod p for as many good primes as
  // possible, (i.e. avoid primes in the denominator of elts of T),
  // collecting the results as <p,charpoly mod p> in the vector v,
  // and continuing until the product of the primes is sufficiently
  // large that an intger t with |t|<=max is uniquely determined
  // by t mod this product.

  lstep:=20; // a prime printer thing, for if Timings is on.
  lb:=lstep;

  while pr le 2.0001*max do // Note that I didn't go for 2 because
                            // of rounding errors: e.g. consider
                            // the following "weird" magma output:
                            // > SquareRoot(2)^2000 + 1 lt 2^1000;
                            // true  
                            // caused of course by rounding errors.

    if Timings then
      if l gt lb then
        print "Up to l =",l,".";
        lb+:=lstep; // for prime printing.
        print "Time so far:",Cputime(tt),"seconds.";
        print "Log(2*max/pr) is",Log((2*max)/pr);
      end if;
    end if;
    C:=CuspidalSubspace(ModularSymbols(N,k,GF(l),1));
    if Dimension(C) eq d then // no spurious torsion, so compute.
      pr*:=l;
      Append(~v,<l,PolynomialRing(Integers())!CharacteristicPolynomial(HeckeOperator(C,p))>);
    else
      if Timings then print "l =",l,"has spurious torsion."; end if;
    end if;
    DisownChildren(C);
    C:=0;
    l:=NextPrime(l);
  end while;

  if Timings then
    print "Looped over sufficiently many primes now. All primes used were less than",l,". Time taken so far:",Cputime(tt),"seconds.";
    print "Now doing all",d,"of the CRTs.";
  end if;

  // v is certainly all the data we need. Now rebuild.
  // The calculation so far used to take about 50% of the total
  // time but from V1.4 we're now probably over half way there,
  // because we're computing the coefficients of the char poly more sensibly.

  // Note that we have computed far too much data to compute,
  // for example, the trace! We have far too many primes in general.
  // Rather surprisingly, tests on CRT seem to indicate it's better
  // to work with lots of small primes in CRT rather than fewer bigger ones.

  maxlist:=[<d-i,maxvec[i]>:i in [1..d]];
    // Note that maxvec[i] is a bound for the coefft of x^(d-i).

  Sort(~maxlist,func<a,b|a[2]-b[2]>); // Sorts the list out into
                                      // order of size of maxvec,
                                      // with smallest first.

  pr:=1; // Start counting again!
  m:=[e[1]:e in v]; // list of primes we used.
  coeffts:=[0:i in [0..(d-1)]]; // coeffts[i] will be coefft of x^(i-1).
  for j:=1 to #v do 
    p:=v[j][1];
    pr*:=p;
    while #maxlist gt 0 and pr gt 2.0001*maxlist[1][2] do // We can compute
                                                          // a coefft now!
      cf:=maxlist[1][1]; // That coefft.
      c:=CRT([Coefficient(v[i][2],cf):i in [1..j]],[m[i]:i in [1..j]]);
      coeffts[cf+1]:=c lt pr/2 select c else c-pr;
      // That last line is because we want the unique soln in [-pr/2,pr/2].
      // Now kill 1st elt of maxlist.
      maxlist:=Reverse(Prune(Reverse(maxlist)));
      if Timings then
        printf "Just evaluated coefft of x^%o and time taken so far is %o.\n",cf,Cputime(tt);
      end if;
    end while;
  end for;

  if #maxlist gt 0 then error("Oops, we missed a coefft."); end if;

  if Timings then
    print "Should be just about done now; time taken so far is",Cputime(tt);
  end if;
  return x^d+Parent(x)!coeffts;
end function;
