/*
 
 10th October, 2014
 
 Copyright (C) 2012 Alan Lauder, 2014 Jan Vonk
 
 Alan Lauder : original implementation of main functions for primes p >= 5
 Jan Vonk    : optimizations and extension to primes p = 2,3. Incorporation of recent speed-ups by Alan Lauder and Steve Donnelly.
 
 ***
 
 /*

----------------------------------------------------------
--|                  The U_p operator                  |--
----------------------------------------------------------

The main function is 

HeckeSeries(p,N,kseq,m: weightbound := 6)

Input: 

p: a prime
N: a positive integer - the level will be Gamma0(N)
kseq: a list [k0,k1, ..., km] of positive integers, all congruenct mod p-1 if p > 3, mod 4 if p = 2 and mod 6 if p = 3.
m: a positive integer
weightbound: (optional) an auxiliary parameter in the routine that computes complimentary subspaces

Output:

The reverse characteristic power series modulo p^m of the U_p operator on the space of p-adic overconvergent modular forms of level N and weight k, for each k in kseq.


Examples:

 > HeckeSeries(5,1,1728,30);
 372529029846191406250*t^6 - 121953464508056640625*t^5 + 448701049074394531250*t^4 + 197944154954912503125*t^3 - 66796265305702220600*t^2 + 100898070553738936224*t + 1
 
 > HeckeSeries(2,5,120,20); 
 q-adic precision NN =  66
 Finding bases for modular forms of low weight.
 Weight k =  2
 Weight k =  4
 Weight k =  6
 Constructing complementary spaces modulo p
 ell, numberoftries:  33 315
 Constructing complementary spaces ...
 385024*t^8 + 158720*t^7 - 96128*t^6 + 360528*t^5 - 478544*t^4 + 302002*t^3 + 416973*t^2 + 1
 
Note: There are several lines of output while the algorithm is running to give more information on what is happening internally. 

*/

Attach("./Documents/Wiskunde/Computations/Weight1.m"); // Steve's new weight 1 code.


// ----------------------------------------------------------
// --|                  The U_p operator                  |--
// ----------------------------------------------------------

//  *** AUXILIARY CODE FOR GENERAL LEVEL ***

// CLASSICAL MODULAR FORMS

// Returns Eisenstein series E_k over given ring S modulo q^NN normalised
// to have constant term a_0 = 1.

function ESeries(k,NN,S)

	R<q> := PowerSeriesRing(S,NN);
	a1 := -S!(2*k/Bernoulli(k));
	Ek := 1 + a1*q;
	for n := 2 to NN-1 do
		coeffn := 0;
		divs := Divisors(n);
		for d in divs do
			coeffn := coeffn + d^(k-1);
		end for;
		Ek := Ek + a1*coeffn*q^n;
	end for;

	return Ek;

end function;

// Returns the power of the Hasse invariant that has been lifted.
// This is always 1, except when p = 2,3.

function ap(p)

	a := p;
	if p eq 2 then
		a := 4;
	end if;
	if p eq 3 then 
		a := 3;
	end if;
	return a;

end function;

// Returns the weight of the Eisenstein series to be used for lifting the Hasse invariant. 
// This is always p-1, unless p = 2,3.

function EisWeight(p)

	return ap(p)*(p-1);

end function;

// Weightbound function default.

function WtB(chi)

	weightbound:=6;
	if Type(chi) ne RngIntElt and (chi(-1) eq -1) then
		weightbound := 4;
	end if;
	
	return weightbound;
	
end function;

// DEGREE BOUND AND CHARACTERISTIC SERIES

// Returns dimension of space of classical modular forms for Gamma_1(N,chi) of weight k.	

function DimensionMNk(chi,k)

	if k gt 0 then
		M:=ModularForms(chi,k); 
		return Dimension(M);
	elif ((k eq 0) and chi eq Parent(chi)!1) then // trivial character
		return 1;
	else // k < 0 or k = 0 with non-trivial character
		return 0;
	end if;

end function;

// Given square matrix A over Z/(p^m) returns det(1 - A t).
// NOTE: A should be defined over IntegerRing(p^m) rather
// than pAdicRing(p,m), otherwise this will take forever.

function ReverseCharacteristicPolynomial(A)
	
	P:=CharacteristicPolynomial(A);
	return ReciprocalPolynomial(P);

end function;

// Check kseq contains weights congruent modulo p-1.

