/* this has the functions we need for the T_l operator */

 /* l is the prime for the T_l operator; it's set by the user (or by 
quaternionic_modular_forms.g originally */

/* if noisyflag is set, the code that computes the coset representatives 
prints its output */

print("Computation of the T_l operator");

noisyflag=1;

if(noisyflag,print("Verbose printing of output for T_l computations"));

if(!isprime(l),error("l must be a prime number"))

if(l==p,error("The primes l and p must be different"))

print("The prime l is ",l)

/* this gives the values of x0 and y0 for the homom between D_l and M_2(\Q_l) */

if(l%4==1,XX_l=-1;y0_l=0);

if(l%8==3,XX_l=-2;y0_l=1);

if(l%8==7,XX_l=-5;y0_l=2); /* NB this may not work */

/* Define isomorphisms th: K(i,j) -> M_2(K) and th1: M_2(K) -> K(i,j) */

th_l(v)=[v[1]+v[2]*x0_l+v[4]*y0_l, v[2]*y0_l-v[3]-v[4]*x0_l; \
v[2]*y0_l+v[3]-v[4]*x0_l, v[1]-v[2]*x0_l-v[4]*y0_l]

th1_l(mm) = [(mm[1,1]+mm[2,2])/2, \
((mm[2,2]-mm[1,1])*x0_l - (mm[1,2]+mm[2,1])*y0_l)/2, \
(mm[2,1]-mm[1,2])/2,((mm[2,2]-mm[1,1])*y0_l+(mm[1,2]+mm[2,1])*x0_l)/2]

bound_l=floor(sqrt(l)+3);

forstep(d0=-bound_l, bound_l, 1/2, forstep(d1=-bound_l, bound_l, 1/2, \
forstep(d2=-bound_l, bound_l, 1/2, forstep(d3=-bound_l, bound_l, 1/2, \
if(d0^2+d1^2+d2^2+d3^2==l, \
mm = subst(th_l([d0,d1,d2,d3]), x0_l, sqrt(XX_l+O(l^ACC))); \
if((mm[1,1]%l == 0) && (mm[2,1]%l == 0), \
d1list_l=concat(d1list_l, [[d0,d1,d2,d3]])))))))

if(flag,stepsize=1/2,stepsize=1);

bigd1list_l=[];

forstep(d0=-bound_l, bound_l, stepsize, forstep(d1=-bound_l, bound_l, stepsize, \
forstep(d2=-bound_l, bound_l, stepsize, forstep(d3=-bound_l, bound_l, stepsize, \
if(d0^2+d1^2+d2^2+d3^2==l , \
bigd1list_l=concat(bigd1list_l, [[d0,d1,d2,d3]]))))))

matunitlist_l = vector(length(unitlist), i, th_l(unitlist[i]));

smallmatunitlist_l = vector(length(smallunitlist), i, th_l(smallunitlist[i]));

/* these test the matrices at l for membership in the respective groups 
U_*(l) */

test_l(mm)= MM = mm; temp = matrix(2,2,i,j, \
valuation(MM[i,j],l)); mmin = vecmin([vecmin(temp[1,]), \
vecmin(temp[2,])]); mmin >=0

testU0_l(n)=m=subst(n,x0_l,sqrt(XX_l+O(l^ACC)));test_l(m)

testU1_l(n)=m=subst(n,x0_l,sqrt(XX_l+O(l^ACC)));test_l(m)&&((m[2,1]%l)==0)&&((m[2,2]%l)==1)

test_at_2(n)=valuation(q_norm(subst(q_sum(n,[-1,0,0,0]),x0_l,XX_l+O(l^ACC))),2)>=3
/* kludge here; checking for membership in 1+m^3 */

test_at_2_general(n,t)=valuation(q_norm(subst(q_sum(n,[-1,0,0,0]),x0_l,sqrt(XX_l+O(l^ACC)))),2)>=t


vtl(t) = if(t!=l,[1,t; 0, l],[l,0;0,1])

