function roots = aberth ( c, max_iter, tol )

%*****************************************************************************80
%
%% aberth() applies the Aberth-Ehrlich algorithm for polynomial roots.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    12 January 2026
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    complex c[d+1]: the polynomial coefficients.  c[0] is the constant.
%
%  Output:
%
%    complex roots[d]: the computed roots.
%
  if ( nargin < 2 )
    max_iter = 100;
  end

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

  verbose = true;

  d = length ( c ) - 1;
%
%  Compute Cauchy bound for maximum root norm.
%
  radius = 1.0 + max ( abs ( c(1:d+1) / c(d+1) ) );
%
%  First estimate is scaled roots of unity.
%
  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

      zi = roots(i);

      [ p, dp ] = poly_and_derivative ( c, zi );

      if ( abs ( p ) < tol )
        continue;
      end
%
%  Aberth correction
%
      sum_term = 0.0;
      for j = 1 : d
        if ( j ~= i )
          sum_term = sum_term + 1.0 / ( zi - roots(j) );
        end
      end

      delta = p / dp;
      roots(i) = zi - delta / ( 1.0 - delta * sum_term );

    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 )
      break;
    end

  end

  return 
end