function IsValidWeightSequence(kseq,p)

	len:=#kseq;
	if len eq 0 then
		return false;
	end if;
	
	result:=true;
	k0:=kseq[1] mod EisWeight(p);
	for i:=2 to len do
		if (kseq[i] mod EisWeight(p)) ne k0 then
			result:=false;
		end if;
	end for;

	return result;

end function;

// *** LEVEL N > 1 ***

// CLASSICAL MODULAR FORMS AND LINEAR ALGEBRA

// Converts a "modular form" to power series in q over Zp/(p^mdash).

function qexpToPS(f,p,mdash,NN)

	Zp:=IntegerRing(p^mdash);
	R<q>:=PowerSeriesRing(Zp,NN);
	
	fPS:=R!0;
	for i:=0 to NN-1 do
		fPS:=fPS + Coefficient(f,i)*q^i;
	end for;
	
	return fPS;

end function;

// Constructs bases of power series expansions for modular forms of level
// Gamma_0(N) and (even) weight k <= weightbound.
// The basis elements are in Zp[[q]]/(p^m,q^NN).
// MAJOR SPEED-UP FROM TRIVIAL CHANGE - SEE COMMENT IN Rankin-1/Tests.m.

function LowWeightBases(chi,p,mdash,NN,weightbound)

	print "q-adic precision NN = ", NN;
	print "Finding bases for modular forms of low weight.";
	generators:=[];
	
	if chi(-1) eq -1 then
		k0 := 1;
	else 
		k0 := 2;
	end if;
	
	for k:=k0 to weightbound by Integers()!(1+(1+chi(-1))/2) do
		print "Weight k = ", k;
		if IsEven(k) and chi(-1) eq -1 then
			M := ModularForms(Modulus(chi),k);
		else
			M := ModularForms(chi,k);
		end if;
		// SetPrecision(M,NN);
		basisweightk:=[];
		b:=Basis(M,NN);
		/* for i:=1 to Dimension(M) do
			f:=M.i;
			Append(~basisweightk,qexpToPS(f,p,mdash));
		end for;
		*/
		for f in b do
			Append(~basisweightk,qexpToPS(f,p,mdash,NN));
		end for;
		Append(~generators,basisweightk);
	end for;

	return generators;

end function;

// Returns "random" bases for spaces of modular forms of level Gamma_0(N) and
// (even) weight k <= weightbound. The basis elements are in Zp[[q]]/(p^mdash,q^NN).
// The random bases are chosen to be ones for modular forms with coefficients in Zp.

function RandomLowWeightBases(N,p,mdash,NN,weightbound)

	LWB:=LowWeightBases(N,p,mdash,NN,weightbound);
	R:=Parent(LWB[1][1]);
	
	RandomLWB:=[];

	for i:=1 to #LWB do // weight k = 2*i
		basisweightk:=LWB[i];
		dimweightk:=#basisweightk;
		coeffsmat:=Random(GL(dimweightk,p));
		randombasisweightk:=[];
		for j:=1 to dimweightk do
			coeffsmatj:=coeffsmat[j];
			f:=R!0;
			for l:=1 to dimweightk do
				f:=f + (Integers()!coeffsmatj[l])*basisweightk[l];
			end for;
			Append(~randombasisweightk,f);
		end for;
		Append(~RandomLWB,randombasisweightk);
	end for;

	return RandomLWB;

end function;



// Computes "random" solution in non-negative integers to
// equation 2 a_2 + 4 a_4 + ... + b a_b = K, where K is even and b is the largest even number <= B.
// It uses a greedy algorithm and the output is not uniformly
// distributed.

function RandomSolutionEven(B,K)

	a:=[];
	Sum := 0;
	assert IsEven(K);
	if IsOdd(B) then
		b := B-1;
	else
		b := B;
	end if;
	
	for i:=b to 4 by -2 do
		ai:=Random(0,Floor(K/i));
		Append(~a,ai);
		Sum := Sum + ai;
		K:=K - ai*i;
	end for;
	a1:=Integers()!(K/2);
	Append(~a,a1);
	arev:=Reverse(a);
	
	return arev, Sum + a1;
	
end function;

// Auxiliary function, adds two lists of integers together.

function AddList(L1,L2)

	L := [];
	for i:= 1 to #L1 do
		Append(~L,L1[i]+L2[i]);
	end for;
	
	return L;

end function;

// Computes "random" solution in non-negative integers to
// equation 2 a_2 + 3 a_3 + ... + B a_B = K.
// It uses a greedy algorithm and the output is not uniformly
// distributed.

