load "localsplit.txt";


//given a number, n, and a prime, p, 
//returns the result of dividing all factors of p from n 
//followed by the highest power of p diving n
function extract_p(n,p)
	n := Integers()!n;
	p := Integers()!p;
	return (n div (p^Valuation(n,p))),Valuation(n,p);
end function;

forward alpha;

//Matrix M, prime p, number t
//diagonal is factored
function local_density(M,p,t)
	e := [];
	l := [];
	k := GetLocal(M,p);
	m := NumberOfRows(k);
	for a in [1..m] do
		e[a],l[a]:=extract_p(k[a][a],p);
	end for;
	answer := alpha(t,p,e,m,l,0);
	return answer;
end function;

function L(l,k,m)
	n:=0;
	Lset := [];
	for a in [1..m] do
		if (l[a]-k lt 0) then
			if ((l[a]-k) mod 2 eq 1) then
				n:=n+1;
				Lset[n] := a;
			end if;
		end if;
	end for;
	return n,Lset;
end function;

function d(k,l,m)
	n:=0;
	for a in [1..m] do
		if (l[a] lt k) then
			n := n + l[a]-k;
		end if;
	end for;
	return (Floor((k + 1/2 * n)));
	
end function;

function v(k,e,p,m,l)
	j,Lset := L(l,k,m);
	r := 1;
	for a in [1..j] do
		r := r * KroneckerSymbol(e[Lset[a]],p);
	end for;
	return KroneckerSymbol(-1,p)^(Floor(j/2))*r;
end function;


//NOTE: this is not actually used in calculations, just for the proof
//function delta(p)
//	p := p mod 4;
//	if (p eq 1) then
//		return 1;
//	end if;
//	if (p eq 3) then
//		return Sqrt(-1);
//	end if;
//end function;


//NOTE: this is not actually used in calculations, just for the proof
function gamma(t,p)
	n,a := extract_p(t,p);
	a := a mod 2;
	if (a eq 0) then
		return 1;
	end if;
	if (a eq 1) then
		return delta(p)*KroneckerSymbol(-n, p);
	end if;
end function;

function f(t,l,m,p)
	n,a := extract_p(t,p);
	if (IsEven(L(l,a+1,m))) then
		return (-1/p);
	end if;
	return KroneckerSymbol(n,p);
end function;

function R(t,p,e,m,l)
	n,a := extract_p(t,p);
	first_term := 0;
	for temp in [1..a] do
		if (IsEven(L(l,temp,m))) then
//		if (IsEven(L(l,a+1,m))) then
			first_term := first_term + v(temp,e,p,m,l) * p^(Integers()!d(temp,l,m));
		end if;
	end for;
	first_term := first_term * (1-1/p);
	second_term := v(a+1,e,p,m,l) * p^(d((a+1),l,m)) * f(t,l,m,p);
	return first_term + second_term;
end function;

function alpha(t,p,e,m,l,l2)
	return 1 + p^(l2) * R(t,p,e,m,l) + (1-1/p) * (l2) * p^(l2);
end function;