# Version 1.0
# 2018-03-12
# Author: Catherine Hsu
# Adapted from Bartosz Naskręcki
#
# This code is written for MAGMA.

#*****************************************************************************
#       Copyright (C) 2016 Catherine Hsu <cathyh@uoregon.edu>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#    This code is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#    General Public License for more details.
#
#  The full text of the GPL is available at:
#
#                  http://www.gnu.org/licenses/
#*****************************************************************************


//We are only considering weight 2 below
function GetFactors(N)
	f:=Factorization(N);
	g:=1;
	for factorf in f do
		g:=g*(factorf[1]+1);
	end for;
	m:= EulerPhi(N)/6;
	m:=Numerator(m)*g;
	s:=Factorization(m);
	r:=[ [f[1],f[2]]: f in s];

	return r;
	
end function; 

function SmallestCoprime(B,N)
	P:=PrimesUpTo(B);
	for p in P do
		if GCD(p,N) eq 1 then
			return p;
		end if;
	end for;
end function;

//Level must be a squarefree integer with an odd number of factors

//name - the name of the file to print congruences data
//The function which uses the valuation approach.
function Congruence(N,name)
	factors:=GetFactors(N);
	divisors:=Divisors(N);
	B:=Ceiling(Index(Gamma0(N*N))/6);
	PrimeHeckeOperators:=PrimesUpTo(B);
	//print(#PrimeHeckeOperators);
	Q := Rationals();
	
R<z> := PowerSeriesRing(Q);
	
AssertAttribute(R, "Precision", B); 
	E<q> := -Eisenstein(2, z)/24;
	for inddiv in [2..#divisors] do
		d:= divisors[inddiv];
		//print(E);
		//print(-MoebiusMu(d)*d*Eisenstein(2, z^d)/24);
		E:= E-MoebiusMu(d)*d*Eisenstein(2, z^d)/24;
	end for;
	//print(E);
	allnewforms:=[];
	for level in divisors do
		M:=level;
		Append(~allnewforms,Newforms(CuspForms(Gamma0(M),2)));
	end for;

	Collection:=[];
	Sum:=[];
	Valuations:=[];
	
	for lfactors in factors do
		l:=lfactors[1];
		e:=lfactors[2];
		//exclude the primes 2,3 for now
	if l gt 3 then
		Append(~Valuations, [l,e]);
		PrintFile(name, "\\vspace{.2cm}");
		PrintFile(name, "$N=" cat IntegerToString(N) cat ",$ $p=" cat IntegerToString(l) cat "$, $\\mathrm{val}_p(\\eta(N))=" cat IntegerToString(e) cat "$:");
		PrintFile(name, "\\[\\begin{array}{|c|c|c|c|c|} \\hline\\mathrm{level} &\\mathrm{depth} & \\mathrm{ramindex} & \\mathrm{resfield} & \\mathrm{conjclass} \\\\\\hline");
		for ind0 in [1..#divisors] do
			M := divisors[ind0];
			S:=allnewforms[ind0];
			print("Newforms of level M=" cat IntegerToString(M) cat " computed");
				
			for ind1 in [1..#S] do
				print("Computing congruence with form i=" cat IntegerToString(ind1) cat ", j=" cat IntegerToString(1));
				//Hardcoded 1
				numconj:=#S[ind1];
				f:=S[ind1][1];
				K:=CoefficientField(f);
				print("Working with number field");
				print(K);
				d:=Discriminant(K);
				if (d le 10^14000) then
					if not(K eq Rationals()) then
						O:=EquationOrder(K);
						O:=pMaximalOrder(O,l);
					end if;
				
					if (K eq Rationals()) then
						O:=Integers();
					end if;
				
					print("pMaximal order computed for M=" cat IntegerToString(M) cat " and l=" cat IntegerToString(l)); 
				
					I:=ideal<O|l>;
			
			
					factI:=Decomposition(O,l);
				
					for ind3 in [1..#factI] do
						ptuple:=factI[ind3];
					
						J:=ptuple[1];
						ramindex:=ptuple[2];
						resfield:=#ResidueClassField(J);

						if (K eq Rationals()) then
							p1:=SmallestCoprime(B,N);
							//print(p1);
							cf:=Coefficient(f,p1);
							nf:=Numerator(cf);
							df:=Denominator(cf);
							ce:=Coefficient(E,p1);
							In:=nf-ce*df;
							Id:=df;
							if In ne 0 then
								max_exponent:=Valuation(In,l)-Valuation(Id,l);
							else 
								max_exponent:= Infinity();
							end if;
							congruenceExp:=max_exponent;
							//print(congruenceExp);
							
							counter:=0;
							for i in PrimeHeckeOperators do
								counter:= counter+1;
								if GCD(i,N) eq 1 then
									cf:=Coefficient(f,i);
									nf:=Numerator(cf);
									df:=Denominator(cf);
									ce:=Coefficient(E,i);
									In:=nf-ce*df;
									Id:=df;
									val:=Valuation(In,l)-Valuation(Id,l);
									tmp_congruenceExp:=Minimum(congruenceExp,val);
									congruenceExp:=tmp_congruenceExp;
									//print(congruenceExp);
									//print(counter);
	
								if (congruenceExp eq 0) then

										break i;
									end if;

								end if;
							end for;
							pExponent:=Valuation(l,l);
						else	
							p1:=SmallestCoprime(B,N);
							cf:=Coefficient(f,p1);
							nf:=Numerator(cf);
							df:=Denominator(cf);
							ce:=Coefficient(E,p1);
							In:=ideal<O|nf-ce*df>;
							Id:=ideal<O|df>;
							max_exponent:=Valuation(In,J)-Valuation(Id,J);
							congruenceExp:=max_exponent;
							for i in PrimeHeckeOperators do
								if GCD(i,N) eq 1 then
									cf:=Coefficient(f,i);
									nf:=Numerator(cf);
									df:=Denominator(cf);
									ce:=Coefficient(E,i);
									In:=ideal<O|nf-ce*df>;
									Id:=ideal<O|df>;
									val:=Valuation(In,J)-Valuation(Id,J);
									tmp_congruenceExp:=Minimum(congruenceExp,val);
									congruenceExp:=tmp_congruenceExp;
									
									if (congruenceExp eq 0) then
	
										break i;
									end if;
								end if;
							end for;
				
							pExponent:=Valuation(I,J);
				
						end if;
						if (congruenceExp gt 0) and IsFinite(congruenceExp) then				
							Append(~Collection,[M,l,2,ind1,1,ind3]);
							Append(~Sum,[M,l,congruenceExp,ramindex,resfield]);
							print("CONGRUENCE FOUND! with r=" cat IntegerToString(congruenceExp));
							list:=[M,l,e,ind1,1,ind3];
							exp:=[ x[2] : x in factI];
							expprint:=&cat[IntegerToString(x) cat " " : x in exp];
						
							if (pExponent gt 1) then
								ram:="t";
							else
								ram:="n";
							end if;
							PrintFile(name, IntegerToString(M) cat "&" cat IntegerToString(congruenceExp) cat "&" cat IntegerToString(ramindex) cat "&" cat IntegerToString(resfield) cat "&" cat IntegerToString(ind1)cat "\\\\");
						end if;
					end for;
				end if; 
			end for;
		
		end for;
		PrintFile(name,"\\hline\\end{array}\\]");
		PrintFile(name,"");
	end if;
	end for;

	return Collection, Sum, Valuations;

end function;