function RandomSolutionOdd(B,K)

	a:=[];
	
	for i:=B to 2 by -1 do
		ai:=Random(0,Floor(K/i));
		Append(~a,ai);
		K:=K - ai*i;
	end for;
	a1:=K;
	Append(~a,a1);
	arev:=Reverse(a);
	
	Sum := 0;
	for i:=1 to Ceiling(#arev/2) do
		Sum := Sum + arev[2*i-1];
	end for;
	
	return arev, Sum;
	
end function;


// Computes "random" combination of previously found weights that sum up to something strictly less than K.
// This should speed up the function RandomNewBasisModp.

function RandomSolutionComp(Comp,K,p,weightbound)

	a:=[];
	Sum := 0;
	l := #Comp;
	I := K;
	for i:= 1 to l do
		if (I-i*EisWeight(p)) gt 20 then
			ai := Random(0,Floor(K/(I-i*EisWeight(p))));
			Append(~a,ai);
			Sum := Sum + ai;
			K := K - ai*(I-i*EisWeight(p));
		else 
			Append(~a,0);
		end if;
	end for;
	arev:=Reverse(a);

	return arev, K, Sum;

end function;


// Computes the code of a random element of weight k, combining
// the functions RandomSolutionOdd/Even and RandomSolutionComp above. 

function RandomSolution(p,chi,k,weightbound,CompSpacesCode,LWBModp)

	RandomCode := [];
	par := 0;
	
	while par eq 0 do
		CompExps,K,sum1:=RandomSolutionComp(CompSpacesCode,k,p,weightbound);
		if IsOdd(k) then
			exps, sum2:=RandomSolutionOdd(weightbound,K);
		else
			exps, sum2:=RandomSolutionEven(weightbound,K);
		end if;
		if (((sum1 + sum2) mod Order(chi)) eq 1) or (Conductor(chi) eq 1) then // We test to see whether the character of the product is chi.
			par := 1;
		end if;	
	end while;
	for j := 1 to #exps do
		for l:=1 to exps[j] do
			Append(~RandomCode,[j,Random(1,#LWBModp[j])]);
		end for;
	end for;
	for j := 1 to #CompExps do
		for l:=1 to CompExps[j] do 
			RandomCode := RandomCode cat CompSpacesCode[j][Random(1,#CompSpacesCode[j])];
		end for;
	end for;
	
	return RandomCode;
	
end function;

// Returns position of the leading entry of a NON-ZERO vector.

function LeadingPosition(v)

	d:=Dimension(Parent(v));	
	i:=1;
	vi:=v[i];
	while vi eq 0 do
		i:=i+1;
		vi:=v[i];
	end while;
	
	return i;

end function;

// Returns leading entry of a vector.

function LeadingEntry(v)

	return v[LeadingPosition(v)];

end function;

// COMPLEMENTARY SPACES.

// Auxiliary function used in main loop of ComplementarySpacesModp

function RandomNewBasisModp(chi,p,k,LWBModp,OldBasisModp,weightbound,CompSpacesCode)
	
	R:=Parent(LWBModp[1][1]); // power series ring mod p,q^elldash.
	
	// Construct TotalBasisModp
	TotalBasisModp:=OldBasisModp; // Recall E_(p-1) = 1 mod p.
	elldash:=Degree(TotalBasisModp); 
	
	// Case k0 + i(p-1) = 0 + 0(p-1) = 0
	if k eq 0 then 
		// empty string <-> 1 for NewBasis.
		TotalBasisModp:=sub< V | V.1 > where V is VectorSpace(GF(p), Degree(TotalBasisModp)); // (1,0,0,...,0)
		return TotalBasisModp, [[]], 0; // March 28: added zero for numberoftries
	end if;
	
	// Case k = k0 + i(p-1) > 0
	di:=DimensionMNk(chi,k);
	diminus1:=DimensionMNk(chi,k-EisWeight(p));
	mi:=di - diminus1;
	
	// Generate mi new forms in weight k.
	NewBasisCode:=[];
	rk:=diminus1; // number of non-zero rows in TotalBasisModp
    numberoftries:=0;
	for i:=1 to mi do // extra forms 
		while (rk lt diminus1 + i) do
			// construct random form of weight k
            numberoftries:=numberoftries+1;
			RandomCode := RandomSolution(p,chi,k,weightbound,CompSpacesCode,LWBModp); // Generate random code in low weight basis.
			TotalBasisi:=R!1;
			TotalBasisiCode:=[]; // empty string <-> 1.
			for j:=1 to #RandomCode do
				TotalBasisi:=TotalBasisi*LWBModp[RandomCode[j][1]][RandomCode[j][2]];
			end for;
			TotalBasisiCode := TotalBasisiCode cat RandomCode;
			Include(~TotalBasisModp, Vector([Coefficient(TotalBasisi,j): j in [0 .. elldash-1]]),~new);
			rk:=Dimension(TotalBasisModp);
		end while;
		Append(~NewBasisCode,TotalBasisiCode); // this one increased dimension.
	end for;
			
	return TotalBasisModp,NewBasisCode,numberoftries;

end function;

// Finds complementary spaces modulo p and returns a list of "codes" describing
// what products of basis forms were chosen.

function ComplementarySpacesModp(chi,p,k0,n,elldash,LWBModp,weightbound)

	print "Constructing complementary spaces modulo p";
	CompSpacesCode:=[];
	ell:=DimensionMNk(chi,k0 + n*EisWeight(p));
	OldBasisModp:=sub<VectorSpace(GF(p),elldash) | >; 
    totalnumberoftries:=0;
	for i:=0 to n do
		TotalBasisModp,NewBasisCodemi,numberoftriesi:=RandomNewBasisModp(chi,p,k0 + i*EisWeight(p),LWBModp,OldBasisModp,weightbound,CompSpacesCode);
		Append(~CompSpacesCode,NewBasisCodemi);
		OldBasisModp:=TotalBasisModp;
        totalnumberoftries+:=numberoftriesi;
	end for;
    print "ell, numberoftries: ", ell, totalnumberoftries;

	return CompSpacesCode;

end function;

// Reduces the basis of low weight forms mod (p,q^elldash).

function LowWeightBasesModp(LWB,p,elldash)

	R:=PowerSeriesRing(GF(p),elldash);
	LWBModp:=[];	
	
	for i:=1 to #LWB do // weight k = 2*i
		LWBModpWeightk:=[];
		for j:=1 to #LWB[i] do
			Append(~LWBModpWeightk,R!LWB[i][j]);
		end for;
		Append(~LWBModp,LWBModpWeightk);
	end for;
	
	return LWBModp;

end function;

// Returns complementary spaces W = [W_0,W_1,...,W_n] as list of 
// basis elements W_i modulo (p^mdash,q^elldashp).

function ComplementarySpaces(chi,p,k0,n,mdash,elldash,elldashp,weightbound)

	// Find q-expansions for k <= weightbound mod (p^mdash,q^elldashp)

	LWB:=RandomLowWeightBases(chi,p,mdash,elldashp,weightbound);
	
	// Find required products mod (p,q^elldash): this was suggested to the
	// author by John Voight.

	LWBModp:=LowWeightBasesModp(LWB,p,elldash);
	CompSpacesCode:=ComplementarySpacesModp(chi,p,k0,n,elldash,LWBModp,weightbound);
    // print "N,p,k0,n,elldash", N,p,k0,n,elldash;
        
	// Reconstruct products mod (p^mdash,q^elldashp)
	W:=[];
	
	// Epm1:=ESeries(p-1,elldashp,IntegerRing(p^mdash)); March 21, 2012
    Zpmdashq:=PowerSeriesRing(IntegerRing(p^mdash),elldashp);

	OldBasis:=[];
	for i:=0 to n do
		CompSpacesCodemi:=CompSpacesCode[i+1];
		Wi:=[];
		for k:=1 to #CompSpacesCodemi do
			// "code" <-> [ ... (k/2,a) ...] where k = weight, a = index of element
			// of weight k chosen, and then one takes the product over the list.
			CompSpacesCodemik:=CompSpacesCodemi[k];			
            Wik:=Zpmdashq!1;
			for j:=1 to #CompSpacesCodemik do 
				wt:=CompSpacesCodemik[j][1];
				ind:=CompSpacesCodemik[j][2];
				Wik:=Wik*LWB[wt,ind];
			end for;
			Append(~Wi,Wik);
		end for;
		Append(~W,Wi); 
	end for;
	
	return W;

end function;


// AUXILIARY CODE: KATZ EXPANSIONS

// Returns basis of Katz Expansions modulo (p^mdash,q^elldashp) as a matrix
// of coefficients. These are (part of) a basis for the space of overconvergent
// modular forms of level N and weight k0.

function HigherLevelKatzExp(chi,p,k0,m,mdash,n,elldash,elldashp,weightbound)

	ordr := ap(p)/(p+1); 
	S := IntegerRing(p^mdash);
	Ep1 := ESeries(EisWeight(p),elldashp,S);
	R := Parent(Ep1);
	q := R.1;

	// construct basis as list of series

	time W := ComplementarySpaces(chi,p,k0,n,mdash,elldash,elldashp,weightbound);
	print "...Complementary spaces.";

	KatzBasis := [];
	Epinv := Ep1^-1;
	Ep1minusj := R!1;
	for j := 0 to n do
		Wj := W[j+1];
		dimj := #Wj;
		for i := 1 to dimj do
			wji := Wj[i];
			b := p^Floor(ordr*j)*wji*Ep1minusj;
			Append(~KatzBasis,b);
		end for;
		Ep1minusj := Ep1minusj*Epinv;
	end for;
			
	// extract basis as matrix
	
	ell:=#KatzBasis;
	M:=ZeroMatrix(S,ell,elldashp);
	for i:=1 to ell do
		for j:=1 to elldashp do
			M[i,j]:=Coefficient(KatzBasis[i],j-1);
		end for;
	end for;
		
	// reduce: already in "almost" echelon form
	
	Mred:=EchelonForm(M);

	// Mred:=MyEchelonFormB(M,p,mdash); // Dec 20: load CCpMF-magma/echelonform.m

	return Mred,Ep1; 

end function;

// Returns ldash, the Sturm bound.

function Computeldash(p,chi,k0,m)

	n:=Floor(((p+1)/EisWeight(p))*(m+1+(ap(p)-1)/(p+1)));
	ldash:=PrecisionBound(ModularForms(chi,k0 + n*EisWeight(p)));

	return ldash;

end function;


// MAIN FUNCTION.

// Returns matrix A modulo p^m in Step 5 of Algorithm 2 for each k in kseq.
// kseq has already been checked to have all weights congruent modulo (p-1).

function HigherLevelUpGj(p,chi,kseq,m,weightbound)// Compute weight bound

	// Step 1
	
	k0:=kseq[1] mod EisWeight(p);
	
	elldash:=Computeldash(p,chi,k0,m);
	elldashp:=elldash*p;
	n:=Floor(((p+1)/EisWeight(p))*(m+1+(ap(p)-1)/(p+1)));	 
	mdash:=m + Ceiling(n*ap(p)/(p+1));
	
	// Steps 2 and 3

	e,Ep1:=HigherLevelKatzExp(chi,p,k0,m,mdash,n,elldash,elldashp,weightbound); 
	ell:=Dimension(Parent(Transpose(e)[1]));
	S:=Parent(e[1,1]);

	// Step 4

	q := Parent(Ep1).1; 
	Ep1p := Evaluate(Ep1,q^p);
	Ep1pinv := Ep1p^-1;
	G := Ep1*Ep1pinv;
	Aseq:=[];
	
	for k in kseq do
		kdiv:=k div EisWeight(p);
		if kdiv ne 0 then
			Gkdiv := G^kdiv;
			u:=Parent(e)!0; // ell x elldashp zero matrix.
		
			for i:=1 to ell do
				ei:=Parent(Gkdiv)!0;
				for j:=1 to elldashp do // extract rows of a as series
					ei:=ei + e[i,j]*q^(j-1);
				end for;
				Gkdivei:=Gkdiv*ei; // act by G^kdiv
				for j:=1 to elldashp do // put back as rows of ak
					u[i,j]:=Coefficient(Gkdivei,j-1);
				end for;
			end for;
		else
			u := e;
		end if;
		// Step 5

		T:=ZeroMatrix(S,ell,elldash);
	
		for i:=1 to ell do
			for j:=1 to elldash do
				T[i,j]:=u[i,p*(j-1) + 1];
			end for;
		end for;
	
		// Step 6: Solve T = AE for A.

		A := ZeroMatrix(S,ell,ell);
		
		E := Parent(T)![[e[j,l]: l in [1 .. elldash]]: j in [1 .. ell]]; // March 25, 2012
    	t0:=Cputime();
    	E_leadingpos:=[Depth(E[j]): j in [1 .. ell]];
	
		for i := 1 to ell do
			Ti := T[i];		
			for j := 1 to ell do
				ej := E[j];
				ejleadpos:=E_leadingpos[j];
				lj:=Integers()!ej[ejleadpos];
				A[i,j] := S!((Integers()!Ti[ejleadpos])/lj);
				Ti := Ti - A[i,j]*ej;
			end for;
		end for;
		
		// Append(~Aseq,MatrixRing(pAdicRing(p,m),ell)!A);
		
		Append(~Aseq,MatrixRing(IntegerRing(p^m),ell)!A);
		
	end for;
		
	return Aseq,e,mdash;


end function;

//  *** LEVEL 1 ***

// Returns dimension of space of modular forms of weight k and level 1.

function DimensionMk(k)

	if ((k mod 2) eq 1) or (k lt 0) then
		return 0;
	end if;
	kmod12 := k mod 12;
	if (kmod12 eq 2) then
		return Floor(k/12);
	else
		return Floor(k/12) + 1;
	end if;

end function;

// Returns basis for W_i (those elements in the Miller basis of weight k which 
// are not in E_(p-1) x {weight (k - (p-1)) forms} and a power of delta
// which can be reused on the next call of function.

function ComputeWi(k,p,delta,deltaj,E4,E6)
	
	// Define a and b
	a := k mod 3;
	b := (k div 2) mod 2;
		
	// Compute dimensions required for Miller basis
	d := DimensionMk(k) - 1;
	e := DimensionMk(k - EisWeight(p)) - 1;
	
	// Construct basis for Wi
	Wi := [];
	for j := e+1 to d do
		// compute aj = delta^j*E6^(2*(d-j) + b)*E4^a
		aj := deltaj*E6^(2*(d-j) + b)*E4^a;
		deltaj := deltaj*delta;
		Wi := Append(Wi,aj);
	end for;

	return Wi,deltaj;

end function;

// AUXILIARY CODE: KATZ EXPANSIONS

// Returns list e of Katz expansions e_(i,s) from Step 3 of Algorithm 1
// and E_(p-1) for use in Step 4.

function KatzExpansions(k0,p,ellp,mdash,n)
	
	S := IntegerRing(p^mdash); // See Footnote (1) below
	Ep1 := ESeries(EisWeight(p),ellp,S);

	E4 := ESeries(4,ellp,S);
	E6 := ESeries(6,ellp,S);
	q := Parent(E4).1;
	delta := Delta(q);
	
	deltaj := Parent(delta)!1;
	e := [];	
	for i := 0 to n do
		Wi,deltaj := ComputeWi(k0 + i*EisWeight(p),p,delta,deltaj,E4,E6);
		for bis in Wi do
			eis := p^Floor(i*ap(p)/(p+1))*Ep1^(-i)*bis;
			e := Append(e,eis);
		end for;
	end for;
	
	return e,Ep1;

end function;	

// MAIN FUNCTIONS

// Returns matrix A modulo p^m from Step 6 of Algorithm 1 for each k on kseq.
// kseq has weights congruent modulo (p-1).

// Notational change from paper: In Step 1 following Wan we defined j by
// k = k_0 + j(p-1) with 0 <= k_0 < p-1. Here we replace j by kdiv so that
// we may use j as a column index for matrices.

function Level1UpGj(p,kseq,m) 
 
	// Step 1

	k0 := kseq[1] mod EisWeight(p);
	n := Floor(((p+1)/EisWeight(p)) * (m+1+(ap(p)-1)/(p+1)));
	ell := DimensionMk(k0 + n*EisWeight(p));
	ellp := ell*p;
	mdash := m + Ceiling(n*ap(p)/(p+1));

	// Steps 2 and 3

	e,Ep1 := KatzExpansions(k0,p,ellp,mdash,n);
	
	// Step 4

	q := Parent(Ep1).1; 
	Ep1p := Evaluate(Ep1,q^p);
	Ep1pinv := Ep1p^-1;
	G := Ep1*Ep1pinv;
	Aseq:=[];
	
	for k in kseq do
		kdiv := k div EisWeight(p);
		Gkdiv := G^kdiv;
		u := [];
		for i := 1 to ell do
			ei := e[i];
			ui := Gkdiv*ei;
			u := Append(u,ui);
		end for;
	
		// Step 5 and computation of T in Step 6

		S := Parent(Coefficient(e[1],0));
		T := ZeroMatrix(S,ell,ell);
	
		for i := 1 to ell do
			for j := 1 to ell do
				T[i,j] := Coefficient(u[i],p*(j-1));
			end for;
		end for;
	
		// Step 6: solve T = AE using fact E is upper triangular.
		// Warning: assumes that T = AE (rather than pT = AE) has 
		// a solution over Z/(p^mdash). This has always been the case in 
		// examples computed by the author, see Note 3.1.

		A := ZeroMatrix(S,ell,ell);
	
		for i := 1 to ell do
			Ti := T[i];		
			for j := 1 to ell do
				ej := Parent(Ti)![Coefficient(e[j],l-1): l in [1 .. ell]];
				lj := Integers()!ej[j];
				A[i,j] := S!((Integers()!Ti[j])/lj);
				Ti := Ti - A[i,j]*ej;
			end for;
		end for;			

		// Append A mod p^m to Aseq.
		
		// Append(~Aseq,MatrixRing(pAdicRing(p,m),ell)!A);
		
		Append(~Aseq,MatrixRing(IntegerRing(p^m),ell)!A);

	end for;

	return Aseq,e,mdash;

end function;

//  *** MAIN CODE FOR GENERAL LEVEL ***:

// Returns matrix A modulo p^m from Step 6 of Algorithm 1 (N = 1) or
// Algorithm 2 (N > 1).

function UpGj(p,chi,kseq,m,weightbound)

	N := Modulus(chi);
	if N gt 1 then
		return HigherLevelUpGj(p,chi,kseq,m,weightbound);
	else
		return Level1UpGj(p,kseq,m);
	end if;

end function;

// Returns reverse characteristic series modulo (p^m,q^(d+1)) of the U_p 
// operator on the space of p-adic overconvergent modular forms of level N and 
// weight k, for each k in kseq. Here d is the optional parameter degreebound. 
// If d = Infinity() then output is correct modulo (p^m).
// The optional parameter weightbound is explained in the preamble.

function HeckeSeries(p,chi,kseq,m:weightbound:=WtB(chi))

	t0 := Cputime();
    oneweight:=false;
    // Convert single weight to list
    if Type(kseq) eq RngIntElt then
        kseq:=[kseq];
        oneweight:=true;
    end if;
	
	if p eq 2 and kseq[1] eq 1 then
		print "WARNING: Possible output for p=2, k=1 should not be trusted!";
	end if;
	
	if Type(chi) eq RngIntElt then
		N := chi;
		chi := DirichletGroup(N)!1;
	else
		N := Integers()!Modulus(chi);
	end if;
	
	// The algorithm may terminate with false output without the following check:
	assert IsValidWeightSequence(kseq,p); // Verify klist is a valid input
	//assert ((N mod p) ne 0); // p does not divide N
	
	// For list of odd weights return all 1s sequence
	if (-1)^kseq[1]*chi(-1) eq -1 then
		return [1: i in [1 .. #kseq]];
	end if;
	
	Zpmt<t>:=PolynomialRing(pAdicRing(p,m));
	Aseq:=UpGj(p,chi,kseq,m,weightbound);
	Pseq:=[];
	for A in Aseq do
		P:=ReverseCharacteristicPolynomial(A);
		Append(~Pseq,Zpmt!P);
	end for;
	
	print "Total time was ", Cputime(t0);

    if oneweight eq true then
        return Pseq[1];
    else
        return Pseq;
    end if;
    
end function;



// ----------------------------------------------------------------------
// --|                          Silly Tests                           |--
// ----------------------------------------------------------------------

// Outputs the valuations of the entries of the U_p matrix in level 1.
// This can be used to verify that our bounds analogous to [Wan, Lemma 3.1] are correct and sharp for p = 2,3.

function ValCheck(p,N,k,m);

	chi := DirichletGroup(N).0;
	U := UpGj(p,chi,[k],m,6)[1];
	ell := Dimension(Parent(U[1]));
	V := ZeroMatrix(Integers(),ell,ell);
	for i := 1 to ell do
		for j := 1 to ell do
			if U[i,j] ne 0 then
				V[i,j] := Integers()!Valuation(Integers()!U[i,j],p);
			end if;
		end for;
	end for;
	
	return V;
	
	
end function;

// Outputs the list of slopes of U_p, which can then be compared to the conjectural formulae by Buzzard if p is Gamma0(N)-regular.

function SlopeCheck(p,chi,L,m:weightbound:=WtB(chi));

	// Convert single weight to list
	if Type(L) eq RngIntElt then
		L:=[L];
	end if;

	if Type(chi) eq RngIntElt then
		N := chi;
		chi := DirichletGroup(N)!1;
	else
		N := Integers()!Modulus(chi);
	end if;
	
	if (-1)^L[1]*chi(-1) eq -1 then
		return 0;
	end if;

	// The algorithm may terminate with false output without the following check:
	assert IsValidWeightSequence(L,p); // Verify klist is a valid input

	Aseq:=UpGj(p,chi,L,m,weightbound);
	for i := 1 to #L do
		print "Weight k = ", L[i];
		A := ChangeRing(Aseq[i],Integers());
		f := CharacteristicPolynomial(A);
		N := NewtonPolygon(f,p);
		V := Vertices(N);
		for i := 1 to #V-1 do
			print (V[i+1][2]-V[i][2])/(V[i][1]-V[i+1][1]), V[i+1][1] - V[i][1];
		end for;
		print 0, Degree(f)-V[#V][1];
	end for;
	return "OK";

end function;




function GetSlopePair(p,NN,L,m);

	// Convert single weight to list
	if Type(L) eq RngIntElt then
		kseq:=[L];
	end if;

	AllSl := [];
	for j := 1 to #L do
		Aseq:=UpGj(p,NN,[L[j]],m,6);
		Sl := [];
		A := ChangeRing(Aseq[1],Integers());
		f := CharacteristicPolynomial(A);
		N := NewtonPolygon(f,p);
		V := Vertices(N);
		for i := 1 to #V-1 do
			s := (V[i+1][2]-V[i][2])/(V[i][1]-V[i+1][1]);
			if s le L[j] then
				Append(~Sl,[s, V[i+1][1] - V[i][1]]);
			end if;
		end for;
		Append(~Sl,[0, Degree(f)-V[#V][1]]);
		Append(~AllSl,[Sl]);
	end for;
	return AllSl;

end function;

// Auxiliary function

function projector_2(A,lambda,p,m)

	dim:=Degree(Parent(A));
	R:=IntegerRing(p^m); // m^2 -> m? does it give same output?
	A:=ChangeRing(A,R);
	l:=R!(Integers()!lambda);
	B:=A - l;
	C,Q:=EchelonForm(B);
	assert (Q*B eq C);	
	S:=Q[dim];

	return Vector([Rationals()!(Integers()!Q[dim,i]): i in [1 .. dim]]);

end function;

// This function calculates the q-expansion of the n-th overconvergent eigenform (for all n in L), where the eigenforms are ordered by increasing slopes. It assumes that every slope occurs with multiplicity one, which is false for primes of higher genus. 

function eigenfunctions(p,N,k,m,L)

	time0 := Cputime();
	LBound := Maximum(L);

	k0 := k mod EisWeight(p);
	n := Floor(((p+1)/EisWeight(p)) * (m+1+(ap(p)*p)/(p+1)));
	ell := DimensionMk(k0 + n*EisWeight(p));
	assert (L[#L] le ell); //You need more precision to find that eigenfunction!
	print "Size of matrix for Up is ", ell;	  

	A,Katz := UpGj(p,N,[k],m,6);
	A := ChangeRing(A[1],Integers());
	qExp := ChangeRing(Parent(Katz[2]),pAdicField(p,m));
	f := ReverseCharacteristicPolynomial(A);
	R := Roots(ReciprocalPolynomial(f), pAdicField(p));
	Eigenfns := [];

	print "Size of matrix for Up is ", Degree(f);	  
	print "Precision OK, total matrix computation time: ", Cputime() - time0;

	for n in L do
		print "The slope ", Valuation(R[Degree(f)+1-n][1]), " is rigorously correct: ", (Valuation(R[Degree(f)+1-n][1]) lt (Degree(f) + n*(n+1)/2)*(p-1)/(p+1) -n);
		if n eq 0 then
			Append(~Eigenfns,1);
			else
			lambda := R[Degree(f)+1-n][1];
			// Compute the projector, the bottleneck.
			t0 := Cputime();
			// P := projector(A,lambda);
			P:=projector_2(A,lambda,p,m);
			t1 := Cputime();
			print "Projector time: ", t1 - t0;

			S := qExp!0;
			for i:=1 to Degree(f) do
				S := S + (pAdicField(p,m)!P[i])*(qExp!Katz[i]);
			end for;
			Append(~Eigenfns,S/Coefficient(S,1));
		end if;
	end for;
	time1 := Cputime();
	print "Total time: ", time1-time0;

	return Eigenfns;

end function;



