function roots = newton_maehly ( c, max_iter, tol )

%*****************************************************************************80
%
%% newton_maehly() applies the Newton-Maehly algorithm to a polynomial.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    15 January 2026
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    complex c[d+1]: the polynomial coefficients.  c[0] is the constant term
%    and c[d] is the coefficients of x**d.
%
%    real tol: convergence tolerance.
%
%    integer maxiter: maximum number of iterations
%
%  Output:
%
%    complex z[d]: the roots.
%
  if ( nargin < 2 )
    max_iter = 100;
  end

  if ( nargin < 3 )
    tol = 1.0e-12;
  end

  verbose = true;

  d = length ( c ) - 1;
%
%  Cauchy bound for root norms.
%
  radius = 1.0 + max ( abs ( c(1:d+1) / c(d+1) ) );
%
%  Initial guess: roots of unity scaled by R
%
  theta = linspace ( 0.0, 2.0 * pi, d + 1 );
  theta = complex ( 0.0, theta(1:d) );
  roots = radius * exp ( theta );
%
%  Repeatedly apply the mapping.
%
  if ( verbose )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Iter  Change  ||Poly||\n' );
    fprintf ( 1, '\n' );
  end

  for iter = 1 : max_iter

    roots_old = roots;

    for i = 1 : d

      [ pz, dpz ] = poly_and_derivative ( c, roots(i) );
      s = 0.0;
      for j = 1 : d
        if ( j ~= i )
          s = s + 1.0 / ( roots(i) - roots(j) );
        end
      end
      denom = dpz - pz * s;

      if ( denom == 0.0 )
        error ( "Zero denominator encountered." );
      end

      roots(i) = roots(i) - pz / denom;

    end
%
%  Convergence check
%
    [ p, dp ] = poly_and_derivative ( c, roots );

    max_change = max ( abs ( roots - roots_old ) );
    max_poly = max ( abs ( p ) );

    if ( verbose )
      fprintf ( 1, '  %4d  %11.5g  %11.5g\n', iter, max_change, max_poly );
    end

    if ( max_change < tol )
      return;
    end

  end

  error ( "Newton-Maehly did not converge." );

  return


