/**************************************************************************

	PARI/GP Module

	PARI/GP Script for caculating the trace of a Hecke-Operator

	(C) Nils-Peter Skoruppa and Reinhard Steffens, 2006

	$Id: TraceHecke.gp,v 1.2 2006/08/06 $

**************************************************************************/
TH_GP = 1;

if( MI_GP != 1, read( "Miscellaneous.gp"));

/************************************************************************
We use Theorem 3 of Oesterle's thesis to calculate the trace T(n) on 
S_k(\Gamma_0(N),\chi) for integers k > 2.
The calculation is done in 4 parts.
A1 is the contribution of scalar matricies
A2 ----------"----------- elliptic matrices
A3 ----------"----------- hyperbolic matrices
A4 ----------"----------- paraolic matrices
*************************************************************************/

TH_TraceHecke(n,k,N,dc)=
/*@@
	Input
		n: positive integer
		k: integer larger than 2 -> weight
		N: positive integer -> level
		dc: Dirichlet Charakter modulo N with dc(-1)=(-1)^k,
		    conductor condDC and order orderDC
		
	
	Output
		the trace of the Hecke-Operator T(n) on S_k(\Gamma_0(N),dc)

*/
{
	local(  A1,A2,A3,A4,
		l1,l2,t_max,summ1,summ2,
		Div,Div2,F,f,discr,b,
		condDC,r,hilf,orderDC  );
	
	orderDC=DC_order(dc);
	r=Mod(x,polcyclo(orderDC)); 
	condDC=DC_conductor(dc);

	
	/*
		Calculation of A1
	*/
	A1=0;
	if(issquare(n)==1,
		A1=n^(k/2-1)*lift(DC_evaluate(dc,r,round(sqrt(n))))*(k-1)/12*TH_phi1(N);
	);
	


	/*
		Calculation of A2
	*/	
	summ1=0;
	t_max=0;
	while(t_max^2-4*n<0,
		t_max=t_max+1;
	);
	t_max=t_max-1;
	for(t=0,t_max,
		F=TH_get_F(t,n);
		Div=divisors(F);
		summ2=0;
		for(i=1,length(Div),
			f=Div[i];
			discr=(t^2-4*n)*f^(-2);
			summ2=summ2 + qfbhclassno(discr) / TH_w(discr)
				*TH_mu(t,f,n,N,dc,r);
		);
		summ1=summ1 - TH_p(k,t,n) * summ2;
	);
	
	A2=summ1;


	/*
		Calculation of A3
	*/
	summ1=0;

	Div=divisors(n);
	for(i=1,length(Div),
		l1=Div[i];
		l2=n/l1;
		if(l1>l2,
			summ2=0;
			Div2=divisors(l1-l2);
			for(j=1,length(Div2),
				b=Div2[j];
				summ2=summ2+eulerphi((l1-l2)/b)
					*TH_mu(l1+l2,b,n,N,dc,r);
			);
		summ1=summ1 - l2^(k-1)/(l1-l2)*summ2;
		);
	);
	A3=summ1;

	
	/*
		Calculation of A4
	*/
	summ1=0;
	
	A4=0;
	if(issquare(n)==1,
		Div=divisors(N);
		for(i=1,length(Div),
			c=Div[i];
			if((N/condDC)%(gcd(c,N/c))==0,
				summ1=summ1+eulerphi(gcd(c,N/c));
			);
		);
		A4=-1/2*n^((k-1)/2)*lift(DC_evaluate(dc,r,round(sqrt(n))))*summ1;	
	);

	return(A1+A2+A3+A4);
}


