intrinsic pPowerTorsion(N::RngIntElt) -> SeqEnum
{The ord_p(#(J_0(N)(F_p))), for p=2,3,5,7,11,13,17,19.  In fact,
this function computes ord_p(f(p+1)), where f is the characteristic
polynomial of T_p acting on S_2(Gamma_0(N)).  The interpretation 
of this quantity when p divides N is left to the user.}

  P := [2,3,5,7,11,13,17,19] ;

  M := CuspidalSubspace(ModularSymbols(N,2,+1));
  if Dimension(M) lt 1 then 
     return [0 : p in P];
  end if;

  ans := [];
  for p in P do
     f := CharacteristicPolynomial(HeckeOperator(M,p) : 
         Al:="Modular", Proof:=false);
     v,_:=Valuation(Evaluate(f,p+1),p);
     Append(~ans,v);
  end for;

  DisownChildren(M);

  return ans;
end intrinsic;

