# Complex Newton's Method

I'm trying to build a complex Newton's method. I've submitted a PR here, where you can see the documentation and tests, as well as get a compiling example. I'd appreciate y'alls help reducing the workload on the maintainer.

This uses the idea that if xn converges to x*, then f(x*) = 0 whenever f' is continuous at x*. If f'(xn) = 0 at some point, it resorts to Muller's method.

#ifndef BOOST_NO_CXX11_AUTO_DECLARATIONS
/*
* Why do we set the default maximum number of iterations to the number of digits in the type?
* Because for double roots, the number of digits increases linearly with the number of iterations,
* so this default should recover full precision even in this somewhat pathological case.
* For isolated roots, the problem is so rapidly convergent that this doesn't matter at all.
*/
template<class Complex, class F>
Complex complex_newton(F g, Complex guess, int max_iterations=std::numeric_limits<typename Complex::value_type>::digits)
{
typedef typename Complex::value_type Real;
using std::norm;
using std::abs;
using std::max;
// z0, z1, and z2 cannot be the same, in case we immediately need to resort to Muller's Method:
Complex z0 = guess + Complex(1,0);
Complex z1 = guess + Complex(0,1);
Complex z2 = guess;

do {
auto pair = g(z2);
if (norm(pair.second) == 0)
{
// Muller's method. Notation follows Numerical Recipes, 9.5.2:
Complex q = (z2 - z1)/(z1 - z0);
auto P0 = g(z0);
auto P1 = g(z1);
Complex qp1 = static_cast<Complex>(1)+q;
Complex A = q*(pair.first - qp1*P1.first + q*P0.first);

Complex B = (static_cast<Complex>(2)*q+static_cast<Complex>(1))*pair.first - qp1*qp1*P1.first +q*q*P0.first;
Complex C = qp1*pair.first;
Complex rad = sqrt(B*B - static_cast<Complex>(4)*A*C);
Complex denom1 = B + rad;
Complex denom2 = B - rad;
Complex correction = (z1-z2)*static_cast<Complex>(2)*C;
if (norm(denom1) > norm(denom2))
{
correction /= denom1;
}
else
{
correction /= denom2;
}

z0 = z1;
z1 = z2;
z2 = z2 + correction;
}
else
{
z0 = z1;
z1 = z2;
z2 = z2  - (pair.first/pair.second);
}

// See: https://math.stackexchange.com/questions/3017766/constructing-newton-iteration-converging-to-non-root
// If f' is continuous, then convergence of x_n -> x* implies f(x*) = 0.
// This condition approximates this convergence condition by requiring three consecutive iterates to be clustered.
Real tol = max(abs(z2)*std::numeric_limits<Real>::epsilon(), std::numeric_limits<Real>::epsilon());
bool real_close = abs(z0.real() - z1.real()) < tol && abs(z0.real() - z2.real()) < tol && abs(z1.real() - z2.real()) < tol;
bool imag_close = abs(z0.imag() - z1.imag()) < tol && abs(z0.imag() - z2.imag()) < tol && abs(z1.imag() - z2.imag()) < tol;
if (real_close && imag_close)
{
return z2;
}

} while(max_iterations--);

// The idea is that if we can get abs(f) < eps, we should, but if we go through all these iterations
// and abs(f) < sqrt(eps), then roundoff error simply does not allow that we can evaluate f to < eps
// This is somewhat awkward as it isn't scale invariant, but using the Daubechies coefficient example code,
// I found this condition generates correct roots, whereas the scale invariant condition discussed here:
// https://scicomp.stackexchange.com/questions/30597/defining-a-condition-number-and-termination-criteria-for-newtons-method
// allows nonroots to be passed off as roots.
auto pair = g(z2);
if (abs(pair.first) < sqrt(std::numeric_limits<Real>::epsilon()))
{
return z2;
}

return {std::numeric_limits<Real>::quiet_NaN(),
std::numeric_limits<Real>::quiet_NaN()};
}
#endif


Usage:

boost::math::tools::polynomial<std::complex<double>> p{{1,0}, {0, 0}, {1,0}};
std::complex<double> guess{1,1};
boost::math::tools::polynomial<std::complex<double>> p_prime = p.prime();
auto f = [&](std::complex<double> z) { return std::make_pair<std::complex<double>, std::complex<double>>(p(z), p_prime(z)); };
std::complex<double> root = complex_newton(f, guess);


For a more intense example usage, see here

I can't see too much that needs improving. I'd suggest that the correction logic be simplified to:

       Complex denom1 = B + rad,
denom2 = B - rad,
correction = (z1-z2)*static_cast<Complex>(2)*C,
denom_corr;
if (norm(denom1) > norm(denom2))
denom_corr = denom1;
else
denom_corr = denom2;

z0 = z1;
z1 = z2;
z2 += correction/denom_corr;


Elsewhere, combine your types and use commas where there are multiple variable declarations. It's a stylistic choice, but I don't like to repeat myself.

And this:

z2 = z2  - (pair.first/pair.second);


should be

z2 -= pair.first/pair.second;


Order of operations still applies so you don't need parens, and use combined operation+assignment.