/* Helper functions */

AsaiLFactor_helper := function(ap, psip, p, w0, kron);
  if kron eq 1 then
    ap1, ap2 := Explode(ap);
    pol1 := Polynomial([psip*p^w0,-ap1,1]);
    pol2 := Polynomial([psip*p^w0,-ap2,1]);
    M1 := MatrixWithGivenCharacteristicPolynomial(pol1);
    M2 := MatrixWithGivenCharacteristicPolynomial(pol2);
    pol := CharacteristicPolynomial(KroneckerProduct(M1,M2));
    pol := Polynomial(Reverse(Eltseq(pol)));
  elif kron eq -1 then
    ap := ap[1];
    pol := Polynomial([1,-ap,p^(2*w0)])*Polynomial([1,0,-p^(2*w0)]);
  else
    ap := ap[1];
    pol := Polynomial([psip*p^w0,-ap,1]);
    M := MatrixWithGivenCharacteristicPolynomial(pol);
    pol := CharacteristicPolynomial(M^2);
    pol := Polynomial(Reverse(Eltseq(pol))) * Polynomial([1,-psip*p^w0]);  
  end if;
  return pol;
end function;


AsaiLFactor := function(f, p, psid);
  // Given a Hilbert newform f and a rational prime, this returns the local
  // factor of the Asai lift of f at p.
  
  F := BaseField(Parent(f));
  DF := Discriminant(F);
  k := Weight(Parent(f));
  k0 := Max(k);
  w0 := k0 - 1;
  ZF := Integers(F);
  psip := KroneckerSymbol(psid, p);
  kron := KroneckerSymbol(DF, p);
  
  if kron eq 1 then
    Ps := [ pp[1] : pp in Factorisation(p*ZF) ];
    ap1 := HeckeEigenvalue(f, Ps[1]);
    ap2 := HeckeEigenvalue(f, Ps[2]);
    pol := AsaiLFactor_helper([ap1,ap2], psip, p, w0, 1);
  elif KroneckerSymbol(DF, p) eq -1 then
    ap := HeckeEigenvalue(f, p*ZF);
    pol := AsaiLFactor_helper([ap], 1, p, w0, -1);  // assumes psid is Dirichlet
  else
    pp := Factorisation(p*ZF)[1][1];
    ap := HeckeEigenvalue(f, pp);
    pol := AsaiLFactor_helper([ap], psid, p, w0, 0);
  end if;
  return pol;
end function;








exnum := 1;   // choose an example, numbered 1 through 16
prec := 20;

/* Initialization */
QQ := Rationals();
PolsQQ<T> := PolynomialRing(QQ);

/* Defining data for each example 
   <num, alphas, z, chid, dF, NN = HMF level, psi = HMF central character, eps = twist, bad factors> 
   all betas are 0, t = 1/z for magma
   L(HGM) = L(chi)*L(f over QQ(sqrt(dF)), Asai, eps) with bad primes specified, 
   HMF level is given by a generator of the level ideal,
   chid is the disc of quadratic field for Dirichlet character chi, so chid = 1 means Riemann zeta  */
_<gr> := NumberField(Polynomial([-1,1,1]));  // golden ratio
        // num   alphas                 z        chid  dF  NN   psi eps  bad factors
