load "p2 case2.txt";
load "localdensity.txt";
load "readable_quaternaries.txt";


mlist:=quatver;

R<q> := PowerSeriesRing(Rationals());
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;


testM:=Matrix([[1,1/2,0,0],[1/2,1,0,0],[0,0,7,0],[0,0,0,7]]);
L := LatticeWithGram(2*testM);
T := R!ThetaSeries(L,100);
T := U(T,2);

D := Determinant(2*testM);
chi := KroneckerCharacter(Integers()!D);
chi := DirichletGroup(84)!chi;
MM := ModularForms(chi);

y := EisensteinProjection(MM!T);
z := CuspidalProjection(MM!T);


R:=RealField(30);
function InfinitePlace(M,m)
  n := NumberOfRows(M);
  D := Determinant(M)/(2^n);
  if n eq 3 then
    omegan := (4/3)*Pi(R);
  end if;
  if n eq 4 then
    omegan := (1/2)*Pi(R)^2;
  end if;
  ret := ((n*omegan)/(2*D^(1/2)))*m^((n-2)/2);
  return ret;
end function;


function get_local_density(M,t,p)
	density := local_density_2(M,t)*InfinitePlace(M,t);
	i := 3;
	while(i le p) do
		density:=density*local_density(M,i,t);
		i:=NextPrime(i);
	end while;
	return density;
end function;

function compare_terms(constant1)
	for i in [1..100] do
		print i, Round(get_local_density(testM,i,2000)/4*(constant1)), Coefficient(y,i)*constant1;
	end for;
return 1;
end function;


function check_list()
for h in [1..#mlist] do
	print "Testing matrix:", h;
	testM:=mlist[h];
	L := LatticeWithGram(2*testM);
	R:=Representatives(Genus(L));
	mass:=&+[1/#AutomorphismGroup(R[i]): i in [1..#R]];
	Eis:=&+[ThetaSeries(R[i],1000)/#AutomorphismGroup(R[i]): i in [1..#R]];
	Eis:=Eis/mass;
	y:=Eis;

	for i in [1..100] do
		if(Abs(get_local_density(testM,i,2000)/4 - Coefficient(y,2*i)) ge 1) then
			print testM, h, i, get_local_density(testM,i,2000)/4, Coefficient(y,2*i);
		end if;
	end for;
end for;
end function;