fact1list_l = [];

/* note: i runs from 0 to length(clist)-1 [cosets]

         t runs from 0 to **l-1** */

/*for(t=0,l-1, for(i=0, length(clist)-1, for(j=1, length(bigd1list_l),\
u = th_l(bigd1list_l[j])*(1/vtl(t)); \
uu = subst(u, x0_l, sqrt(XX_l + O(l^ACC))); \
if(test_l(uu), fact1list_l = concat(fact1list_l, \
[[t,i,j,[quatinv(bigd1list_l[j]), u, vtl(t)]]]); break(1)))))*/


unitlist_H_l= vector(length(unitlist_H), i, subst(unitlist_H[i],x0,x0_l));

/* this code factors c_i v_t^-1 as d c_j u */

/* fl_l=[];for(a=0,length(clist)-1,for(b=0,l,test_list=[];for(i=1,length(bigd1list_l),\
current=subst(th_l(bigd1list_l[i])*vtl(b)^-1,x0_l ,-sqrt(XX_l +O(l^ACC)));\
current2=th_l(bigd1list_l[i])*vtl(b)^-1;\
if(testU0_l(current)&&testU0_l(current^-1),test_list=concat(test_list,[[bigd1list_l[i],current2*vtl(b)]])));\
for(i=1,length(test_list),for(j=0,length(clist)-1,\
if(testUstar(myc(j)^-1*th(test_list[i][1])*myc(a)),\
if(test_list[i][2][2,2]==test_list[i][2][2,2],\
if(noisyflag,print([a,b,i,j]));fl_l=concat(fl_l,[[a,b,i,j,myc(j)^-1*th(test_list[i][1])*myc(a)]]);break(2))))))) 
*/

/*fl_l=[];for(b=0,l,for(i=1,length(bigd1list_l),if(test_at_2(bigd1list_l[i]),print(i),)))

my_quat=[1, -2, 0, 0];

fl_l=[];for(b=0,l,fl_l=concat(fl_l,[[0,0,0,0,qprod(my_quat,th1_l(vtl(b)^-1))]])) */

list_of_good_quats=[];for(i=1,length(bigd1list_l),if(test_at_2_general(bigd1list_l[i],3),list_of_good_quats=concat(list_of_good_quats,[bigd1list_l [i]]),))

fl_l=[];for(b=0,l,for(i=1,length(list_of_good_quats),mess=th_l(list_of_good_quats[i])*vtl(b)^-1;if(testU0_l(mess),;fl_l=concat(fl_l,[[0,0,0,0,th_l(list_of_good_quats[i])]]),)))

/* these functions create the matrix of T_l */

getNewEpsilonMatrices_l(a,b)=temp=[]; \
for(i=1,length(fl_l),if(fl_l[i][1]==a&&fl_l[i][4]==b,temp=concat(temp,[fl_l[i][5]])));temp

h_ij_new_l(i,j,k,acc)=subst(genfn_new(getNewEpsilonMatrices_l(i,j),k),x0_l,sqrt(XX_l+O(l^acc)))

getRow_new_l(acc,k,dim,row)=tempm=(my_hm_new2(h_ij_new_l(row,0,k,acc),dim));\
for(i=1,length(clist)-1,tempm=concat(tempm,my_hm_new(h_ij_new_l(row,i,k,acc),dim)));\
mattranspose(tempm)

Tlmatrix(acc,k,dim)=tempp=getRow_new_l(acc,k,dim,length(clist)-1);\
forstep(i=length(clist)-2,0,-1,print(i);tempp=concat(getRow_new_l(acc,k,dim,i),tempp));\
mattranspose(tempp)

charpowerseries_l(acc,k,dim)=matdet(Tlmatrix(acc,k,dim)-BIGT*matid(dim*length(clist)))

findSlopes_l(acc,k,dim)=newtonpoly(charpowerseries_l(acc,k,dim),p)

charpowerseries_l(50,6,5)