//q-Expansions of Hilbert Modular Forms (HMF):
//1. Hilbert Series of even HMF
//2. Eisenstein series of weight 2 in q-Expansion
//3. q-Expansion of Hecke eigenform
//4. Some more functions
//5. q-Expansions of generators of HMF and basis of space of HMF of given weight (in progress)
//for real quadratic number fields with narrow class number one
//2016 David Drsthoff. 
//david.dursthoff@rwth-aachen.de

forward HilbertSeries;
forward HilbertEisensteinSeries;
forward SiegelZeta;
forward HilbertSigma;
forward MyqExpansion;
forward MyHeckeEigenvalue;
forward Monic;
forward q1q2Valuation;
forward q1q2Coefficient;
forward MyGaussian;
forward q1q2Sort;
forward MyProduct;
forward HilbertModularForms;
forward HeckeConjugate;

/////////////////////////////////////
//1. Hilbert series, see Gerhard van der Geer, Hilbert Modular Surfaces, Sringer 1988, p. 188
function HilbertSeries( D, precision ) //D>1 square-free. F=Q[sqrt(D)] real quadratic number field. Returns Hilbert series 1 + \sum_{k=1}^\infty dim M_2k( SL_2(Z_K) ), including non-sym. HMF!
 F:=QuadraticField(D);
 _<t>:=PolynomialRing(Rationals());
 N:=ideal<Integers(F)|1>;
 D:=Discriminant(Integers(F));
 h:=ClassNumber(F);
 if D eq 5 then
  HSeries := (1+t^10)/(1-t)/(1-t^3)/(1-t^5);
  return PowerSeriesRing(Rationals(),precision)!HSeries, HSeries; 
 end if;
 if D eq 12 then
  HSeries:= (1-t^6)*(1-t^8)/(1-t)/(1-t^2)/(1-t^3)/(1-t^4);
  return PowerSeriesRing(Rationals(),precision)!HSeries, HSeries; 
 end if;
 if D eq 8 then
  a3plus := 1;
 else
  if D mod 3 eq 0 then
   a3plus:=1/2*h*(-3*D);
  else
   if (D mod 9 eq 3 ) then
    a3plus:=4*h*(-D/3);
   else
    a3plus:=3*h*(-3*D);
   end if;
  end if;
 end if;  
 d1:=h+Dimension(HilbertCuspForms(F,N,[2,2]));
 d2:=h+Dimension(HilbertCuspForms(F,N,[4,4]));
 d3:=h+Dimension(HilbertCuspForms(F,N,[6,6]));
 chiYGamma := d1-ClassNumber(F)+1; //because dim M_2 = b_1 + 2, b_1 = chi(YGamma)+h-3
 zetaFm1 := 1/4*(d3-2*d2+chiYGamma+h-2/3*a3plus); 
 chi := d3-3*d2+3*chiYGamma+2*h-a3plus;
 b:=[ 1 , chiYGamma+h-3 , 4*zetaFm1-chi-1/3*a3plus-h+3, 4*zetaFm1+2/3*a3plus-2, 4*zetaFm1-chi-1/3*a3plus-h+3, chiYGamma+h-3 ,1 ];
 HSeries:= &+[ b[i+1]*t^i : i in [0..6] ] / (1-t)^2 / (1-t^3);
 return PowerSeriesRing(Rationals(),precision)!HSeries, HSeries; 
end function; 


/////////////////////////////////////
//2. The Eisenstein series of weight 2, see Rolf Mueller, Hilbertsche Modulformen und Modulfunktionen zu Q[sqrt(8)], Mathematische Annalen 266 (1984)
function HilbertEisensteinSeries( D , powerbound )//D>1 square-free. Computes Hilbert Eisenstein Series up to q1^powerbound.
 F:=QuadraticField(D); 
 zeta:=SiegelZeta(F);
 d:=Discriminant(Integers(F)); 
 PR<q1,q2>:=PolynomialRing(F,2);
 if (Floor(d-Sqrt(d)))/2 in Integers() then  //(n *c - m)/2 is power of q1, (n c - m)/2 should be positive. c is a constant
  c:=d - Floor(d-Sqrt(d)) +2; 
 else
  c:=d - Floor(d-Sqrt(d))+1;
 end if; 
 kappa:=Integers()!4*zeta^-1; 
 
 g2:=PR!1;
 for n in [1..powerbound] do//computing coefficients
  if n*D/2 in Integers() then
   g2+:= kappa * HilbertSigma(F,[0,n]) * q1^n * q2^(Integers()!(c*n/2));
  end if;
  for m in [m : m in [1..Floor(n*Sqrt(d))] | (m-n*d)/2 in Integers()] do
   g2+:= kappa * HilbertSigma(F,[m,n]) * q1^n * q2^(Integers()!(c*n/2-m/2) ) + kappa * HilbertSigma(F,[m,n]) * q1^n * q2^(Integers()!(c*n/2+m/2));
  end for;
 end for;
 return g2;
