/*

Using this script you may be able to prove that a Hecke zeta function has a simple zero
close to a given point or that it has no zero there.
Only Hecke zeta functions associated to characters of the ordinary ideal class group H(K)
of an algebraic number field K are considered.

We make use of the some effective estimates for the derivatives of the functions in question.
Let d_K denote the absolute value of the field's discriminant,
n the degree, s = sigma + it the complex variable.
Then we have

|d/ds (s-1)zeta_K(s)| <= 4/3 max(d_K pi^{-n}, 2^n) (|t|+3)^{n+1},
|d^2/ds^s (s-1)zeta_K(s)| <= 4 max(d_K pi^{-n}, 2^n) (|t|+3)^{n+1}, 0.25 <= sigma <= 0.75,

and, for all non-principal characters chi of H(K),

|d/ds zeta(s, chi)| <= 4/3 max(d_K pi^{-n}, 2^n) (|t|+3)^n,
|d^2/ds^s zeta(s, chi)| <= 4 max(d_K pi^{-n}, 2^n) (|t|+3)^n, 0.25 <= sigma <= 0.75.

We evaluate a function in question and its derivatives using the ComputeL package
and then check the criteria for existence or non-existence of a zero.
The location of a possible zero must be known in advance.

Script parameters:
fpol	-	generating polynomial of the field K
chi		-	a character of H(K), given as a vector coordinates in the cyclic decomposition
ImRho   -   imaginary part of suspected zero

You can also set chi=0 to specify the principal character whithout computing
the class group.


This script is suitable for batch processing.
To use it, first set the parameters, like this:

fpol = x^2+65;
chi = [1,0];
ImRho = 1.05325893699922446326153868;

Then read this file into Pari and, finally, compute results, as in the example
in the comments at the bottom.

*/






/*

Main functions.

*/

/*

Try to prove that zeta_K(s) = 0 for an s close to rho.
Parameters:
 zkr, zkpr - zeta_K(rho), zeta_K'(rho)
 t - the imaginary part of the suspected zero; rho = 1/2 + it
 n - the degree of the field
 C - max(d_K pi^{-n}, 2^n)

The return value is a vector:
 [10, anything] means we could not conclude anything
 [x, r] means: if x < 1 then there is a simple zero
        within the distance r from rho.

*/

ProveZeroZetaK (zkr, zkpr, t, n, C) = 
{
	local (fr, fpr, fppBound, r, x);
	fr = abs(zkr*(-0.5+t*I)); \\ The absolute value of f(s) = (s-1)zeta_K(s) at rho.
	fpr = abs(zkpr*(-0.5+t*I) + zkr); \\ |f'(rho)|
	fppBound = 4 * C * (abs(t) + 3.25)^(n+1);  \\ Bound for f''(s) within 0.25 from rho.

	if (fpr == 0,
		return ([10, 0])
		);
	r = 2*fr/fpr;
	if (r > 0.25,
		return ([10, 0])
		);
	x = 2*fr*fppBound/(fpr*fpr);
	if (x >= 1,
		return ([10, 0])
		);
	return ([x, r]);
}

MinZeroDistanceZetaK (zkr, t, n, C) = 
{
	local (fr, fpBound);
	fr = abs(zkr*(-0.5+t*I)); \\ The absolute value of f(s) = (s-1)zeta_K(s) at rho.
	fpBound = 4/3 * C * (abs(t) + 3.25)^(n+1); \\ Bound for f'(s) within 0.25 from rho.
	return (fr/fpBound);
}

ProveZeroZetaChi (zcr, zcpr, t, n, C) = 
{
	local (fr, fpr, fppBound, r, x);
	fr = abs(zcr); \\ The absolute value of zeta(rho, chi).
	fpr = abs(zcpr); \\ |zeta'(rho, chi)|
	fppBound = 4 * C * (abs(t) + 3.25)^n;  \\ Bound for f''(s) within 0.25 from rho.

	if (fpr == 0,
		return ([10, 0])
		);
	r = 2*fr/fpr;
	if (r > 0.25,
		return ([10, 0])
		);
	x = 2*fr*fppBound/(fpr*fpr);
	if (x >= 1,
		return ([10, 0])
		);
	return ([x, r]);
}

MinZeroDistanceZetaChi (zcr, t, n, C) = 
{
	local (fr, fpBound);
	fr = abs(zcr); \\ The absolute value of zeta(rho, chi).
	fpBound = 4/3 * C * (abs(t) + 3.25)^n; \\ Bound for zeta'(s, chi) within 0.25 from rho.
	return (fr/fpBound);
}


