/*

This file defines several functions, useful when studying zeros of complex functions.

Three methods are implemented:
- linear regression
- looking at the variation of the logarithm of the function (the principle of argument)
- integration of the logarithmic derivative

Please note that the results of the methods presented here may only provide
"experimental evidence", not proofs, for the existence of non-existence
of zeros inside a given contour.

*/


/*
General helper functions.
*/

relativechange (a, b) = 
{
	local (mag);
	mag = abs(a) + abs(b);
	if (mag == 0,
		1,
		abs(a-b) / mag
		)
}

magnitude (x) =
{
	if (x == 0,
		-9999,
		round (log (abs (x)))
		);
}


/*

1. Linear regression.
Using this method one can detect and locate simple zeros, lack of zeros
close to a given point, and rule out the case of multiple zeros.

We check if the equality
(1)   f(s) = a1 (s-s0) + a0
holds (approximately) in a certain neighbourhood of s0.
If it does, the location of a zero close to s0 (if any) can be estimated from a0 and a1.

Usage:
lrdata = lrinit (s0); \\ Initialize intermediate data. Put a number or a var in place of s0.
lrdata = lraddpoint (lrdata, x, y); \\ Add a point (x, y), i.e. (x, f(x)), to the data.
lrdata = lraddpoint (lrdata, ..., ...); \\ Typically the points x will be randomly chosen close to and/or around s0.
...
lrsimplezeroreport (lrsimplezero (lrdata)) \\ Print the results, in a descriptive form.

The function lrsimplezero(lrdata) returns a vector:
 [corr^2, N, rho, mistake-likelihood]
where:
 corr^2 is the squared correlation coefficient (it should be very close to 1
        if we should accept that (1) is true),
 N      is the number of zeros (on the assumption of (1)), can be 0 or 1;
        we put N = 1 if the distance of the estimated location of a zero from
        s0 is not greater than that of the most distant of the data points,
 rho    is the estimated location of a zero if N == 1, undefined otherwise,
 mistake-likelihood reflects the uncertainty of the decision between N = 0 and
        N = 1; equals +/- magnitude(estimated distance / max data point distance);
        low negative value is good.

*/

lrinit (s0) =
{
	vector (8, X, if (X==1, s0, 0));
}

lraddpoint (corrdata, x, y) =
{
	x = x - corrdata[1];
	corrdata[2] = corrdata[2] + 1;
	corrdata[3] = corrdata[3] + x;
	corrdata[4] = corrdata[4] + y;
	corrdata[5] = corrdata[5] + norm(x);
	corrdata[6] = corrdata[6] + norm(y);
	corrdata[7] = corrdata[7] + y*conj(x);
	corrdata[8] = max (corrdata[8], norm (x));
	corrdata;
}

lrsimplezero (corrdata) =
{
	local (a0, a1, zeroinside);
	corrdata[3] = corrdata[3] / corrdata[2];
	corrdata[4] = corrdata[4] / corrdata[2];
	corrdata[5] = corrdata[5] / corrdata[2] - norm(corrdata[3]);
	corrdata[6] = corrdata[6] / corrdata[2] - norm(corrdata[4]);
	corrdata[7] = corrdata[7] / corrdata[2] - corrdata[4]*conj(corrdata[3]);
	
	a1 = corrdata[7]/corrdata[5];
	a0 = corrdata[4] - a1*corrdata[3];
	zeroinside = norm(a1)*corrdata[8]/norm(a0);
	if (zeroinside > 1,
		return ([
			norm (corrdata[7])/(corrdata[5] * corrdata[6]),
			1,
			corrdata[1] - a0/a1,
			-magnitude(zeroinside)
			])
		,
		return ([
			norm (corrdata[7])/(corrdata[5] * corrdata[6]),
			0,
			9999.99,
			magnitude(zeroinside)
			])
		);
}

lrsimplezeroreport (lrsimplezerodata) =
{
	local (mystr);
	mystr = concat (
		concat ("Linear regression (corr square ", lrsimplezerodata[1]),
		concat (", uncertainty ", lrsimplezerodata[4]));
	mystr = concat (mystr, ") ");
	if (lrsimplezerodata[2] == 1,
		mystr = concat (
			concat (mystr, "detects a simple zero at "),
			concat (lrsimplezerodata[3], "."))
		,
		mystr = concat (mystr, "detects no zero.")
		);
}