end function;

SiegelZeta:=function(F)  //F tot. real number field. Returns Zeta function over F evaluated at -1, accorduing to Siegel, see van de Geer, p.20.
 d:=Discriminant(Integers(F));
 bound:=Floor( Sqrt(d-4));
 zetaF:=0;
 for n in [-bound..bound] do
  x:=((d-n^2)/4);
  if x in Integers() then 
   Factors:=Factorisation(Integers()!x); 
   sigma1:=1;                    
   for I in Factors do
    InnerSum:=1 + &+[I[1]^i : i in [1..I[2]]];
    sigma1*:=InnerSum; 
   end for;
   zetaF+:=sigma1;
  end if;
  end for;
 return zetaF/60;
end function;

HilbertSigma:=function(F,mn); // mn=[m,n]. Return sigma_(k-1)(ZF,nu), see Mueller, p.87.
 ZF:=Integers(F); N:=ideal< ZF | 1>;
 if Discriminant(Integers(F))/4 in Integers() then
  Nu:=ideal<ZF | mn[1]/2 + mn[2]*F.1>; 
 else 
  Nu:=ideal<ZF | mn[1]/2 + mn[2]/2*F.1 >;
 end if;
 Factors:=Factorisation(Nu); 
 sigma:=1;                    
 for I in Factors do
  InnerSum:=1+ &+[Abs(Norm(I[1]^i)) : i in [1..I[2]]];
  sigma*:=InnerSum; 
 end for;
 return sigma;
end function; 


/////////////////////////////////////
//3. q-Expansion of Hecke eigenform, see MAGMA'S Hilbert modular forms package by Dembele and Voight
function MyqExpansion(D,k, eigenform, powerbound)  //returns q1q2 expansion of eigenform up to q1^powerbound 
 F:=QuadraticField(D);
 d:=Discriminant(Integers(F));
 if (Floor(d-Sqrt(d))) mod 2 eq 0  then
  c:=d - Floor(d-Sqrt(d))+2;
 else
  c:=d - Floor(d-Sqrt(d))+1;
 end if; //(n *c - m)/2 is power of q1, (n c - m)/2 should be positive.
 if d mod 4 eq 0 then
  eps:=1;       
 else 
  eps:=2;
 end if; //Discriminant determs formula for elements of inverse different
 PR1<q1,q2>:=PolynomialRing(F,2);
 s:=PR1!0;
 for n in [1..powerbound] do
  if n*d mod 2 eq 0  then
   ideal:=ideal<Integers(F) | n/eps*F.1>;
   a:=MyHeckeEigenvalue(F,k,eigenform,ideal);
   if a notin Integers() and not Parent(a) subset BaseRing(PR1) then
    PR1<q1,q2>:=PolynomialRing(Parent(a), 2);
   end if;
   s:=PR1!s + PR1!a* PR1.1^n * PR1.2^(Integers()!(c*n/2));
  end if;
  for m in [m : m in [1..Floor(n*Sqrt(d))] | (m-n*d)/2 in Integers()] do
   ideal:=ideal<Integers(F) | m/2+n/eps*F.1>;
   a:=MyHeckeEigenvalue(F,k,eigenform,ideal);
   if a notin Integers() and not Parent(a) subset BaseRing(PR1) then
    PR1<q1,q2>:=PolynomialRing(Parent(a), 2);
   end if;
   s:=PR1!s + PR1!a* PR1.1^n * PR1.2^(Integers()!(c*n/2+m/2) );
   
   ideal:=ideal<Integers(F) | -m/2+n/eps*F.1>; 
   a:=MyHeckeEigenvalue(F,k,eigenform,ideal);  //equal to a above iff s is symmetric
   if a notin Integers() and not Parent(a) subset BaseRing(PR1) then
    PR1<q1,q2>:=PolynomialRing(Parent(a), 2);
   end if;
   s:=PR1!s + PR1!a* PR1.1^n * PR1.2^(Integers()!(c*n/2-m/2)); 
  end for;
 end for;
 return Monic(s); 
