/* The following routines have been written by Alvaro Lozano-Robledo */
/* Please notify me (by email) about questions, comments and corrections */
/* Feel free to use and modify them as you please! */

/* Recall that bernfrac(k) produces the kth (classical) Bernoulli number */
/* The following routines calculate: */
/*	a) Generalized Bernoulli numbers for quadratic Dirichlet characters modulo N */
/*	b) Values of Dirichlet L-functions for quadratic Dirichlet characters modulo N */
/*	c) Values of Dedekind Zeta Functions for real quadratic number fields K=Q(\sqrt{N}) */

/* genbern(k,N) is the kth generalized Bernoulli number for the quadratic Dirichlet character modulo N */

genbern(k,N)=
{local(SOL,expan,chi);
default(seriesprecision,k+5);
expan=t*exp(t);
for(i=2,N-1,
	chi=kronecker(N,i);
	expan=expan+chi*t*exp(i*t);
);
expan=expan/(exp(N*t)-1);
SOL=polcoeff(expan,k,t)*k!;
SOL
}

/* genbernn(k,N) finds the first k generalized Bernoulli number */

genbernn(k,N)=
{local(SOL,expan,chi);
SOL=[];
default(seriesprecision,k+5);
expan=t*exp(t);
for(i=2,N-1,
	chi=kronecker(N,i);
	expan=expan+chi*t*exp(i*t);
);
expan=expan/(exp(N*t)-1);
for(j=1,k,
	SOL=concat(SOL,[polcoeff(expan,j,t)*j!]);
);
SOL
}


/* lfunct(s,N) calculates L(s,chi) where s is a negative integer and chi is a quadratic character modulo N */

lfunct(s,N)=
{local(SOL);
SOL=-genbern(1-s,N)/(1-s);
SOL
}

/* zdedek(k,N) calculates Z(1-k,K) where k is a positive integer and K is a real quadratic number field of */
/*  discriminant N */

zdedek(k,N)=
{local(SOL);
SOL=genbern(k,N)*bernfrac(k)/(k)^2;
SOL
}

/* zdedekn(m,N), m positive, calculates all values Z(1-j,K) where j, 2 <= j <= m*/

zdedekn(m,N)=
{local(SOL,genb);
genb=genbernn(m,N);
SOL=[0];
for(j=2,m,
	SOL=concat(SOL,[genb[j]*bernfrac(j)/(j)^2]);
);
SOL
}

/* writezdedekn(m,N,file), m positive, calculates all values Z(1-j,K) where j, 2 <= j <= m*/  
/*    and prints to file */ 

writezdedekn(m,N,file)=
{local(Z,genb);
genb=genbernn(m,N);
Z=[0];
for(j=2,m,
	Z=concat(Z,[genb[j]*bernfrac(j)/(j)^2]);
);
write(file,"print(\" \");");
write(file,"print(\" Vector Z uploaded to memory, length(Z)=\",", length(Z),");"); 
write(file,"print(\" Z[j] returns the value Z(1-j,K) where K is the real\" );");
write(file,"print(\" quadratic number field of discriminant \",", N,");");
write(file,"print(\" and Z(s,K) is the Dedekind zeta function of K \" );");
write(file,"print(\" \");");
write(file,"Z=",Z,";");
Z  
}

/* writezdedeknn(m,N,file), m positive, calculates all values Z(1-j,K) where j, 2 <= j <= m, */  
/*  and here N is a vector of discriminants, and prints to file */ 

writezdedeknn(m,N,file)=
{local(Z,genb,Ni);
write(file,"N=",N,";");
write(file,"print(\" Vectors ZN uploaded to memory, for discriminants N in the list \",N );");
write(file,"print(\" \");");
write(file,"print(\" ZN[j], j<=",m,", returns the value Z(1-j,K) where K is the real\" );");
write(file,"print(\" quadratic number field of discriminant N \");");
write(file,"print(\" and Z(s,K) is the Dedekind zeta function of K \" );");
write(file,"print(\" \");");
for(i=1,length(N),
	Ni=N[i];
	genb=genbernn(m,Ni);
	Z=[0];
	for(j=2,m,
		Z=concat(Z,[genb[j]*bernfrac(j)/(j)^2]);
	);
	write(file,"Z",Ni,"=",Z,";");
	write(file," ");
);
Z="Done.";
Z  
}