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


OUT:= recformat <
  der               : ModTupFld,
  inn               : ModTupFld,
  ca                : Mtrx,
  cb                : Mtrx,
  cu                : Mtrx,
  ta                : Mtrx,
  tb                : Mtrx,
  tu                : Mtrx,
  weight            : SeqEnum,
  level             >;



DIM:=function(a1,b1,k,l) 

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

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




sl:=a1+b1*z;        /* this is the prime for the congruence subgroup */

pl:=Norm(sl,Integers());     /* norm of ls over the integers, the level */

cl:=GF(pl)! -(a1/b1) ;   /* this is image of w under reduction modulo pl */


L,ff:=ResidueClassField(O,sl*O);




ul:=map<GF(pl)->GF(pl) | x :-> -1/x>;          /* this is the permutation of B on the transversal for the congruence subgroup in PSL*/


vl:=map<GF(pl)->GF(pl) | x:-> x+cl >;        /* this is the permutation of U on the transversal for congruence subgroup in PSL*/

jl:=map<GF(pl)->GF(pl) | x :-> -x>;          /* this is the permutation of J on the transversal for the congruence subgroup in PSL*/



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


cl:=Ml(cl);


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;

A:=Matrix(K,2,2,[1,0,1,1]);
B:=Matrix(K,2,2,[0,-1,1,0]);
U:=Matrix(K,2,2,[1,0,w,1]);  
U0:=Matrix(K,2,2,[1,0,-w,1]);  
J:=Matrix(K,2,2,[-1,0,0,1]);

AB2:=A*B*A*B;     
AB1:=A*B;        

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

/* the followings are actions of the above matrices on Symm */

SBk:=Symm(k,B);
SBl:=Symm(l,B);

SAk:=Symm(k,A);
SAl:=Symm(l,A);

SAB1k:=Symm(k,AB1);   
SAB1l:=Symm(l,AB1);

SAB2k:=Symm(k,AB2);    
SAB2l:=Symm(l,AB2);

SUk:=Symm(k,U);
SU0l:=Symm(l,U0);

SJk:=Symm(k,J);
SJl:=Symm(l,J);

CAB:=TensorProduct(SAB1k,SAB1l);   /* action of g */

CAB2:=TensorProduct(SAB2k,SAB2l);    /* action of g^2 */

CB:=TensorProduct(SBk,SBl);          /* action of B */
CU :=TensorProduct(SUk,SU0l);         /* action of U */
CA:=TensorProduct(SAk,SAl);
CJ:=TensorProduct(SJk,SJl);  


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

TA:=N;
for i in [0..level-2] do
TA:=InsertBlock(TA,IDD,(i+1)*(k+1)*(l+1)+1,(i)*(k+1)*(l+1)+1);
end for;
TA:=InsertBlock(TA,CA^level,1,(level-1)*(k+1)*(l+1)+1);
TA:=InsertBlock(TA,CB*CA*CB,(level)*(k+1)*(l+1)+1,(level)*(k+1)*(l+1)+1);

/* action of A on the induced module */



TB:=InsertBlock(N,IDD,1,(level)*(k+1)*(l+1)+1);
for i in [1..level-1] do
TB:=InsertBlock(TB,CA^(-Ml(ul(i)))*CB*CA^i,(Ml(ul(i)))*(k+1)*(l+1)+1,(i)*(k+1)*(l+1)+1);
end for;
TB:=InsertBlock(TB,IDD,(level)*(k+1)*(l+1)+1,1);


/* action of B on the induced module */





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

TJ:=N;
for i in [0..level-1] do
TJ:=InsertBlock(TJ,CJ,Ml(jl(i))*(k+1)*(l+1)+1,(i)*(k+1)*(l+1)+1);
end for;
TJ:=InsertBlock(TJ,CJ,(level)*(k+1)*(l+1)+1,(level)*(k+1)*(l+1)+1);




TAB2:=TA*TB*TA*TB;     
TAB:=TA*TB; 

E0:=VerticalJoin(N,TB+ID);
E0:=VerticalJoin(E0,N);


E1:=VerticalJoin((TB*(TAB2+TAB+ID)),(TAB2+TAB+ID));
E1:=VerticalJoin(E1,N);


E2:=VerticalJoin((TA^-1*(ID-TU^-1)),N);
E2:=VerticalJoin(E2,((TA^-1-ID)*TU^-1));


E3:=VerticalJoin(N,((TU^-1*TB*TU+TU)*(TB*TU^-1*TB*TU+ID)));
E3:=VerticalJoin(E3,((ID-TU^-1*TB*TU)*(TB*TU^-1*TB*TU+ID)));

DER:= Kernel(E1) meet Kernel(E2);     
DER:=DER meet Kernel(E3); 
DER:=DER meet Kernel(E0);               /* space of derivations */


F:=HorizontalJoin((ID-TA),(ID-TB));
F:=HorizontalJoin(F,(ID-TU));

INN:=Image(F);    


J1:=VerticalJoin(ID+TA^-1*TJ,N);
J1:=VerticalJoin(J1,N);

J2:=VerticalJoin(N,ID-TJ);
J2:=VerticalJoin(J2,N);

J3:=VerticalJoin(N,N);
J3:=VerticalJoin(J3,ID+TU^-1*TJ);

JJ:=HorizontalJoin(J1,J2);
JJ:=HorizontalJoin(JJ,J3);


PGL1:=Kernel(JJ) meet DER;
PGL2:=Kernel(JJ) meet INN;

W0:=VerticalJoin(ID,N);
W0:=VerticalJoin(W0,N);
W1:=VerticalJoin(N,N);
W1:=VerticalJoin(W1,ID);

WW:=HorizontalJoin(W0,W1);

WDER:=Kernel(WW) meet DER;
WINN:=Kernel(WW) meet INN;

full:=Dimension(DER)-Dimension(INN);
cusp:=Dimension(WDER)-Dimension(WINN);


plus:=Dimension(PGL1)-Dimension(PGL2);
cuspplus:=Dimension(PGL1 meet WDER)-Dimension(PGL2 meet WINN);

Data:=rec<OUT|
der:=DER,
inn:=INN,
ca:=CA,
cb:=CB,
cu:=CU,
ta:=TA,
tb:=TB,
tu:=TU,
level:=sl,
weight:=[k,l]>;


print "full,cusp",[full,cusp], "plus,cusp-plus", [plus,cuspplus] ;

return Data;
end function;