function b = polynomial_root_bound ( c )

%*****************************************************************************80
%
%% bisect() carries out the bisection method.
%
%  Discussion:
%
%    The method is due to Cauchy, and the bound is called the Cauchy bound.
%
%    MATLAB made the WRONG choice in indexing polynomial coefficients
%    with c(1) multiplying the highest order term!
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 December 2023
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John D Cook,
%    Bounding complex roots by a positive root,
%    https://www.johndcook.com/blog/2023/12/02/cauchy-radius/
%    02 December 2023.
%
%  Input:
%
%    complex c(n+1): the coefficients of the polynomial.
%    c(1) multiplies x^n and c(n+1) is the constant term (multiplying x^0).
%
%  Output:
%
%    real b: a bound on the norm of all roots of the polynomial.
%
  n = length ( c ) - 1;
%
%  If the constant term is zero, the problem reduces to finding roots of p(z)/z.
%  Shift coefficients, repeatedly if necessary.
%
  while ( true )
%
%  Polynomial could be identically zero.
%
    if ( n + 1 <= 0 )
      b = 0.0;
      return
    end

    if ( c(n+1) ~= 0.0 )
      break;
    end
%
%  Drop the last coefficient, which is zero and reduce the order.
%
    n = n - 1
    c = c(1:n+1);

  end
%
%  Define associated real polynomial q(x):
%    q(x) = |c1| x^(n) - |c2| x^(n-1) - |c3| x^(n-2) - ... - |cn+1|
%
  q = zeros ( n + 1, 1 );

  q(1) = abs ( c(1) );
  for i = 2 : n + 1
    q(i) = - abs ( c(i) );
  end
%
%  Determine change of sign interval for Q(x).
%
  xneg = 0.0;
  qneg = q(n+1);

  xpos = 1.0;
  qpos = polyval ( q, xpos );
  while ( qpos <= 0.0 )
    xpos = xpos * 2.0;
    qpos = polyval ( q, xpos );
  end
%
%  Use bisection to find X such that Q(X)=0.
%
  tol = 10.0 * eps * abs ( xpos - xneg );
  [ xneg, xpos, it ] = bisect ( xneg, xpos, tol, @(x) polyval(q,x) );

  b = 0.5 * ( xneg + xpos );

  return
end
function [ a, b, it ] = bisect ( a, b, tol, f )

%*****************************************************************************80
%
%% bisect() carries out the bisection method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A, B, the endpoints of a change of sign interval.
%
%    real TOL, the interval size tolerance.  Once |B-A| < TOL,
%    the iteration will stop.
%
%    pointer value=F(x), a pointer to a function.  This can be
%    an anonymous function or the name of a MATLAB M-file function.
%
%  Output:
%
%    real A, B, the new endpoints that constitute an
%    change of sign interval no larger than TOL.
%
%    integer IT, the number of bisections.
%
  fa = f ( a );
  fb = f ( b );
  it = 0;

  while ( tol < abs ( b - a ) )
    c = ( a + b ) / 2.0;
    fc = f(c);
    it = it + 1;

    if ( 100 <= it )
      error ( 'bisect(): too many iterations!\n' );
    end

    if ( sign ( fc ) == sign ( fa ) )
      a = c;
      fa = fc;
    else
      b = c;
      fb = fc;
    end
  end

  return
end

