/* **********************************************************************
   ** hecke.m                                                          **
   **                                                                  **
   ** Computing Bianchi modular forms and the Hecke action  for d=-2   **
   **                                                                  **
   ** Mehmet Haluk Sengun                                              **
   ** version of 13/08/2008                                            **
   **                                                                  **
   ********************************************************************** */

HECKE:=function(a,b,Data) 

K<w>:=QuadraticField(-2);

O:=MaximalOrder(K);
z:=O!w;

k:=Data`weight[1];
l:=Data`weight[2];

TA:=Data`ta;
TB:=Data`tb;
TU:=Data`tu;
CA:=Data`ca;
CB:=Data`cb;
CU:=Data`cu;

DER:=Data`der;
INN:=Data`inn;

s:=a+b*z;         /* this is the prime for the Hecke Operator*/

cs:=a-b*z;        /* this is the conjugate of s, also given by p/s */



p:=Norm(s,Integers());     /* norm over the integers */


c:=GF(p)! -(a/b) ;   /* this is image of w under reduction modulo p */


sl:=Data`level;
pl:=Norm(sl,Integers());  


Ml:=map<GF(pl) -> Integers() | x :-> x>;      /* this is to coerce ul(i)'s and vl(i)'s into Integers */

   

L,ff:=ResidueClassField(O,sl*O);
ms:=GF(pl)!ff(s);
mcs:=GF(pl)!ff(cs);

u:=map<GF(p)->GF(p) | x :-> -1/x>;          /* this is the permutation of B on the transversal */

v:=map<GF(p)->GF(p) | x:-> x-c >;        /* this is the permutation of U on the transversal */



ZZ:=map<GF(pl)->GF(pl) | x:-> ms*x>;



M:=map<GF(p) -> Integers() | x :-> x>;      /* this is to coerce u(i)'s and v(i)'s into Integers */



c:=M(c);


R<x,y>:=PolynomialRing(K,2);

Symm:=function(k,T)   /* returns the action of the matrix T on Symm(k) as a matrix */
  /* T:=Determinant(T)*T^-1; */
  ST:=Matrix(K,k+1,k+1,[<i,j,0>: i,j in [1..k+1]]);
  for i in [0..k] do Q:=(T[1,1]*x+T[1,2]*y)^(k-i)*(T[2,1]*x+T[2,2]*y)^(i);
  for j in [0..k] do ST[i+1,j+1]:=MonomialCoefficient(Q,x^(k-j)*y^(j));
  end for; end for;
  return ST;
end function;

D:=Matrix(K,2,2,[1,0,0,s]);  

D0:=Matrix(K,2,2,[1,0,0,cs]);   /* this is conjugate of \alpha^{\iota} */

SS:=MatrixAlgebra(K,(k+1)*(l+1));    
IDD:=SS!1; NN:=SS!0;

SDk:=Symm(k,D);
SD0l:=Symm(l,D0);


CD :=TensorProduct(SDk,SD0l);        /* action of delta inverse on the module */

level:=pl;
SSS:=MatrixAlgebra(K,(level+1)*(k+1)*(l+1));
N:=SSS!0;
ID:=SSS!1;

RRR:=[];
for i in [0..level-1] do
Append(~RRR,Integers()!Eltseq(s*i-Ml(ZZ(i)))[1]);
end for;

III:=[];
for i in [0..level-1] do
Append(~III,Integers()!Eltseq(s*i-Ml(ZZ(i)))[2]);
end for;

TD:=N;
for i in [0..level-1] do
TD:=InsertBlock(TD,CA^RRR[i+1]*CU^III[i+1]*CD,Ml(ZZ(i))*(k+1)*(l+1)+1,(i)*(k+1)*(l+1)+1);
end for;
TD:=InsertBlock(TD,CB*CD*CB,(level)*(k+1)*(l+1)+1,(level)*(k+1)*(l+1)+1);




TAA:=function(m)          

                         
W:=N;
if m lt 0 then            
for j in [0..(-m-1)] do
W:=W+TA^j;
end for;
W:= -W*TA^(m);
end if;

if m gt 0 then
for j in [0..(m-1)] do
W:=W+TA^j;
end for;
end if;

if m eq 0 then
W:=N;
end if;

return W;
end function;

/*=================================================================
FUNCTION 2    uses TU's
this will return 1+U+..+U^(k-1) or -U^k*(1+U+...+U^(-k-1)) if k is negative and 0 if k=0
=====================================*/
TUU:=function(m)          

                         
W:=N;
if m lt 0 then            
for j in [0..(-m-1)] do
W:=W+TU^j;
end for;
W:= -W*TU^m;
end if;

if m gt 0 then
for j in [0..(m-1)] do
W:=W+TU^j;
end for;
end if;

if m eq 0 then
W:=N;
end if;

return W;
end function;


/*============================================
FUNCTION  3  

word decomposition function, will work only for Euclidean fields 
this will be run in a for cycle with index i

======================================================*/
G:=function(i)               
S:=[];   /* this will return the decomposition */
         

B:=Matrix(O,2,2,[0,1,-1,0]);  /* we need to redefine B inside Matrix(O) */  

v:=O!M(u(i));

X:=Matrix(O,2,2,[v,(1+i*v)/s,-s,-i]);      /* this is the i'th matrix of the image of B under the transfer map*/

if Abs(Norm(X[1,2])) gt Abs(Norm(X[2,2])) then
X:=B*X;
Append(~S,O.1);   /* this will notify us that there's a B in the beginning of the decomposition, note that only in this case O.1=1, will be the first element of the sequence*/

end if;

while Norm(X[1,2])ne 0 do    /* this is the Euclidean Division */
r:=X[2,2] mod X[1,2];
q:=(X[2,2]-r)/X[1,2];
Append(~S,q);               
Q:=Matrix(O,2,2,[1,0,-q,1]);
X:=B*Q*X;
end while;

if X[2,2] eq 1 then Append(~S,X[2,1]);
end if;

if X[2,2] eq -1 then Append(~S,-X[2,1]);
end if;

return S;
end function;   

/* ==================================================== */
                     

RR<X,Y,Z>:=PolynomialRing(SSS,3);  /* X,Y,Z represent f(A), f(B), f(U), coefficients are matrices */

POLY:=RR!0; 
 
MOLY:=TAA(-2*b)*TU^a*TD*X+TUU(a)*TD*Z;  /* we need MOLY for image of U */


/* -----------------------BEGINNING of THE FOR CYCLE FOR IMAGE OF B ---------------------------*/


for i in [1..(p-1)] do      /* we will get the summation that defines the Hecke operator's image on B */

      
t:=#G(i);


 /* =================================================================
FUNCTION 4       (uses function1 G)

returns the real part of q_l, the l'th element of the word decomposition set S of the i'th matrix above 
============================================================= */
R:=function(l)     

h:=Eltseq(G(i)[l])[1];
h:=Integers()!h;
return h;
end function;

/* =================================================================
FUNCTION 5         (uses function1 G )

returns the imaginary part of q_l 
============================================================== */
I:=function(l)          

h:=Eltseq(G(i)[l])[2];
h:=Integers()!h;
return h;
end function;

/*===============================================================


this will return f(M_k) that appears in the image of B under the Hecke operator
t will be #Euc but we have to be careful about the first element of G(i)
==========================================================*/

               
                             

            
 
   


if G(i)[1] ne 1 then    /* this is the case where the decomposition does NOT start with a B*/

   g:=TAA(R(1))*TU^I(1)*X+TUU(I(1))*Z;    /*this is to define s_0 */

   for m in [2..t] do
   g:= (MonomialCoefficient(g,X)*TB*TA^R(m)*TU^I(m)+TAA(R(m))*TU^I(m))*X+
       (MonomialCoefficient(g,Y)*TB*TA^R(m)*TU^I(m)+TA^R(m)*TU^I(m))*Y+
       (MonomialCoefficient(g,Z)*TB*TA^R(m)*TU^I(m)+TUU(I(m)))*Z;
   end for;



else
   g:=TAA(R(2))*TU^I(2)*X+TA^R(2)*TU^I(2)*Y+TUU(I(2))*Z;    /*this is to define s_0 */
   for m in [3..t] do
   g:= (MonomialCoefficient(g,X)*TB*TA^R(m)*TU^I(m)+TAA(R(m))*TU^I(m))*X+
       (MonomialCoefficient(g,Y)*TB*TA^R(m)*TU^I(m)+TA^R(m)*TU^I(m))*Y+
       (MonomialCoefficient(g,Z)*TB*TA^R(m)*TU^I(m)+TUU(I(m)))*Z;
   end for;


end if;

   



POLY:=  POLY+MonomialCoefficient(g,X)*TD*TB*TA^i*X+
        MonomialCoefficient(g,Y)*TD*TB*TA^i*Y+ 
        MonomialCoefficient(g,Z)*TD*TB*TA^i*Z; 




/*this is the image of f(B) */
end for;


/* ---------------------------------BEGINNING OF THE FOR CYCLE FOR IMAGE OF U ----------------------------------*/







 /* ===============================================================
FUNCTION 6       

returns the real part of the upper corner of the k'th summand in image of U 
======================================================*/
RE:=function(m)     

h:=Eltseq((M(v(m))-m+z)/s)[1];
h:=Integers()!h;
return h;
end function;

/* ===============================
FUNCTION 7         
returns the imaginary part 
==================================== */
IM:=function(m)          

h:=Eltseq((M(v(m))-m+z)/s)[2];
h:=Integers()!h;
return h;
end function;

for j in [0..p-1] do

MOLY:=MOLY+(TAA(RE(j))*TU^IM(j)*TB*TD*TB*TA^j)*X+(TA^RE(j)*TU^IM(j)*TB+ID)*(TD*TB*TA^j)*Y+(TUU(IM(j))*TB*TD*TB*TA^j)*Z;
end for;

/*------------------------------END OF THE FOR CYCLE FOR IMAGE OF U -------------------------------------- */








H11:=TAA(a)*(TU^-b*TB*TD*TB+TU^b*TD);                 /* this will be the block at index (1,1) */
H21:=(TA^a*TU^-b*TB+ID)*(TD*TB);      
H31:=TUU(-b)*TB*TD*TB+TUU(b)*TD;


C1:=VerticalJoin(H11,H21);
C1:=VerticalJoin(C1,H31);



H12:=MonomialCoefficient(POLY,X);
H22:=MonomialCoefficient(POLY,Y);
H32:=MonomialCoefficient(POLY,Z);

C2:=VerticalJoin(H12,H22);
C2:=VerticalJoin(C2,H32);      /* this is the second column of the matrix of the Hecke operator on E^3 */




H13:=MonomialCoefficient(MOLY,X);
H23:=MonomialCoefficient(MOLY,Y);
H33:=MonomialCoefficient(MOLY,Z);
C3:=VerticalJoin(H13,H23);
C3:=VerticalJoin(C3,H33);           /* this is  column 3 of the matrix that gives the action of the Hecke opearator given by a+bw */


                                  

HO:=HorizontalJoin(C1,C2);
HO:=HorizontalJoin(HO,C3); 




V:=VectorSpace(K,3*(level+1)*(k+1)*(l+1));
DDER:=sub<V|Basis(DER)>;
BB:=ExtendBasis(DDER,V);
V:=VectorSpaceWithBasis(BB);

H:=Coordinates(V,BB[1]*HO);
HH:=V![H[t]:t in [1..3*(level+1)*(k+1)*(l+1)]];

for t in [2..3*(level+1)*(k+1)*(l+1)] do
H:=Coordinates(V,BB[t]*HO);
e:=V![H[r]:r in [1..3*(level+1)*(k+1)*(l+1)]];
HH:=VerticalJoin(HH,e);
end for;

HH:=Transpose(HH);

for i in [1..3*(level+1)*(k+1)*(l+1)-Dimension(DER)] do
HH:=RemoveRowColumn(HH,Dimension(DER)+1,Dimension(DER)+1);
end for;

IINN:=sub<DDER|Basis(INN)>;
CC:=ExtendBasis(IINN,DDER);
DDER:=VectorSpaceWithBasis(CC);

Y:=Coordinates(DDER,ReduceVector(IINN,CC[Dimension(INN)+1]*HO));
YY:=VectorSpace(K,Dimension(DER))![Y[t]:t in [1..Dimension(DER)]];

 for y in [Dimension(INN)+2..Dimension(DER)] do
 Y:=Coordinates(DDER,ReduceVector(IINN,CC[y]*HO));
 e:=VectorSpace(K,Dimension(DER))![Y[t]:t in [1..Dimension(DER)]];
 YY:=VerticalJoin(YY,e);
 end for;

 YY:=Matrix(YY);

 if Dimension(DER)-Dimension(INN) eq 1 then

 for i in [1..Dimension(INN)] do
 YY:=RemoveColumn(YY,1);
 end for;
 
 else

 YY:=Transpose(YY);
 for i in [1..Dimension(INN)] do
 YY:=RemoveRow(YY,1);
 end for;
 end if;

 EV:=SetToSequence(Eigenvalues(YY));
 ES:=[**];   /* creates a list */
 
 for i in [1..#EV] do
 Append(~ES,Eigenspace(YY,EV[i][1]));
 Append(~ES,EV[i][1]);
 end for;



 return FactoredCharacteristicPolynomial(YY),Eigenvalues(YY),ES;


 end function;