/*

2. Variation of the logarithm.
Using this method one can compute the number of zeros of a function meromorphic
inside a given contour, regular along the contour. Zeros are counted with
multiplicities, poles treated as zeros with negative multiplicities.
The value of log (f(s)) is examined for subsequent points along the closed curve.
You must supply enough points. The greatest distance between subsequent
logarithms mod 2*Pi*I is computed - it should be very close to 0 if the result
is to be believed.

Usage:
lvdata = lvinit (fs); \\ Initialize intermediate data. Specify the value at the final point along the contour.
lvdata = lvaddpoint (lvdata, y); \\ Specify the value at subsequent points along the contour.
lvdata = lvaddpoint (lvdata, ...); \\ The points must be given in order.
lvdata = lvaddpoint (lvdata, ...); \\ The contour has to be closed.
lvdata = lvaddpoint (lvdata, ...); \\ It is necessary to repeat the initial point at the end, but not at the beginning.
...
lvreport (lvresult (lvdata)); \\ Print the results, in a descriptive form.

The function lvresult(lvdata) returns a vector:
 [N, magnitude (max-distance/Pi), number of points]
where:
 N            is the number of zeros,
 max-distance is the greatest distance mod 2*Pi*I between logarithms of
              subsequent values.

*/

lvinit (fs) =
{
	[log (fs), 0, 0, 0];
}

lvaddpoint (lvdata, y) =
{
	local (dlog);
	y = log (y);
	dlog = (y - lvdata[1])/(Pi*I);
	if (real(dlog) > 1,
		dlog = dlog - 2;
		lvdata[2] = lvdata[2] - 1;
		);
	if (real(dlog) < -1,
		dlog = dlog + 2;
		lvdata[2] = lvdata[2] + 1;
		);
	lvdata[3] = max (lvdata[3], abs(dlog));
	lvdata[1] = y;
	lvdata[4] = lvdata[4] + 1;
	lvdata;
}

lvresult (lvdata) =
{
	[lvdata[2], magnitude (lvdata[3]), lvdata[4]];
}

lvreport (lvresultdata) = 
{
	concat (
		concat (
			concat ("Argument of the logarithm (", lvresultdata[3]),
			concat (" points, max change ", lvresultdata[2])
			),
		concat (
			concat (") indicates ", lvresultdata[1]),
			" zero(s)."
			)
	);
}


/*

3. Integrate the logarithmic derivative along a contour.
Using this method one can compute the number of zeros of a function meromorphic
inside a given contour, regular along the contour. Zeros are counted with
multiplicities, poles treated as zeros with negative multiplicities.
The logarithmic derivative f'(s)/f(s) is integrated over the closed curve.
You must supply enough points for the result to be precise.
The number of zeros is computed as:
 N = round (integral/2PiI).
Several indicators may help you decide if you used enough points and were not
too close to a zero or pole:
 Nerr           the fractional part of the estimate; Nerr = abs (integral/2PiI - N)
 DerivativeErr  max relative change between subsequent values of f(s) or f'(s),
                computed as change(x, y) = |x-y|/(|x|+|y|)
 ZeroProximity  test if the border is getting too close to a zero;
                equals max|f(s)| / min|f(s)|.

Usage:
ldidata = ldiinit (fs, fps); \\ Initialize with the values of f(s) and f'(s) at the final point along the curve.
ldidata = ldiaddpoint (ldidata, fs, fps, ds); \\ Specify the values of f(s) and f'(s) at subsequent points along the curve.
ldidata = ldiaddpoint (ldidata, ...); \\ The points must be given in order.
ldidata = ldiaddpoint (ldidata, ...); \\ The contour has to be closed.
...
ldireport (ldiresult (ldidata)); \\ Print the results, in a descriptive form.

The function ldiresult(ldidata) returns:
 [N, magnitude(Nerr), magnitude(DerivativeErr), magnitude(ZeroProximity)].
The first two diagnostic returns should best be negative (say: -50) and the
last one - not too high.

*/

ldiinit (fs, fps) =
{
	\\ The ldata vector contains:
	\\ [last f(s), last f'(s), integral, max rel. change of f(s) or f'(s), min |f(s)|, max |f(s)|]
	[fs, fps, 0, 0, abs(fs), abs(fs)];
}

ldiaddpoint (ldidata, fs, fps, ds) =
{
	ldidata[3] = ldidata[3] + fps*ds/fs;
	ldidata[4] = max (ldidata[4], max (relativechange (ldidata[1], fs), relativechange (ldidata[2], fps)));
	ldidata[5] = min (ldidata[5], abs(fs));
	ldidata[6] = max (ldidata[6], abs(fs));
	ldidata[1] = fs;
	ldidata[2] = fps;
	ldidata;
}

ldiresult (ldidata) =
{
	local (N);
	ldidata[3] = ldidata[3] / (2*Pi*I);
	N = round (ldidata[3]);
	[N, magnitude(ldidata[3] - N), magnitude(ldidata[4]), magnitude(ldidata[6]/ldidata[5])];
}

ldireport (ldiresultdata) =
{
	concat (concat (
		concat (
			concat ("Logarithmic derivative integral (Nerr ", ldiresultdata[2]),
			concat (", max change ", ldiresultdata[3])
			),
		concat (
			concat (", zero prox. ", ldiresultdata[4]),
			concat (") indicates ", ldiresultdata[1])
			)
		), " zero(s).");
}