Z:=Integers();
R<q>:=PowerSeriesRing(Z);
M:=Matrix([[1,0],[0,2]]);
L:=LatticeWithGram(M);
R<q>:=PowerSeriesRing(Z);
T:=R!ThetaSeries(L,100);

function U(f,d)
	qq:=Parent(f).1;
	start:=Valuation(f);
	en:=AbsolutePrecision(f)-1;
	ret:=0;
	for n in [start..en] do
		if (n mod d eq 0) then
			ret := ret + Coefficient(f,n)*qq^(Floor(n/d));
		end if;
	end for;
	ret := ret + BigO(qq^(Ceiling(AbsolutePrecision(f)/d)));
	return ret;
end function;

forward represents;

minrept:=function(mlist,min,max)
	tlist:=[];
	for a in [1..#mlist] do
		L:=LatticeWithGram(2*mlist[a]);
		T:=R!ThetaSeries(L,max*2);
		T:=U(T,2);
		Append(~tlist,T);
	end for;
	y:=1;
	oldmin:=min;
	while(y le #mlist) do
		while(Coefficient(tlist[y],min) eq 0) do
			min:=min+1;
			c:=true;
			while c do
				if IsPrime(min) then
					min:=min+1;
				else if IsSquarefree(min) then
					c:=false;
				else
					min:=min+1;
				end if;
				end if;
			end while;
			if(min ge max) then
				return -1;
			end if;
		end while;
		if(oldmin ne min) then
			oldmin:=min;
			y:=1; 
		else
			y:=y+1;
		end if;
	end while;
	return min;
end function;

minrep:=function(mlist,min,max)
	y:=2;
	while (y le #mlist) do
		k:=represents(mlist[1],mlist[y],min,max);
		if(k eq -1) then
			return -1;
		else if(k eq min) then
			y:=y+1;
		else 
			min;
			min:=k;
			y:=2;
		end if;
		end if;
	end while;
	return k;
end function;	

represents:=function(S,D,min,max)
	L1:=LatticeWithGram(2*S);
	T1:=R!ThetaSeries(L1,max*2);
	T1:=U(T1,2);
	L:=LatticeWithGram(2*D);
	T:=R!ThetaSeries(L,max*2);
	T:=U(T,2);
	for n in [min..max] do
		if Coefficient(T,n) ne 0 then
			if Coefficient(T1,n) ne 0 then
				if IsSquarefree(n) then
					if IsPrime(n) then
					else
						return n;
					end if;
				end if;
			end if;
		end if;
	end for;
	return -1;
end function;	

truant:=function(M,max)
	L:=LatticeWithGram(2*M);
	T:=R!ThetaSeries(L,max*2);
	T:=U(T,2);
	min:=-7;
	for n in [1..max] do
		if Coefficient(T,n) eq 0 then
			if Gcd(3,n) ne 3 then
				min:=n;
				break;
			end if;
		end if;
	end for;
	return min;
end function;

function cleanrepeats(mlist)
	rep:=false;
	for i in [1..#mlist-1] do
		if IsIsometric(LatticeWithGram(2*mlist[i]),LatticeWithGram(2*mlist[#mlist])) then 
			rep := true;
		 	break;
		end if;
	end for;
	return rep;
end function;;

//escalate every form in mlist, useTruant defaults to first truant, but can be changed to second
//to accomodate the 10-14 switch
function escalate(mlist : useTruant:=1, donelist:=[[]:b in [1..99999]])
	whentoprint:=0;
	numprint:=0;
	possible:={i:i in [1..453]};
	donelistRE:=[[]:b in [1..99999]];
	notruant:=[];
	esclist:=[];
	rankequal:=[];
	j:=1;
	while (j le #mlist) do
		M:=mlist[j];
		width:=NumberOfColumns(M);
		if (width ge 2) then
			missing:=truant(M,1000);
		else
			missing:=2;
		end if;
		while(missing lt 0) do
			Append(~notruant,M);
			j:=j+1;
			if j ge #mlist then
					return esclist,rankequal,notruant;
			end if;
			M:=mlist[j];
			missing:=truant(M,1000);
		end while;
		if (useTruant ge 2) then
			missing:=useTruant;
		end if;
		M:=ChangeRing(M,Rationals());
		N:=DiagonalJoin(M,Matrix(Rationals(),[[missing]]));
		N:=ChangeRing(N,Rationals());
		width:=NumberOfColumns(N);
		if width eq 4 then
			limits1:=Floor(2*Sqrt(N[1][1]*missing));
			limits2:=Floor(2*Sqrt(N[2][2]*missing));
			limits3:=Floor(2*Sqrt(N[3][3]*missing));
			for a in [-limits1..limits1] do
				N[1][4]:=a/2;
				N[4][1]:=a/2;
				for b in [-limits2..limits2] do
					N[2][4]:=b/2;
					N[4][2]:=b/2;
					for c in [0..limits3] do
						N[3][4]:=c/2;
						N[4][3]:=c/2;
						if(IsPositiveSemiDefinite(N)) then
							#esclist;
							P,temp,rank:=LLLGram(N);
							if(rank eq 3) then
								N2:=Submatrix(P,1,1,rank,rank);
								P2:=LLLGram(N2);	
								d:=Z!Determinant(2*P2);
								Append(~donelistRE[d],P2);
								if(cleanrepeats(donelistRE[d])) then
									Remove(~donelistRE[d],#donelistRE[d]);
								else
									Append(~rankequal,P2);
								end if;
							else
								d:=Z!Determinant(2*P);
								Append(~donelist[d],P);
								if(cleanrepeats(donelist[d])) then
									Remove(~donelist[d],#donelist[d]);
								else
									Append(~esclist,P);
								end if;
							end if;
						end if;
					end for;
				end for;
			end for;
		end if;
		if width eq 2 then
			limits1:=Floor(2*Sqrt(N[1][1]*missing));
			for a in [0..limits1] do
				N[1][2]:=a/2;
				N[2][1]:=a/2;
				if(IsPositiveSemiDefinite(N)) then
					P,temp,rank:=LLLGram(N);
					if(rank eq 1) then
						N2:=Submatrix(P,1,1,rank,rank);
						P2:=LLLGram(N2);	
						d:=Z!Determinant(2*P2);
						Append(~donelist[d],P2);
						if(cleanrepeats(donelist[d])) then
							Remove(~donelist[d],#donelist[d]);
						else
							Append(~rankequal,P2);
						end if;
					else
						d:=Z!Determinant(2*P);
						Append(~donelist[d],P);
						if(cleanrepeats(donelist[d])) then
							Remove(~donelist[d],#donelist[d]);
						else
							Append(~esclist,P);
						end if;
					end if;
				end if;
			end for;
		end if;
		if width eq 3 then
			limits1:=Floor(2*Sqrt(N[1][1]*missing));
			limits2:=Floor(2*Sqrt(N[2][2]*missing));
			for a in [-limits1..limits1] do
				N[1][3]:=a/2;
				N[3][1]:=a/2;
				for b in [0..limits2] do
					N[2][3]:=b/2;
					N[3][2]:=b/2;
					if(IsPositiveSemiDefinite(N)) then
						P,temp,rank:=LLLGram(N);
						if(rank eq 2) then
							N2:=Submatrix(P,1,1,rank,rank);
							P2:=LLLGram(N2);	
							d:=Z!Determinant(2*P2);
							Append(~donelist[d],P2);
							if(cleanrepeats(donelist[d])) then
								Remove(~donelist[d],#donelist[d]);
							else
								Append(~rankequal,P2);
							end if;
						else
							d:=Z!Determinant(2*P);
							Append(~donelist[d],P);
							if(cleanrepeats(donelist[d])) then
								Remove(~donelist[d],#donelist[d]);
							else
								Append(~esclist,P);
							end if;
						end if;
					end if;
				end for;
			end for;
		end if;
		j:=j+1;
	end while;
	if useTruant ge 2 then
		return esclist,rankequal,notruant,donelist;
	else
		return esclist,rankequal,notruant;
	end if;
end function;	

function notrepresented(mlist,max)
	notreplist:=[[]:b in [1..#mlist]];
	y:=1;
	while(y le #mlist) do
		i:=1;
		L:=LatticeWithGram(2*mlist[y]);
		T:=R!ThetaSeries(L,max*2);
		T:=U(T,2);
		while(i le max) do
			if Coefficient(T,i) eq 0 then
				Append(~notreplist[y],i);
			end if;
			i:=i+1;
			if Gcd(i,3) eq 3 then
				i:=i+1;
			end if;
		end while;
		y:=y+1;
	end while;
	return notreplist;
end function;









// Identify which forms have an undesirably long list of unrepresented numbers
formlist:=[
Matrix([[1,0,0,-1/2],[0,1,0,-1/2],[0,0,2,1],[-1/2,-1/2,1,5]]),
Matrix([[1,0,-1/2,0],[0,2,1,0],[-1/2,1,5,1],[0,0,1,4]]),
Matrix([[1,0,-1/2,0],[0,2,-1,0],[-1/2,-1,3,1],[0,0,1,4]]),
Matrix([[1,0,-1/2,0],[0,2,1,0],[-1/2,1,5,3/2],[0,0,3/2,9]]),
Matrix([[1,0,-1/2,-1/2],[0,2,1,-1],[-1/2,1,5,-2],[-1/2,-1,-2,9]]),
Matrix([[1,0,-1/2,0],[0,1,1/2,0],[-1/2,1/2,7,-1],[0,0,-1,10]]),
Matrix([[1,0,-1/2,0],[0,1,1/2,0],[-1/2,1/2,7,3],[0,0,3,10]]),
Matrix([[1,0,-1/2,-1/2],[0,1,1/2,-1/2],[-1/2,1/2,7,-3/2],[-1/2,-1/2,-3/2,7]]),
Matrix([[1,0,-1/2,-1/2],[0,1,-1/2,1/2],[-1/2,-1/2,3,1/2],[-1/2,1/2,1/2,7]]),
Matrix([[1,0,0,-1/2],[0,1,0,-1/2],[0,0,7,7/2],[-1/2,-1/2,7/2,12]]),
Matrix([[1,0,0,0],[0,1,0,0],[0,0,7,7/2],[0,0,7/2,14]])
];

function checkBadEscalation(form)
	for k in [1..#formlist] do
		if(IsIsometric(LatticeWithGram(2*form),LatticeWithGram(2*formlist[k]))) then
			return k;
		end if;
	end for;
	return -1;
end function;


function badEscalation(M)
	badescs:=[];
	width:=NumberOfColumns(M);
	missing:=truant(M,71);
	if(missing lt 0) then
		return badescs;
	end if;
	M:=ChangeRing(M,Rationals());
	N:=DiagonalJoin(M,Matrix(Rationals(),[[missing]]));
	N:=ChangeRing(N,Rationals());
		
	limits1:=Floor(2*Sqrt(N[1][1]*missing));
	limits2:=Floor(2*Sqrt(N[2][2]*missing));
	limits3:=Floor(2*Sqrt(N[3][3]*missing));
	for a in [-limits1..limits1] do
		N[1][4]:=a/2;
		N[4][1]:=a/2;
		for b in [-limits2..limits2] do
			N[2][4]:=b/2;
			N[4][2]:=b/2;
			for c in [0..limits3] do
				N[3][4]:=c/2;
				N[4][3]:=c/2;
				if(IsPositiveSemiDefinite(N)) then
"checking form: ";
N;
					P,temp,rank:=LLLGram(N);
					if rank ne 3 then
						d:=Z!Determinant(2*P);
						bad:=checkBadEscalation(P);
						if bad ne -1 then
							Append(~badescs,checkBadEscalation(P));
						end if;
					end if;
				end if;
			end for;
		end for;
	end for;
	return badescs;
end function;

function IDbadEscalations(tern_list)
	j:=1;
	badterns:=[];
	badescs:=[];
	while (j le #tern_list) do
		badEscList := badEscalation(tern_list[j]);
		if (#badEscList ne 0) then
			Append(~badterns,tern_list[j]);
			Append(~badescs,badEscList);
		end if;
		j:=j+1;
	end while;	
	return badterns, badescs;
end function;

// Writes list of matrices passed in M to a file specified by name.
// This file is then loadable by Magma.
function writeMtoFile(M,name)
	fprintf name,"mlist:=[";
	for i in [1..#M] do
		if i le 1 then
			fprintf name,"Matrix([[";
		else
			fprintf name,",Matrix([[";
		end if;
		for j in [1..NumberOfRows(M[i])] do
			for k in [1..NumberOfColumns(M[i])] do
				fprintf name,"%o",M[i][j][k];
				if k eq 4 then
					fprintf name,"]";
				else
					fprintf name,",";
				end if;
			end for;
			if j eq 4 then
				fprintf name,"])\n";
			else
				fprintf name,",[";
			end if;
		end for;
	end for;
	fprintf name,"];";
	return 1;
end function;


// Checks two lists of matrices for isometric matrices
function compareEsclists(list1, list2)
	Same:=[];
	Diff:=[];
	found:=0;
	for k in [1..#list1] do
		j:=1;
		while (j le #list2) do
			if(IsIsometric(LatticeWithGram(2*list2[j]),LatticeWithGram(2*list1[k]))) then
				Append(~Same,list1[k]);
				found:=1;
				j:=#list2+1;
			end if;
			j:=j+1;
		end while;
		if found le 0 then
			Append(~Diff,list1[k]);
		end if;
		found:=0;
	end for;
	return Same,Diff;
end function;