/*
  This code can be loaded, or copied and pasted, into Magma.
  It will load the data associated to the BMF, including
  the field, level, and Hecke and Atkin-Lehner eigenvalue data.
  At the *bottom* of the file, there is code to recreate the
  Bianchi modular form in Magma, by creating the BMF space
  and cutting out the corresponding Hecke irreducible subspace.
  From there, you can ask for more eigenvalues or modify as desired.
  It is commented out, as this computation may be lengthy.
*/

P<x> := PolynomialRing(Rationals());
g := P![6, -1, 1];
F<a> := NumberField(g);
ZF := Integers(F);

NN := ideal<ZF | {9, 3*a + 6}>;

primesArray := [
[2,a],
[2,a-1],
[3,a],
[3,a-1],
[13,a+4],
[13,a-5],
[23,a+11],
[5,5],
[29,a+10],
[29,a-11],
[31,a+7],
[31,a-8],
[41,a+15],
[41,a-16],
[47,a+13],
[47,a-14],
[7,7],
[59,a+26],
[59,a-27],
[71,a+20],
[71,a-21],
[73,a+29],
[73,a-30],
[101,a+24],
[101,a-25],
[11,11],
[127,a+27],
[127,a-28],
[131,a+48],
[131,a-49],
[139,a+16],
[139,a-17],
[151,a+42],
[151,a-43],
[163,a+62],
[163,a-63],
[167,a+77],
[167,a-78],
[173,a+45],
[173,a-46],
[179,a+65],
[179,a-66],
[193,a+19],
[193,a-20],
[197,a+59],
[197,a-60],
[211,a+35],
[211,a-36],
[223,a+36],
[223,a-37],
[233,a+103],
[233,a-104],
[239,a+53],
[239,a-54],
[257,a+78],
[257,a-79],
[269,a+120],
[269,a-121],
[271,a+46],
[271,a-47],
[277,a+81],
[277,a-82],
[17,17],
[307,a+49],
[307,a-50],
[311,a+70],
[311,a-71],
[317,a+75],
[317,a-76],
[331,a+44],
[331,a-45],
[347,a+164],
[347,a-165],
[349,a+158],
[349,a-159],
[353,a+37],
[353,a-38],
[19,19],
[397,a+84],
[397,a-85],
[409,a+28],
[409,a-29],
[439,a+118],
[439,a-119],
[443,a+51],
[443,a-52],
[449,a+169],
[449,a-170],
[461,a+163],
[461,a-164],
[463,a+74],
[463,a-75],
[487,a+112],
[487,a-113],
[491,a+153],
[491,a-154],
[499,a+31],
[499,a-32],
[509,a+110],
[509,a-111]];
primes := [ideal<ZF | {F!x : x in I}> : I in primesArray];

heckePol := x;
K := Rationals(); e := 1;

heckeEigenvaluesList := [*
1,
-1,
-1,
0,
2,
2,
0,
6,
2,
-2,
8,
8,
6,
-6,
0,
0,
-2,
4,
-4,
-8,
8,
-6,
-6,
6,
-6,
-22,
16,
16,
12,
-12,
-4,
-4,
0,
0,
4,
4,
8,
-8,
-18,
18,
-20,
20,
-18,
-18,
26,
-26,
-4,
-4,
8,
8,
10,
-10,
8,
-8,
14,
-14,
-18,
18,
16,
16,
22,
22,
-18,
-12,
-12,
-16,
16,
-2,
2,
20,
20,
12,
-12,
-14,
-14,
-30,
30,
22,
14,
14,
-10,
-10,
-8,
-8,
36,
-36,
18,
-18,
18,
-18,
-32,
-32,
-8,
-8,
-12,
12,
-12,
-12,
34,
-34
*];

heckeEigenvalues := AssociativeArray();
for i in [1..#heckeEigenvaluesList] do
    heckeEigenvalues[primes[i]] := heckeEigenvaluesList[i];
end for;

ALEigenvalues := AssociativeArray();
ALEigenvalues[ideal<ZF | {3, a}>] := 1;
ALEigenvalues[ideal<ZF | {a + 2, 3}>] := -1;

// EXAMPLE:
// pp := Factorization(2*ZF)[1][1];
// heckeEigenvalues[pp];

print "To reconstruct the Bianchi newform f, type
  f, iso := Explode(make_newform());";

function make_newform();
 M := BianchiCuspForms(F, NN);
 S := NewSubspace(M);
 // SetVerbose("Bianchi", 1);
 NFD := NewformDecomposition(S);
 newforms := [* Eigenform(U) : U in NFD *];

 if #newforms eq 0 then;
  print "No Bianchi newforms at this level";
  return 0;
 end if;

 print "Testing ", #newforms, " possible newforms";
 newforms := [* f: f in newforms | IsIsomorphic(BaseField(f), K) *];
 print #newforms, " newforms have the correct Hecke field";

 if #newforms eq 0 then;
  print "No Bianchi newform found with the correct Hecke field";
  return 0;
 end if;

 autos := Automorphisms(K);
 xnewforms := [* *];
 for f in newforms do;
  if K eq RationalField() then;
   Append(~xnewforms, [* f, autos[1] *]);
  else;
   flag, iso := IsIsomorphic(K,BaseField(f));
   for a in autos do;
    Append(~xnewforms, [* f, a*iso *]);
   end for;
  end if;
 end for;
 newforms := xnewforms;

 for P in primes do;
  if Valuation(NN,P) eq 0 then;
   xnewforms := [* *];
   for f_iso in newforms do;
    f, iso := Explode(f_iso);
    if HeckeEigenvalue(f,P) eq iso(heckeEigenvalues[P]) then;
     Append(~xnewforms, f_iso);
    end if;
   end for;
   newforms := xnewforms;
   if #newforms eq 0 then;
    print "No Bianchi newform found which matches the Hecke eigenvalues";
    return 0;
   else if #newforms eq 1 then;
    print "success: unique match";
    return newforms[1];
   end if;
   end if;
  end if;
 end for;
 print #newforms, "Bianchi newforms found which match the Hecke eigenvalues";
 return newforms[1];

end function;