dat := [ [* 1,  [1/2,1/2,1/2,1/2,1/2], -1/2^2,     1,  5,   4,  1,  1,   [<2,4,1-4*T>] *], // 1
         [* 2,  [1/2,1/2,1/2,1/2,1/2], -2^2,       1,  5,   4,  1,  1,   [<2,4,1-4*T>] *], // 2
         [* 3,  [1/2,1/2,1/2,1/2,1/2], -1/2^10,    1,  41,  1,  1,  1,   [<2,0,(1-4*T)*(1+5*T+10*T^2+80*T^3+256*T^4)>] *], // 3
         [* 4,  [1/2,1/2,1/2,1/2,1/2], -2^10,      1,  41,  1,  1,  1,   [<2,0,(1-4*T)*(1+5*T+10*T^2+80*T^3+256*T^4)>] *], // 4
         [* 5,  [1/2,1/2,1/2,1/3,2/3], (3/4)^3,    1,  37,  1,  1,  1,   [<2,0,(1-4*T)*(1+6*T-96*T^3-256*T^4)>, <3,0,(1-9*T)*(1+11*T+84*T^2+891*T^3+6561*T^4)>] *], // 5
         [* 6,  [1/2,1/2,1/2,1/3,2/3], -3^3,       1,  28,  8,  1,  -4,  [<2,4,(1-4*T)*(1-16*T^2)>, <3,0,(1-9*T)*(1+8*T+30*T^2+648*T^3+6561*T^4)>] *], // 6
         [* 7,  [1/2,1/5,2/5,3/5,4/5], -5^5/2^8,   1,  69,  1,  1,  1,   [<3,1,(1-9*T)*(1+9*T)*(1+5*T+81*T^2)>, <5,0,(1-25*T)*(1+4*T-70*T^2+2500*T^3+390625*T^4)>, <23,1,(1-529*T)*(1+529*T)*(1-470*T+279841*T^2)>] *],  // 7
         [* 8,  [1/2,1/2,1/2,1/4,3/4], 1/2^4,      1,  60,  4,  3,  1,   [<2,2,(1-4*T)^2*(1 + 6*T + 16*T^2)>, <3,1,(1-9*T)*(1+9*T)*(1+2*T+81*T^2)>, <5,1,(1-25*T)*(1+25*T)*(1-2*T+625*T^2)>] *], // 8
         [* 9,  [1/2,1/3,2/3,1/4,3/4], -1/(2^4*3), 1,  12,  81, 1,  1,   [<2,2,(1-4*T)^3*(1+4*T)>, <3,5,(1-9*T)^2*(1+9*T)>] *], // 9
         [* 10, [1/2,1/3,2/3,1/4,3/4], -3^3/2^4,   1,  172, 4,  1,  1,   [<2,2,(1-4*T)^3*(1+4*T)>, <3,0,(1-9*T)*(1+14*T+102*T^2+1134*T^3+6561*T^4)>, <43, 1, (1-43^2*T)^2*(1+946*T+3418801*T^2)> ] *], // 10
         [* 11, [1/2,1/3,2/3,1/6,5/6], -(3/4)^6,   1,  193, 1,  1,  1,   [<2,0,(1-4*T)*(1-16*T^2)^2>] *], // 11
         [* 12, [1/2,1/3,2/3,1/6,5/6], (3/5)^6,    1,  76,  8,  1,  -4,  [<2,4,(1-4*T)*(1-16*T^2)>, <3,0,(1-9*T)*(1+5*T-24*T^2+405*T^3+6561*T^4)>] *], // 12
         [* 13, [1/2,1/3,2/3,1/6,5/6], -1/80^3,    5,  129, 0,  1,  0,   [] *], // 13
         [* 14, [1/2,1/4,3/4,1/6,5/6], -1/2^10,    12, 492, 0,  0,  0,   [] *], // 14
         [* 15, [1/2,1/8,3/8,5/8,7/8], 1/7^4,      28, 168, 0,  0,  0,   [] *], // 15
         [* 16, [1/2,1/2,1/2,1/3,2/3], 3^3*gr^15,  0,  0,   0,  0,  0,   [] *]  // 16
       ]; 

assert exnum notin [7,13,14,15,16];  // not known

// Grab data for chosen example
assert exnum eq dat[exnum][1];
alphas := dat[exnum][2];
t := 1/dat[exnum][3];  // t = 1/z for Magma
chid := dat[exnum][4];
if chid eq 1 then zetaorder := 1; else zetaorder := 0; end if;
F<b> := QuadraticField(dat[exnum][5]);
ZF := Integers(F);
dF := Discriminant(Integers(F));
NNgen := dat[exnum][6];
NN := NNgen*ZF;
psid := dat[exnum][7];
epsd := dat[exnum][8];
badprimes := dat[exnum][9];

