R<x> := PolynomialRing(Integers());
C := HyperellipticCurve(x^3-1, x^6-x^5+5*x^4+6*x^2+5*x+1);
A := ModularSymbols("245H",+1);
S<y> := FieldOfFractions(PolynomialRing(Rationals()));

for p in [n : n in [2..500] | IsPrime(n) and n ne 5 and n ne 7] do
   Zeta := Numerator(ZetaFunction(BaseExtend(C,GF(p))));
   f := CharacteristicPolynomial(Tn(A,p));
   Hecke := Numerator(Evaluate(S!f, 1/y + p*y));
   "Zeta for p =", p, " is:";
   Zeta;
   if Eltseq(Zeta) ne Eltseq(Hecke) then
      print "PROBLEM with p = ", p;
   end if;
end for;
