function roots = wdk ( c, tol, max_iter )

%*****************************************************************************80
%
%% wdk() applies the Weierstrass-Durand-Kerner 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 )
    tol = 1.0e-12;
  end

  if ( nargin < 3 )
    max_iter = 1000;
  end

  verbose = true;

  d = length ( c ) - 1;
%
%  Cauchy bound for initial radius
%
  R = 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 = R * exp ( theta );
%
%  Repeatedly apply the mapping.
%
  if ( verbose )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Iter  Change  ||Poly(root)||\n' );
    fprintf ( 1, '\n' );
  end

  for iter = 1 : max_iter

    roots_old = roots;

    for i = 1 : d
      zi = roots_old(i);
      denom = 1.0;
      for j = 1 : d
        if ( i ~= j )
          denom = denom * ( zi - roots(j) );
        end
      end

      roots(i) = zi - poly_eval ( c, zi ) / denom;

    end

    max_change = max ( abs ( roots - roots_old ) );
    max_poly = max ( abs ( poly_eval ( c, roots ) ) );
%
%  Optional report.
%
    if ( verbose )
      fprintf ( 1, '  %4d  %11.5g  %11.5g\n', iter, max_change, max_poly );
    end
%
%  Exit the loop if the roots haven't changed much.
%
    if ( max_change < tol )
      break;
    end

  end

  return
end