TH_get_F(t,n)=
/*@@
	Input
		t: integer such that t^2-4n<0
		n: positive integer
		
	Output
		the integer F sucht that (t^2-4n)*F^(-2) is the discriminant
		of an imaginary quadratic field K

*/
{
	local(discr,F,Div);
	
	discr=t^2-4*n;
	F=1;
	Div=factor(discr);
	for(i=1,length(Div[,1]),
		if(Div[i,2]>=2,
			if(Div[i,2]%2==0,
				F=F*Div[i,1]^(Div[i,2]/2);
			,	F=F*Div[i,1]^((Div[i,2]-1)/2);
			);	
		);	
	);
	discr=discr/(F^2);
	if(discr%4==1,return(F),return(F/2));
}


TH_p(k,t,n)=
/*@@
	Input
		k: integer > 2 ->weight
		t: integer
		n: positive integer

	Output
		the number (r-conj(r))^(-1)*(r^(k-1)-conj(r)^(k-1)) where r and
		conj(r) are the roots of x^2-tx+n
		this equals the coefficient of x^(k-2) in the Taylor-Expansion
		of 1/(1-t*x+n*x^2) around 0 

*/
{
	if( ("t_INT" != type(k)) || (k < 2), error( k" must be a nonnegative
	integer larger then 2"));
	if( k == 2, return( 1));
	if( k == 3, return( t));
	
	return( t*TH_p(k-1, t, n) - n*TH_p( k-2, t, n));
}


TH_phi1(N)=
/*@@
	Input
		N: positive integer

	Output
		the number N* \prod_{p|N}(1+1/p)
		
*/
{
	local(div,product);
	
	div=factorint(N);
	product=1;
	for(i=1,length(div[,1]),
		product=product*(1+1/div[i,1]);
	);
	return(N*product);
}

TH_init_E(f,N,t,n)=
/*@@
	Input
		f: integer, divisor of F
		N: positive integer -> level
		t: integer st t^2-4*n<0
		n: integer

	Output
		the vector of equivalence classes x mod N such that
		gcd(fN,N^2) devides x^2-tx+n
		The equivalence classes are represented by Numbers between 0
		and N-1

*/
{
	local(counter,E_vec_pre,E_vec,Div);
	
	counter=0;
	E_vec_pre=vector(N);
	Div=N*gcd(f,N);
	
	for(i=0,N-1,
		if((i^2-t*i+n)%Div==0,
			E_vec_pre[counter+1]=i;
			counter+=1;
		);
	);

	E_vec=vector(counter);

	for(j=1,counter,
		E_vec[j]=E_vec_pre[j];
	);
	return(E_vec);
}


TH_mu(t,f,n,N,dc,r)=
/*@@
	Input
		t: integer
		f: integer, divisor of F
		n: integer
		N: positive integer -> level
		dc: Dirichlet Character modulo N
		r:  an element of an abelian group whose order divides the order of dc, typically a POLYMOD defining an e-th root of unity (where e is the order of dc)	

	Output
		the number TH_phi1(N) * TH_phi1(N/gcd(N,f))^(-1) * \sum_{x \in
		E} dc(x)

*/
{
	local(summ,wert,E);

	if(((t^2-4*n)/(f^2))%4==2||(((t^2-4*n)/(f^2))%4==3),error("Error in function mu"));
	
	E=TH_init_E(f,N,t,n);
	if(length(E)==0,return(0));

	summ=0;	
	
	for(i=1,length(E),
		summ=summ+lift(DC_evaluate(dc,r,E[i]));
	);
	   
	wert=TH_phi1(N)*TH_phi1(N/gcd(N,f))^(-1)*summ;
	return(wert);
}


TH_w(k)=
/*@@
	Input
		k: negative integer congruent to 0 or 1 modulo 4

	Output
		number of units in the imaginary quadratic field K with
		discriminant k

*/
{
	/* The formula for w wich is used here is from Cohen's book 
	   "Computational Algebraic Number Theory" */

	local(number);

	if(k>=0,error("TH_w(k): k has to be negative"));
	if(k%4==2 || k%4==3,error("TH_w(k): k is not congruent to 0 or 1 mod 4"));

	if(k<-4, number=2,
		 if(k==-4, number=4,
			   if(k==-3, number=6;)
		 );
	);

	return(number);

}