/*********************************************************/
/*                       MyMF.gp                         */
/*							 */
/*                 von R.Steffens			 */
/*             geschrieben am: 20.10.05			 */
/*             letzte Änderung: 23.11.05		 */
/*********************************************************/


MyMF(N,n=20,r=1,s=1)=
{
	v=vector(n,k,0);
	M=makeSL2Z(r,s);
	print(M);
	if(kronecker(-4*N,M[1,1])*chi(M[1,1])!=0,
	for(D=0,n-1,
		v[D+1]=Koeffizient(N,D/N,M);
	);
	);
	return(v);
}


Koeffizient(N,D,M=[1,0;0,1])=
{
	KF=0;
	Sb=floor(sqrt(D)*abs(M[1,1]*M[2,2]+M[1,2]*M[2,1])+2*D*N*(abs(M[1,2]*M[2,2])+abs(M[1,1]*M[2,1])))+20;

	if((type(D)=="t_INT")&&(issquare(D)==1),flg=1,flg=0);
	forstep(b=-Sb,Sb,1,
		if(flg==0,
			div=divisors(N*b^2-N*D);
			for(j=1,2*length(div),
				if(j%2==0,a=div[j/2],a=-div[(j+1)/2]);
				c=(b^2-D)/(4*a);
				KF=KF+ComputeKF(N,a,b,c,M);
			);
		);
		if(flg==1,
			if(b^2-D != 0,
				div=divisors(N*b^2-N*D);
				for(i=1,2*length(div),
					if(i%2==0,a=div[i/2],a=-div[(i+1)/2]);
					c=(b^2-D)/(4*a);
					KF=KF+ComputeKF(N,a,b,c,M);	
				);
			,
				/*a=0;*/
				Sc=floor(abs(M[1,1]*M[1,2])*sqrt(D)+M[1,1]^2)+10;
				forstep(c=-Sc,Sc,1/(4*N),
					KF=KF+ComputeKF(N,0,b,c,M);
				);
			);
		);
	);
	return(KF);
}


ComputeKF(N,a,b,c,M)=
{
	
	aM=M[1,1]^2*a+M[2,1]*M[1,1]*b+M[2,1]^2*c;
	cM=M[1,2]^2*a+M[2,2]*M[1,2]*b+M[2,2]^2*c;

	if(aM*cM<0,SM=sign(aM),
	if(aM*cM>0,SM=0,
	if((cM==0)&&(0<aM)&&(aM<1),SM=1/2-aM,
	if((aM==0)&&(0<cM)&&(cM<1),SM=-1/2+cM,SM=0))));
		
	return(SM*kronecker(-4*N,4*N*c)*chi(4*N*c));
}

makeSL2Z(a,c)=
{
	if(gcd(a,c)==1,A=[a,-bezout(a,c)[2];c,bezout(a,c)[1]],return("ERROR")); 
	return(A);
}


koblitz(b,a)=
{
	if(gcd(b,a)!=1,return(0));
	if(issquarefree(b)==0,return(kronecker(b,a)),
	m=b%4;
	if(m==1,
		return(kronecker(a,abs(b)));
	,	
		if(m==3,return(kronecker(-1,a)*kronecker(a,abs(b))),
			if((b/2)%4==1,return(kronecker(2,a)*kronecker(a,b/2)),
				return(kronecker(2,a)*kronecker(-1,a)*kronecker(a,abs(b)/2))
			);
		);
	);
	);
}

shimura(a,b)=
{
	if(gcd(a,b)!=1, return(0));
	if((b>0)&&(isprime(b)==1)&&(b!=2),return(kronecker(a,b)));
	if((a==0)&&((b==1)||(b==-1)),return(1));
	if((b=-1),if(a>=0,return(1),return(-1)));
	div=divisors(b);
	ret=1;
	for(i=1,length(div),
		ret=ret*kronecker(a,div[i]);
	);
	if(b>0,return(ret),return(sign(a)*ret));	
}


Hecke(N,n=20,r=1,s=1,p=3)=
{
	A=makeSL2Z(r,s);
	v=vector(n,k,0);
	print(A);
        
	if(isprime(p)==0,print(p" ist nicht prim");return(0));
	if(N%p==0, print(p" teilt "N);return(0));
	
	
	for(D=0,n-1,
		/*Hilfsvariablen zurücksetzen        */
		HF1=0;
		HF2=0;
		HF3=0;
			/*a_{n/p^2}=0 falls p^2 n nicht teilt      */
		if((N*D)%(p^2)==0,
			HF1=Koeffizient(N,D/(N*p^2),A);
			s=p*HF1,
			s=0;
		);
			
		/*Berechnung von a_{n*p^2}                 */
		HF2=Koeffizient(N,D/N*p^2,A);
		
		/*Berechnung von a_{n}                     */
		HF3=Koeffizient(N,D/N,A); 
		
		/*Zusammensetzen der berechneten Werte zum Koeffizienten */
		v[D+1]=HF2+chi(p)*kronecker(-D,p)*HF3+s*chi(p^2);
	);
	return(v);
}


chi(n)=
{
return(kronecker(1,n));
}