/* Code written by Dan Jacobs, (c) 2002.
Copied out by LJPK, 2004. 
Extensions written by LJPK, 2004.
*/

q=5;
p=factor(q)[1,1];

l=3;

flag=0; /* this is 1 for classic, 0 for extra level structure */
flag_small=0; /* 0 for everything, 1 just for the squares */

verbose_flag=1; /* 0 for less verbose, 1 for everything */

print("******************************")
print("Program for calculating the action of U_p on quaternionic modular 
forms")
print("Original written by Dan Jacobs, (c) 2002")
print("Extended by L J P Kilford, (c) 2004-2005")
print("******************************")

print();
print("Version 0.1.2");
print();

/* version 0.1.2 because we now have T_l working! At least for standard 
level structure at 2 anyway */

if(verbose_flag,print("Verbose mode on");print())

print("The prime p is ",p)

if(flag,print("with standard level structure at 2"),\
print("with extra level structure at 2"))

if(flag_small,print("of the smaller level U_{sq}(",q,")"),\
print("of level U_1(",q,")"))
print()

/* This allocates the values of x0 and y0 for the homomorphism between 
D_q and M_2(\Q_q) */

if(p%4==1,XX=-1;y0=0);

if(p%8==3,XX=-2;y0=1);

if(p%8==7,XX=-5;y0=2); /* NB this may not work */

ACC=40;

if(verbose_flag,print("The p-adic accuracy is ",ACC),)

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

th(v)=[v[1]+v[2]*x0+v[4]*y0, v[2]*y0-v[3]-v[4]*x0; v[2]*y0+v[3]-v[4]*x0, 
v[1]-v[2]*x0-v[4]*y0]

th1(mm) = [(mm[1,1]+mm[2,2])/2, ((mm[2,2]-mm[1,1])*x0 - 
(mm[1,2]+mm[2,1])*y0)/2, 
(mm[2,1]-mm[1,2])/2,((mm[2,2]-mm[1,1])*y0+(mm[1,2]+mm[2,1])*x0)/2]

/* Now we use the new coset functions from new_quat_stuff */

\r auxiliary_quaternion_algebra.g

unitlist_H = vector(length(unitlist_H), i, th(unitlist_H[i]));

/* Choose whether we want small or big space of forms */

if(flag_small==0,clist=findCosetMatrices(),clist=findCosetMatrices_new());

print("The cosets are ",clist)
print(" ")
/* Inverting the quaternion [q0,q1,q2,q3] */

quatinv(v) = vv = v; N = sum(ii=1,4,v[ii]^2); for(ii=2,4,vv[ii] *= -1/N); 
vv[1] /= N; vv

/* Finding elements d1 of O_D with N(d1)=p and d1 == [0, *; 0, *] mod p 
*/

d1list=[];

bound=floor(sqrt(p)+3);

forstep(d0=-bound, bound, 1/2, forstep(d1=-bound, bound, 1/2, \
forstep(d2=-bound, bound, 1/2, forstep(d3=-bound, bound, 1/2, \
if(d0^2+d1^2+d2^2+d3^2==p, \
mm = subst(th([d0,d1,d2,d3]), x0, sqrt(XX+O(p^ACC))); \
if((mm[1,1]%p == 0) && (mm[2,1]%p == 0), \
d1list=concat(d1list, [[d0,d1,d2,d3]])))))))

bigd1list=[];

forstep(d0=-bound, bound, 1/2, forstep(d1=-bound, bound, 1/2, \
forstep(d2=-bound, bound, 1/2, forstep(d3=-bound, bound, 1/2, \
if(d0^2+d1^2+d2^2+d3^2==p, \
bigd1list=concat(bigd1list, [[d0,d1,d2,d3]]))))))


/* forstep(d0=-2, 2, 1/2, forstep(d1=-2, 2, 1/2, \
forstep(d2=-2, 2, 1/2, forstep(d3=-2, 2, 1/2, \
if(d0^2+d1^2+d2^2+d3^2==3, print([d0,d1,d2,d3])))))) */

/* Calculating the units of O_D and their matrix versions */

unitlist=[];

forstep(a=-2, 2, 1/2, forstep(b=-2, 2, 1/2, \
forstep(cc=-2, 2, 1/2, forstep(d=-2, 2, 1/2, \
if(a^2+b^2+cc^2+d^2==1, \
unitlist = concat(unitlist, [[a,b,cc,d]]))))))

smallunitlist=[];

forstep(a=-2, 2, 1, forstep(b=-2, 2, 1, \
forstep(cc=-2, 2, 1, forstep(d=-2, 2, 1, \
if(a^2+b^2+cc^2+d^2==1, \
smallunitlist = concat(smallunitlist, [[a,b,cc,d]]))))))

/*forstep(a=-2, 2, 1/2, forstep(b=-2, 2, 1/2, \
forstep(cc=-2, 2, 1/2, forstep(d=-2, 2, 1/2, \
if(a^2+b^2+cc^2+d^2==1, \
print([a,b,cc,d])))))) */

matunitlist = vector(length(unitlist), i, th(unitlist[i]));

smallmatunitlist = vector(length(smallunitlist), i, 
th(smallunitlist[i]));

/* "test": matrices for integral entries */

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

testU0(n)=m=subst(n,x0,sqrt(XX+O(p^ACC)));test(m)

testU1(n)=m=subst(n,x0,sqrt(XX+O(p^ACC)));test(m)&&((m[2,1]%q)==0)&&((m[2,2]%q)==1)

testUsq(n)=m=subst(n,x0,sqrt(XX+O(p^ACC)));test(m)&&((m[2,1]%q)==0)&&issquare(m[2,2]+O(p^2))

/* Step 1 Factorisation; write c_i v_t^{-1} = du */

fact1list = [];

vtp(t) = [p,0; q*t, 1]

myc(iii) = clist[iii+1]

fact1list = [];

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

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

for(t=0,p-1, for(i=0, length(clist)-1, for(j=1, length(bigd1list),\
u = th(bigd1list[j])*myc(i)*(1/vtp(t)); \
uu = subst(u, x0, sqrt(XX + O(p^ACC))); \
if(test(uu), fact1list = concat(fact1list, \
[[t,i,j,[quatinv(bigd1list[j]), u, vtp(t)]]]); break(1)))))

/* Testing the factorisations */

/*print("Testing the factorisations")

for(i=1, length(fact1list), mmmm=fact1list[i]; \
print1(Mod(myc(mmmm[2])*(1/vtp(mmmm[1]))-\
th(mmmm[4][1])*mmmm[4][2], x0^2+XX) "; "))*/

print("Testing the factorisations ...")

testflag=1;

for(i=1, length(fact1list), mmmm=fact1list[i]; \
testflag=testflag*(0==(myc(mmmm[2])*(1/vtp(mmmm[1]))-subst(th(mmmm[4][1])*mmmm[4][2], 
\
x0,sqrt(XX+O(p^ACC)))))\
)

if(testflag==1,print("The factorisations are OK"),error("The 
factorisations are not OK"))

/*for(i=1, length(fact1list), mmmm=fact1list[i]; \
print(myc(mmmm[2])*(1/vtp(mmmm[1]))-subst(th(mmmm[4][1])*mmmm[4][2], \
x0,sqrt(XX$
)*/


/*print(subst(th(mmmm[4][1])*mmmm[4][2], x0,sqrt(XX+p^ACC)))))*/


/* Finding alpha & c_i st c_i^{-1} alpha^{-1} u \in U_1(q) */

/*D was [-1/2, 1/6, 1/6, -1/6] in the original case  --- let's do this in 
general*/

D = fact1list[1][4][1];

fl=[];

if(verbose_flag,print("Functions for testing matrices"),)

mysearch(mlist,myx)=mflag=0;for(i=1,length(mlist),if(mlist[i]==myx,mflag=1));mflag

/* changed to matunitlist for q=9 */

/*getFL(start)=output=[];for(aa=0,length(clist)-1,for(bb=1,length(matunitlist),\
vp=(1/myc(aa))*(1/matunitlist[bb])*start;\
vvpp=subst(vp,x0,sqrt(XX+O(p^ACC))); vvpp %=q; \
if(testU1(vvpp),output=([matunitlist[bb],myc(aa),vp]))));output*/

testUstar(vvpp)=temp=0;if(flag_small==0,temp=testU1(vvpp),temp=testUsq(vvpp));temp

getFL(start)=output=[];for(aa=0,length(clist)-1,for(bb=1,length(unitlist_H),\
vp=(1/myc(aa))*(1/unitlist_H[bb])*start;\
vvpp=subst(vp,x0,sqrt(XX+O(p^ACC))); vvpp %=q; \
if(testUstar(vvpp),output=([unitlist_H[bb],myc(aa),vp]))));output

getFLtest(start)=output=[];for(aa=0,length(clist)-1,for(bb=1,length(matunitlist),\
vp=(1/myc(aa))*(1/matunitlist[bb])*start;\
vvpp=subst(vp,x0,sqrt(XX+O(p^ACC))); vvpp %=q; \
if(testU1(vvpp),output=concat(output,[[matunitlist[bb],myc(aa),vp]]))));output


if(verbose_flag,print("Further factorisations"),)

/*for(i=1,length(fact1list), up =fact1list[i][4][2]; for(ii=0, 
length(clist)-1, \
for(jj=1, length(matunitlist), \
vp=(1/myc(ii))*(1/matunitlist[jj])*up; \
vvpp=subst(vp,x0,sqrt(XX+O(p^ACC))); vvpp %=p; \
if((vvpp[2,1] ==0) && (vvpp[2,2] == 1), \
innerthing=subst(th(D)*matunitlist[jj]*myc(ii)*vp,x0,sqrt(XX+O(p^ACC)))*fact1list[i][4][3];\
bigflag=mysearch(clist,innerthing);\
if(bigflag,\
fl=concat(fl, [[D, matunitlist[jj], myc(ii), vp, 
fact1list[i][4][3]]]);print([i,ii,jj]);break(2))))))*/

for(i=1,length(fact1list),temp=getFL(fact1list[i][4][2]);\
if(length(temp)!=0,fl=concat(fl,[concat([fact1list[i][4][1],temp[1]],[temp[2],temp[3],fact1list[i][4][3]])])))

/*print("debug 2")*/

/* Output of the decomposition test */

/* for(i=1, length(f1), print(f1[i])) */

if(verbose_flag,print("Finding the matrices for the generating 
functions"),)

/* finds the matrices for the generating function */

getNewEpsilonMatrices(a,b)=temp=[]; \
for(i=1,p,if(fl[length(clist)*(i-1)+a+1][3]==myc(b),temp=concat(temp,[fl[length(clist)*(i-1)+a+1][4]*fl[length(clist)*(i-1)+a+1][5]])));temp

/* generating function for the components of the U_p operator */

genfn(matrix1,k)=temp=0; \
for(i=1,length(matrix1),temp=temp+(matrix1[i][2,1]*x+matrix1[i][2,2])^k/((matrix1[i][2,1]*x+matrix1[i][2,2])*(matrix1[i][2,1]*x+matrix1[i][2,2]-matrix1[i][1,1]*x*y-matrix1[i][1,2]*y)));temp

genfn_new(matrix1,k)=temp=[]; \
for(i=1,length(matrix1),temp=concat(temp, \
[(matrix1[i][2,1]*x+matrix1[i][2,2])^k/((matrix1[i][2,1]*x+matrix1[i][2,2])*(matrix1[i][2,1]*x+matrix1[i][2,2]-matrix1[i][1,1]*x*y-matrix1[i][1,2]*y))]));temp

genfnXX(matrix1,k)=temp=0;for(i=1,length(matrix1),temp=temp+(matrix1[i][2,1]*A+matrix1[i][2,2])^k/((matrix1[i][2,1]*A+matrix1[i][2,2])*(matrix1[i][2,1]*A+matrix1[i][2,2]-matrix1[i][1,1]*A*B-matrix1[i][1,2]*B)));temp

if(verbose_flag,print("Defining the U_p operator"),);



h_ij_new(i,j,k,acc)=subst(genfn_new(getNewEpsilonMatrices(i,j),k),x0,sqrt(XX+O(p^acc)))

h_ij_real(i,j,k,acc)=subst(genfn_new(getNewEpsilonMatrices(i,j),k),x0,sqrt(XX))

matrix_coeff(mess,i,j)=tempmess=mess+O(y^(j+1))+O(x^(i+1));polcoeff(polcoeff(tempmess,j,y),i,x)

my_hm_new2(hij,dim)=temp=matrix(dim,dim,i,j,0);for(k=1,length(hij),temp=temp+matrix(dim,dim,i,j,matrix_coeff(hij[k],i-1,j-1)));temp

my_hm_new(hij,dim)=temp=matrix(dim,dim,i,j,0);for(k=1,length(hij),temp=temp+matrix(dim,dim,i,j,polcoeff(polcoeff(hij[k]+O(x^(dim+1)),i-1,x)+O(y^(dim+1)),j-1,y)));temp

my_hm_real(hij,dim)=temp=matrix(dim,dim,i,j,0);for(k=1,length(hij),temp=temp+matrix(dim,dim,i,j,polcoeff(polcoeff(hij[k]+O(x^(dim+1)),i-1,x)+O(y^(dim+1)),j-1,y)));temp


/* creating the matrix for the U_p operator */

/* first we make the rows, because it was too hairy to do it all at once 
*/

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

getRow_real(acc,k,dim,row)=tempm=(my_hm_real(h_ij_real(row,0,k,acc),dim));\
for(i=1,length(clist)-1,tempm=concat(tempm,my_hm_real(h_ij_real(row,i,k,acc),dim)));\
mattranspose(tempm)

/* then we put them together to make the whole matrix */

/* bigUp_new(acc,k,dim)=tempp=getRow_new(acc,k,dim,length(clist)-1);\
forstep(i=length(clist)-2,0,-1,print(i);tempp=concat(getRow_new(acc,k,dim,i),tempp));\
mattranspose(tempp)

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

findSlopes(acc,k,dim)=newtonpoly(matdet(bigUp_new(acc,k,dim)-BIGT*matid(dim*length(clist))),p)

charpowerseries(acc,k,dim)=matdet(bigUp_new(acc,k,dim)-BIGT*matid(dim*length(clist)))

charpowerseries_real(acc,k,dim)=matdet(bigUp_real(acc,k,dim)-BIGT*matid(dim*length(clist)))*/

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

charpowerseries(acc,k,dim)=matdet(Upmatrix(acc,k,dim)-BIGT*matid(dim*length(clist)))

findSlopes(acc,k,dim)=newtonpoly(charpowerseries(acc,k,dim),p)

/* this loads the new code for T_l */

\r functions_to_compute_T_l.g

findSlopes(100,3,2)