// Check the congruences

_t0 := Cputime();

if not assigned maxtrace then
  maxtrace := 96;
end if;

///// Handle command line parameters /////

if assigned check_trace and Type(check_trace) eq MonStgElt then
  mintrace := StringToInteger(check_trace);
  maxtrace := StringToInteger(check_trace);
end if;

if assigned maxtrace and Type(maxtrace) eq MonStgElt then
  maxtrace := StringToInteger(maxtrace);
end if;
 
if assigned mintrace and Type(mintrace) eq MonStgElt then
  mintrace := StringToInteger(mintrace);
end if;

if not assigned mintrace then
  mintrace := 1;
end if;

SetColumns(0); // avoid breaking lines

/////

F<sqrt5> := QuadraticField(5);
AssignNames(~F, ["α"]);
ZF := Integers(F);
p79 := (1+4*sqrt5)*ZF;
assert Norm(p79) eq 79;

H22 := HilbertCuspForms(F, p79, [2,2]);
S22 := NewformDecomposition(NewSubspace(H22));
assert #S22 eq 1;
assert Dimension(S22[1]) eq 1;
f79 := Eigenform(S22[1]);

QO := QuaternionOrder(H22); // use same order for caching
H24 := HilbertCuspForms(F, p79, [2,4] : QuaternionOrder:=QO);
S24 := NewformDecomposition(NewSubspace(H24));
assert #S24 eq 2;
assert Dimension(S24[1]) eq 1;
assert Dimension(S24[2]) eq 3;
h79 := Eigenform(S24[2]);
g79 := Eigenform(S24[1]);

K := BaseField(h79);
AssignNames(~K, ["β"]);
ZK := Integers(K);

// the unique prime above sqrt5
lambda := Factorization(sqrt5 * ZK)[1,1];
assert lambda^3 eq sqrt5 * ZK;
_, red5 := ResidueClassField(lambda);
assert Codomain(red5) eq GF(5);

// the ramified prime above 2
q1 := Factorization(2 * ZK)[1,1];
q2 := Factorization(2 * ZK)[2,1];
assert q1^2 * q2 eq 2*ZK;
_, red2 := ResidueClassField(q1);
assert #Codomain(red2) eq 4;

polF := MinimalPolynomial(F.1); AssignNames(~polF, ["X"]);
polK := MinimalPolynomial(K.1); AssignNames(~polK, ["X"]);

if mintrace eq 1 then
print  "";
printf "F = Q(%o) where %o is a fixed root of %o\n", F.1, F.1, polF;
printf "N is the ideal of F generated by <%o>, of norm %o\n",
       F!Generators(p79)[1], Norm(p79);
print  "";
print  "f₇₉ is the HMF of weight (2,2) level N with eigenvalues in Q (dimension 1)";
print  "";
print  "h₇₉ is the HMF of weight (2,4) level N with eigenvalues in K (dimension 3)";
print  "g₇₉ is the HMF of weight (2,4) level N with eigenvalues in F (dimension 1)";
print  "";
printf "K = F(%o) where %o is a fixed root of %o\n", K.1, K.1, polK;
print  "";
printf "Check congruences with maxtrace := %o\n", maxtrace;
print  "";
print "       | μₚ(f₇₉)     | μₚ(f₇₉)  |";
print "       |   μₚ(h₇₉)   |   μₚ(h₇₉)|";
print "       |     μₚ(g₇₉) |          |";
print "       |             |          |";
print " Nm(p) |    mod λ    |   mod q₁ |   μₚ(f₇₉)  μₚ(h₇₉)";
print "       |             |          |";
print "";
end if;

// List of Fourier indices with trace "a"
function coefficients_with_trace(a)
    // x = a*(1+sqrt5)/2 + b*sqrt5/5 is in sqrt5^(-1) iff a,b in Z
    // The bounds for b ensure x is totally positive
    return [ a*(1+sqrt5)/2 + b*sqrt5/5
           : b in [ Ceiling((-5-Sqrt(5))*a/2)
                   .. Floor((-5+Sqrt(5))*a/2) ] ];
end function;

// All primes required to check by trace
function primes_by_trace(maxtrace)
    coefs := [ coefficients_with_trace(a) : a in [1..maxtrace] ];
    norms := [ { Integers() | Norm(c*sqrt5) : c in cs } : cs in coefs ];
    prims := [ { p : p in PrimeDivisors(n), n in ns} : ns in norms ];
    all := {};
    for i in [1..#prims] do
        prims[i] diff:= all;
        all join:= prims[i];
    end for;
    return [ [pp : pp in PrimeIdealsOverPrime(F,p), p in Sort(Setseq(ps))]
           : ps in prims ];
end function;

primes := primes_by_trace(maxtrace);

for trace in [mintrace..maxtrace] do
_t1 := Cputime();

if #primes[trace] eq 0 then
printf "trace: %o (no new prime needed)\n", trace;
_t1 := _t0;//continue;
else
norms := { Norm(p) : p in primes[trace] };
printf "trace: %o (needs %o new primes with norms %o - %o)\n",
       trace, #primes[trace], Minimum(norms), Maximum(norms);
end if;

for p in primes[trace] do

  n := Norm(p);

  if not assigned dont_compute_eigenvalues then
    mu_f := HeckeEigenvalue(f79, p);
    mu_h := HeckeEigenvalue(h79, p);
    mu_g := HeckeEigenvalue(g79, p);
  else
    mu_f := 0; mu_h := 0; mu_g := 0;
  end if;

  mu_g_x := Integers()!mu_g@red5;
  if n mod 5 eq 4 and mu_g_x ne 0 then mu_g_x -:= 5; end if;

  printf "%5o  |  %o   %o  %2o  |   %o  %o   | %7o    %o\n",
    n,  mu_f@red5, mu_h@red5, mu_g_x,
        mu_f@red2, mu_h@red2,
        mu_f, mu_h;

  // Check congruence μₚ(h₇₉) ≡ μₚ(f₇₉) (mod λ)
  if p ne sqrt5*ZF then
    assert mu_h@red5 eq mu_f@red5;
  end if;

  // Check congruence μₚ(h₇₉) ≡ μₚ(f₇₉) (mod q₁)
  if p ne 2*ZF then
    assert mu_h@red2 eq mu_f@red2;
  end if;

  // Check congruence μₚ(h₇₉) ≡ 0 (mod 2)
  if p ne p79 then
    assert mu_h@red2 eq 0;
  end if;

  // Check congruence ah ≡ twist(ag) (mod 5)
  assert mu_h@red5 eq (mu_g*n)@red5;

end for;
  printf "time for trace %o: %o seconds\n\n", trace, Cputime(_t1);
end for;

if mintrace lt maxtrace then
if mintrace eq 1 then
printf "Computed with maxtrace := %o\n", maxtrace;
else
printf "Computed with mintrace := %o, maxtrace := %o\n", mintrace, maxtrace;
end if;
print "";
printf "Total time: %o seconds\n", Cputime(_t0);
end if;

quit;