end function;

function MyHeckeEigenvalue(F,k, eigenform , ideal ) // eigenform f, ideal<>1 in Z_F. Returns Hecke eigenvalue of f, see 11.2.. 
//Acording to [Dembele, Shimura78 (prop. 2.1)] right??
 if ideal eq ideal<Integers(F)|1> then
  return 1;  // Because MAGMA gives the norminated HEF - but HeckeEigenvalue doesnot work for ideals notequal prim ideals
 end if;
 factors:=Factorization(ideal);
 if #factors eq 1 then
  P:=factors[1];
  if P[2] eq 1 then
   return HeckeEigenvalue( eigenform , P[1] );
  else
   return HeckeEigenvalue(eigenform,P[1]) * MyHeckeEigenvalue(F,k,eigenform,P[1]^(P[2]-1)) - Norm(P[1])^(k-1) * MyHeckeEigenvalue(F,k,eigenform, P[1]^(P[2]-2));
  end if;
 end if;
 return &*[ MyHeckeEigenvalue(F,k, eigenform, a[1]^(a[2]) ) : a in factors ];
end function;

Monic:=function(f) //returns monic form of f
 if f eq 0 then return f; end if;
 ab:=q1q2Valuation(f);
 return f/q1q2Coefficient(f,ab);
end function;