// Define HGM
H := HypergeometricData(alphas, [0,0,0,0,0]);
LHt := LSeries(H, t : ZetaOrder := zetaorder, BadPrimes := badprimes, Precision := prec);
if zetaorder eq 1 then
  LHtprim := LHt/TateTwist(RiemannZeta(),-2);
else
  LHtprim := LHt/TateTwist(LSeries(KroneckerCharacter(chid)),-2);
end if;

Conductor(LHtprim);
CFENew(LHtprim);  // check functional equation for the HGM

/* The Euler factors for the HGM. */
pps := PrimesUpTo(200,F);
ps := Sort({@ Generators(pp meet Integers())[1] : pp in pps | Gcd(Norm(pp),dF*NNgen*Denominator(&*alphas)*Numerator(t)*Denominator(t)) eq 1 @});
Lps_full := [EulerFactor(LHt,p) : p in ps];
 
/* The Euler factors for the (usually irreducible) degree 4 part of the HGM. */
Lps := [ PolsQQ | EulerFactor(LHt, p)/(1 - KroneckerSymbol(chid,p)*p^2*T) : p in ps ];

/* Now "undo" the tensor product to find the minimal polynomials for HMF Hecke eigenvalues */
pols := [];
for m := 1 to #ps do
  if KroneckerSymbol(dF, ps[m]) ne 1 then continue; end if;  // this only works for split primes
  // We should implement the inert prime case, it's even easier; and also good ramified primes
  rts := [ 1/r[1] : r in Roots(Lps[m], pAdicField(ps[m], 200)) ];
  if not Sort([ Valuation(s) : s in rts ]) eq [0,1,3,4] then continue; end if;  // avoid supersingular
  psidm := KroneckerSymbol(psid,ps[m]);
  if #rts ne 0 then
    rts0 := Roots(Polynomial([ -psidm*rts[1]*rts[2]/ps[m], 0, 1 ]));  // here's where we change character
    if #rts0 ne 0 then
      alpha := rts0[1][1];
      pol := PowerRelation(alpha + psidm*ps[m]^3/alpha, 4);
      Append(~pols, <ps[m], pol*Evaluate(pol,-Parent(pol).1)>);  // only know it up to sign and conjugation
    end if;
  end if;
end for;

/* Local Euler factors for the Asai lift. */
// SetVerbose("ModFrmHil", true);
S := NewSubspace(HilbertCuspForms(F, NN, [2,4]));
Vfs := NewformDecomposition(S : Dimensions := [1,2]);  
  // too expensive to do in general
  // N.B. dimension means *relative* dimension, which changes based on the weight!  Not sure what that means in general
// fs := [* Eigenform(f) : f in Vfs *];   // giving error?!
Vfs_passed := [* *];
for Vf in Vfs do
  Vf_passed := true;
  for pol in pols do
    p := pol[1];
    for ppfact in Factorization(pol[1]*ZF) do
      Tpp := CharacteristicPolynomial(HeckeOperator(Vf,ppfact[1]));
      if Gcd(pol[2],Tpp) eq 1 then
        Vf_passed := false;
        break;
      end if;
    end for;
    if not Vf_passed then
      break;
    end if;
  end for;
  if Vf_passed then
    Append(~Vfs_passed, Vf);
  end if;
end for;
print Vfs_passed;
fs_passed := [* Eigenform(f) : f in Vfs_passed *]; 

/* Comparing the Tate twist of the Asai lift to the irreducible 4-part of the HGM. */
for f in fs_passed do
  AsaiLocalFacts := [ PolsQQ | AsaiLFactor(f, p, psid) : p in ps ]; 
  TT_Asai_lifts := [ PolsQQ | Evaluate(AsaiLocalFacts[l], KroneckerSymbol(epsd,ps[l])*T/ps[l]) : l in [1..#ps] ];
  if exnum eq 9 then
    TT_Asai_lifts := [ Evaluate(TT_Asai_lifts[l], KroneckerSymbol(epsd,ps[l])*T/ps[l]) : l in [1..#ps] ];
  end if;
  print TT_Asai_lifts eq Lps, [TT_Asai_lifts[l] eq Lps[l] : l in [1..#ps]];
end for;

