/*

This Pari script may help find zeros of the Dedekind zeta function.
First, we read some helper functions and define a few more.
Examples of using them are in the comments below.

*/


read("zeros-search.gp");


/*
Count the number of zeros within a rectangle a < Re s < b, c < Im s < d.

Parameters:
a, b, c, d
m, n - the number of points on the rectangle sides (horizontal and vertical, resp.)
       where the function is to be evaluated.

Returns:
[lr, lv]
where each component is the result of a corresponding method in "zeros-search.gp".
*/

ZerosCount (nfz, a, b, c, d, m, n) = 
{
	local (s0, ds, scurrent, Zs);
	local (lrdata, lvdata);

	scurrent = a + (c + (d-c)/(2*n))*I;
	Zs = zetak (nfz, scurrent);
	
	lrdata = lrinit (a + (b-a)/2 + (c + (d-c)/2)*I);
	lvdata = lvinit (Zs);
	
	ds = (b-a)/m;
	s0 = a + c*I - ds/2;
	for (i=1, m,
		scurrent = s0 + i*ds;
		Zs = zetak (nfz, scurrent);
		lrdata = lraddpoint (lrdata, scurrent, Zs);
		lvdata = lvaddpoint (lvdata, Zs);
	);

	ds = (d-c)*I/n;
	s0 = b + c*I - ds/2;
	for (i=1, n,
		scurrent = s0 + i*ds;
		Zs = zetak (nfz, scurrent);
		lrdata = lraddpoint (lrdata, scurrent, Zs);
		lvdata = lvaddpoint (lvdata, Zs);
	);

	ds = (a-b)/m;
	s0 = b + d*I - ds/2;
	for (i=1, m,
		scurrent = s0 + i*ds;
		Zs = zetak (nfz, scurrent);
		lrdata = lraddpoint (lrdata, scurrent, Zs);
		lvdata = lvaddpoint (lvdata, Zs);
	);

	ds = (c-d)*I/n;
	s0 = a + d*I - ds/2;
	for (i=1, n,
		scurrent = s0 + i*ds;
		Zs = zetak (nfz, scurrent);
		lrdata = lraddpoint (lrdata, scurrent, Zs);
		lvdata = lvaddpoint (lvdata, Zs);
	);
	[lrsimplezero(lrdata), lvresult(lvdata)];
}

NZ (nfz, a, b, n) = 
{
	local (r);
	r = (b - a)/2;
	if (r > 0.25, r = 0.25);
	ZerosCount (nfz, 0.5 - r, 0.5 + r, a, b, n, n);
}

showresult (nfz, a, b, n) = 
{
	local (results);
	results = NZ (nfz, a, b, n);
	print (lrsimplezeroreport (results[1]));
	print (lvreport (results[2]));
}



/*

Examples of computing the zeros are given below.
They are all commented out - you can try out some or all of them at your convenience.

*/



/*
fpol = x^2 - 26;
nfz = zetakinit (fpol);
showresult (nfz, 1.37058396, 1.370583969, 10000);

\\ Linear regression (corr square 0.9999999999999999839834395309, uncertainty -9) detects
\\ a simple zero at 0.5000000000000000000044752319 + 1.370583964578187003117066175*I.
\\ Argument of the logarithm (40000 points, max change -10) indicates 1 zero(s).
\\ 1.370583964578 - sig. digits
*/


/*
fpol = x^2+65;
nfz = zetakinit (fpol);
showresult (nfz, 1.05325893699922443, 1.053258936999224496, 10000);

\\ Linear regression (corr square 0.9999999999999992475, uncertainty -10) detects
\\ a simple zero at 0.5000000000000000000000000000 + 1.053258936999224463261538680*I.
\\ Argument of the logarithm (40000 points, max change -10) indicates 1 zero(s).
*/


/*
fpol = x^2 + 9982;
nfz = zetakinit (fpol);
showresult (nfz, 0.276597017487181000, 0.276597017487181176, 10000);

\\ Linear regression (corr square 0.9999999999999992190, uncertainty -13) detects
\\ a simple zero at 0.5000000000000000000000000000 + 0.2765970174871810881810812134*I.
\\ Argument of the logarithm (40000 points, max change -10) indicates 1 zero(s).
\\ ***   last result computed in 31mn, 59,680 ms (1.6 GHz Intel Pentium 4).
*/


/*
fpol = x^3 - x^2 + 7*x + 8;
nfz = zetakinit (fpol);
showresult (nfz, 1.350474195, 1.350474196123217711143088340, 10000);

\\ Linear regression (corr square 0.9999999999999999980161159056, uncertainty -86) detects
\\ a simple zero at 0.5000000000000000000000000000 + 1.350474195561608855571544170*I.
\\ Argument of the logarithm (40000 points, max change -10) indicates 1 zero(s).
*/


/*
allocatemem;
fpol = x^3 - x^2 - 97*x - 384;
nfz = zetakinit (fpol);
showresult (nfz, 0.430639281, 0.4306392814897862936621532496, 1000);

\\ Linear regression (corr square 0.9999999999999999995405959897, uncertainty -88) detects
\\ a simple zero at 0.5000000000000000000000000000 + 0.4306392812448931468310766248*I.
\\ Argument of the logarithm (4000 points, max change -7) indicates 1 zero(s).
\\ ***   last result computed in 1hr, 5mn (1.6 GHz Intel Pentium 4).
*/