q1q2Valuation:=function(f);  //f in P. Return Order of f
 if f eq Parent(f)!0 then 
  return [-1,-1]; 
 end if;
 m:=Monomials(f)[#Monomials(f)];
 return [Degree(Evaluate(m,[Parent(f).1,1])),Degree(Evaluate(m,[1,Parent(f).2]))];
end function;

q1q2Coefficient:=function( f, ab ) // f in C[[q1,q2]], return coefficient of q1^a q2^b
 if ab eq [-1,-1] then
  return 0;
 else
  return Coefficient( Coefficient(f , Parent(f).1, ab[1]) , Parent(f).2, ab[2]);
 end if;
end function;

/////////////////////////////////////
//4. Some other useful procedures / functions

procedure MyGaussian(~list) //returns basis f_1 < f_2 < ... of same dimension as <list>
 i:=1;
 while i le #list do  //1. reducing steps
  for j in [1..(i-1)] cat [(i+1)..#list] do
   list[i]:= Monic( list[i]-q1q2Coefficient(list[i],q1q2Valuation(list[j]))*list[j]);
  end for;
  if list[i] eq 0 then
   Remove(~list,i);
  else
   i:=i+1;
  end if;
 end while;
  i:=1;
 while i le #list do  //2. reducing steps
  for j in [1..(i-1)] cat [(i+1)..#list] do
   list[i]:= Monic( list[i]-q1q2Coefficient(list[i],q1q2Valuation(list[j]))*list[j]);
  end for;
  if list[i] eq 0 then
   Remove(~list,i);
  else
   i:=i+1;
  end if;
 end while;
 q1q2Sort(~list);
end procedure;

q1q2Sort:=procedure(~list); //list of HMFs. Returns sorted list Basis by valuation.  (bubblesort)
  i:=#list;
  while i gt 1 do
    for j in [1..(i-1)] do
      if q1q2Valuation(list[j]) gt q1q2Valuation(list[j+1]) then
        a:=list[j]; list[j]:=list[j+1]; list[j+1]:=a;
      end if;
    end for;  
    i:=i-1;
  end while;
end procedure; 

MyProduct:=function(list,powerbound);
 if #list eq 0 then
  return Parent(list[1])!0;
 end if;  //return *&list bis q1^b
 Coeffs:=Coefficients(&*[f : f in list],Parent(list[1]).1)[1..Minimum(powerbound+1,#Coefficients(&*[f : f in list],Parent(list[1]).1))];
 if #Coeffs eq 0 then
  return Parent(list[1])!0;
 else
  return &+[ Coeffs[i]*Parent(list[1]).1^(i-1) : i in [1..#Coeffs]];
 end if;
end function;

/////////////////////////////////////
//5. Generators and Basis of Hilbert modular forms. 
//!!!Works for Q[sqrt(5)] and Q[sqrt(2)]. Not for ather fields.  
function HilbertModularFormsRing( D ,  weightbound : powerbound:=3) //F=Q[sqrt(D)] real quadratic number field, powerbound is bound for power of q1. Returns minimal generator list (up to weight weightbound) of the algebra of (symmetric) HMF and a basis of the space of (sym.) HMF of weight weightbound. Only for h_F^+=1.
 local F,PR,ZF,d,N,g2,Gen,Bas,i,j;
 F:=QuadraticField(D); 
 ZF:=Integers(F);
 d:=Discriminant(ZF); zeta:=SiegelZeta(F);   //invariants
 N:=ideal< ZF| 1>; //level 1 (at least right now)
 PR<q1,q2>:=PolynomialRing(F,2);  //algebra of q-extensions, actually it is a power series ring. 
 K:=F; //field expension over which cusp forms are definend, is overwritten later  
 
 print "Compute Eisenstein series of weight 2";
 g2:=PR!HilbertEisensteinSeries( D , powerbound ); //computes ES of weight 2, given in q1q2 expansion
 Gen:=[<g2,2>]; Bas:=[g2]; //first generator of HMF ring. No cups forms of weight 2.
 
 for k in [4..weightbound by 2] do  //calculate Basis of M_k in q-expansions, find new generators of wight k
  print "Compute weight", k;
  A:=[a : a in CartesianProduct([[0..Floor(k/g[2])] : g in Gen]) | &+[a[i]*Gen[i,2] : i in [1..#Gen] ]  eq k and a[1] lt k/2  ];  //known cusp forms weight k) (linear independent if generators algebraic independent)
  Bas:=[ MyProduct( [PR!Gen[i,1]^a[i] : i in [1..#a]], powerbound) : a in A];
  MyGaussian(~Bas); //Gaussian algorithm to get nice basis, especially independent.
  MyGaussian(~Bas);
  if #Bas lt #A then
   print "Found relation";
  end if;
  dimknown:=#Bas; //dim of known part of S_k 
  //TODO: Hilbert series.
  S:=HilbertCuspForms(F,N,[k,k]);  //space of HCF by MAGMA's HMF Package
  dim:=Dimension(S);
  if dimknown eq dim then
   print "No new generators of weight", k;
   print "dim M_" cat IntegerToString(k) cat "=", 1+#Bas;
  else  //new generators appear
   print "New generator(s) of weight", k;
   for Hauptraum in NewformDecomposition(S) do  //increase field extension, new cusp forms probably have coefficients not in Q. Usually cancel out later.
    K:=CompositeFields( K, HeckeEigenvalueField(Hauptraum) )[1];
   end for;
   PR<q1,q2>:=PolynomialRing(K,2);  //overwrite polynomial ring  
   NewS:=[ PR!MyqExpansion(F,k,eigenform,powerbound) : eigenform in Eigenforms(S)]; 
   DeleteHeckePrecomputation(QuaternionOrder(S));
   ClearStoredModularForms(F);  //we donot need the HMF MAGMA package in weight k anymore 
  
   Bas:=[ PR!f : f in Bas]; Gen:=[ <PR!f[1],f[2]> : f in Gen]; //q expansions over new given field
   for s in Bas do  //Gaussian on new cusp forms with old cusp forms. So we maximize valuation of new cusp forms.
    NewS:=[ Monic(snew-q1q2Coefficient(snew,q1q2Valuation(s))*s) : snew in NewS ];  
   end for;
   MyGaussian(~NewS);  //maximize valuation of new cusp forms, remove zero forms
   MyGaussian(~NewS); //to be sure.

   if dimknown + #NewS ne dim then //testing
    print "Something went wrong. Did not found all HMFs. \n Extensions are not implemented.";
   end if;
   print "dim M_" cat IntegerToString(k) cat "=", 1+#Bas+#NewS;
   Gen cat:= [<snew,k> : snew in NewS ];  //append new cusp forms 
  end if;//dimknown eq dim
  if k eq weightbound then  //(gaussian on old basis elements)
   Bas cat:=NewS;
   Append(~Bas, MyProduct( [g2 : i in [1..k/2]], powerbound ) );
   MyGaussian(~Bas);
   MyGaussian(~Bas);
  end if;
 end for;  //k
 return Gen, Bas;
end function;