/*

Helper functions to output results.

*/

NiceReport (chi, bnf, zr, zpr, t) =
{
	local (n, C, zeroexists, zerodistance);
	print ("rho = 0.5 + ", ImRho*I);
	n = poldegree (bnf.pol);
	C = max (abs(bnf.disc)/Pi^n, 2^n);
	if (chi == vector (matsize (bnf.cyc)[2]),
		zeroexists = ProveZeroZetaK (zr, zpr, t, n, C);
		zerodistance = MinZeroDistanceZetaK (zr, t, n, C);
	,
		zeroexists = ProveZeroZetaChi (zr, zpr, t, n, C);
		zerodistance = MinZeroDistanceZetaChi (zr, t, n, C);
	);
	if (zeroexists[1] > 1,
		print ("There is no zero within ", zerodistance, " from rho.");
		return()
		);
	if (zeroexists[2] >= zerodistance,
		print ("The criteria are met with ", zeroexists[1], " < 1.");
		print ("There is a simple zero in a distance ", zerodistance);
		print ("to ", zeroexists[2], " from rho.");
		return()
		);
	print ("Apparently the criteria for zero-existence are met: ", zeroexists[1], " < 1.");
	print ("However, the zero distances are in contradiction:");
	print (zerodistance, " < ", zeroexists[2]);
	print1 ("It should never be possible - either the estimates must be wrong or there was a computational error. ");
	print ("Please report the problem to the author.");
}

ReportZetaKZero (filename, chi, bnf, zr, zpr, t) =
{
	local (n, C, zeroexists, zerodistance);
	n = poldegree (bnf.pol);
	C = max (abs(bnf.disc)/Pi^n, 2^n);
	if (chi != vector (matsize (bnf.cyc)[2]),
		write (filename, "Error: the character is not principal!");
		return();
		);
	zeroexists = ProveZeroZetaK (zr, zpr, t, n, C);
	zerodistance = MinZeroDistanceZetaK (zr, t, n, C);
	write (filename, "zerocond = ", zeroexists[1], "; \\\\ Should be <1. Close to 0 if we are close to a zero of zeta_K(s).");
	write (filename, "zerowithin = ", zeroexists[2], "; \\\\ The upper bound for the distance of a simple zero of zetak.");
	write (filename, "contradictioncheck = ", zerodistance/zeroexists[2], "; \\\\ Min dist./max dist. Must be <1.");
}

ReportZetaChiNonZero (filename, chi, bnf, zr, t, mindistance) =
{
	local (n, C, zerodistance);
	n = poldegree (bnf.pol);
	C = max (abs(bnf.disc)/Pi^n, 2^n);
	if (chi == vector (matsize (bnf.cyc)[2]),
		write (filename, "Error: the character is principal!");
		return();
		);
	zerodistance = min (mindistance, MinZeroDistanceZetaChi (zr, t, n, C));
	write (filename, "zerodistance = ", zerodistance, "; \\\\ Lower bound for the distance to the nearest zero.");
}



\\ Initialization...

read("computel");
read("zetachi.gp");
bnf     = bnfinit(fpol);
Fdegree = poldegree(fpol);
disc    = bnf.disc;
r1      = bnf.sign[1];
r2      = bnf.sign[2];
gammaV  = concat(vector(r1+r2,X,0),vector(r2,X,1));
conductor = abs(disc);                   \\ exponential factor
weight    = 1;                           \\ L(scurrent)=sgn*L(weight-scurrent)
differentclass = bnfisprincipal (bnf, bnf.diff, 2);   \\ The class of the different.

if (chi==0,	chi = vector (matsize (bnf.cyc)[2]));

sgn       = conj (chareval (chi, differentclass, bnf.cyc));    \\ sign in the functional equation
Lpoles    = if (chi == vector (matsize (bnf.cyc)[2]), [1], []);  \\   No pole if chi is not chi0.
dzc       = dirzetachi (bnf, chi, cflength());     \\ coefficients a(k) in L(scurrent)
initLdata("dzc[k]",,"conj(dzc[k])");         \\ initialize L-series 


\\ Example of results output. To try it out, remove the comment delimiters below
\\ and define parameters at the beginning.

/*
print ("");
print ("fpol = ", fpol);
print ("chi = ", chi);
print ("func. eq. discrepancy: ", errprint(checkfeq()));
print ("");

rho = 0.5 + ImRho*I;
NiceReport (chi, bnf, L(rho), L(rho,,1), ImRho);
*/