/* INSTRUCTIONS FOR APPROXIMATING THE RESIDUE OF THE DEDEKIND ZETA FUNCTION
   ========================================================================

Use the following steps to compute the Euler product of a number field of
degree 4 (modify it yourself for other degrees, and email it to the
assistants of the course if you do):

1. Open the online Magma calculator

     http://magma.maths.usyd.edu.au/calc/  

(or start Magma if it is installed on your computer).

2. Paste the content of this file in there.

3. Modify the final lines by inserting your own polynomial and by correcting
   for the singular and ramified primes.

4. Click Submit

*/

// Code written by Reinier Broker for Mastermath Algebraic Number Theory 2004.

E := function(p,f)

	F := FiniteField(p);
	S<y> := PolynomialRing(F);
	g := S ! f;

	n := #Roots(g);

	if n eq 4 then 
		return (1-(1/p))^3;
	end if;
	if n eq 2 then
		return (1-(1/p))*(1-(1/p)^2);
	end if;
	if n eq 1 then
		return (1-(1/p)^3);
	end if;

	d := Discriminant(f);

	if n eq 0 and KroneckerSymbol(d,p) eq 1 then
		return (1+(1/p))*(1-(1/p)^2);
	end if;

	return 1+(1/p)+(1/p^2)+(1/p^3);
end function;
	

Ep := function(b,f)
// f should have degree 4
// ramifying primes are not considered here (modify yourself!)

	T := 1.0;

	for i in [1..b] do
		if IsPrime(i) then
			T := T*E(i,f)^(-1);
		end if;
	end for;

	return T;
end function;

P<x>:=PolynomialRing(Integers());

// Insert your polynomial f here:

f := x^4 - 2*x^2 + 3*x - 7;

// The function Ep(N, f) gives the Euler product, but the Euler factors are
// correct only for primes that are regular and unramified in ZZ[x]/(f),
// where they are given by the formula on page 75 of the notes.
// We have to correct for the other primes.

// Here are some (corrected) entries of the table on page 75, where the
// unique prime number that ramifies in ZZ[x]/(f) is 98443, which in O_K
// has exactly 3 primes of degree 1 over it.

Ep(100, f);
Ep(500, f);
Ep(2000, f);
Ep(10000, f);
Ep(20000, f);
Ep(100000, f) * E(98443, f) * (1-1/98443)^-2;

// In the last line, the incorrect Euler factor E(98443, f)^-1 is removed
// and the correct Euler factor ( (1-1/p)^3 / (1-1/p) )^-1=(1-1/p)^-2 is
// inserted.

// Good luck!
// Marco Streng, 2012