
function newton_maehly_test ( )

%*****************************************************************************80
%
%% newton_maehly_test() tests newton_maehly().
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    15 January 2026
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../newton_maehly' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'newton_maehly_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version: %s\n', version ( )  );
  fprintf ( 1, '  newton_maehly() uses the Newton Maehly algorithm\n' );
  fprintf ( 1, '  to compute all roots of a polynomial.\n' );

  newton_maehly_test01 ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'newton_maehly_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  timestamp ( );
  rmpath ( '../newton_maehly' );

  return
end
function newton_maehly_test01 ( )

%*****************************************************************************80
%
%% newton_maehly_test01() considers degree 4 polynomial with a complex pair of roots.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    15 January 2026
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'newton_maehly_test01():\n' );
  fprintf ( 1, '  Try newton_maehly() on a degree four polynomial with a pair of complex roots.\n' );

  c = [ -56.0, -75.0, -26.0, 0.0, 1.0 ];

  poly_print ( c, '  Polynomial to be analyzed:' );

  roots = newton_maehly ( c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed roots:\n' );
  disp ( roots );

  exact = [ ...
    0.5 * (   5.0 + sqrt ( 57.0 ) ), ...
    0.5 * (   5.0 - sqrt ( 57.0 ) ), ...
    0.5 * ( - 5.0 + sqrt ( -3.0 ) ), ...
    0.5 * ( - 5.0 - sqrt ( -3.0 ) ) ];
  fprintf ( 1, '  Exact roots:\n' );
  disp ( exact );

  return
end
function poly_print ( c, label )

%*****************************************************************************80
%
%% poly_print() prints a polynomial.
%
%  Discussion:
%
%    The power sum form of a polynomial is:
%
%      p(x) = c(0) + c(1) * x + ... + c(n-1) * x^(n-1) + c(n) * x^(n)
%
%    Because MATLAB doesn't allow 0 indexing, we use an adjusted formula:
%
%      p(x) = c(1) + c(2) * x + ... + c(n) * x^(n-1) + c(n+1) * x^(n)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    23 October 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      if ( 2 <= i )
        fprintf ( 1, '%+14f * x^%d\n', c(i+1), i );
      elseif ( i == 1 )
        fprintf ( 1, '%+14f * x\n', c(i+1) );
      elseif ( i == 0 )
        fprintf ( 1, '%+14f\n', c(i+1) );
      end

    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 November 2022
%
%  Author:
%
%    John Burkardt
%
  t = now ( );
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

