# include <cmath>
# include <complex>
# include <cstdlib>
# include <iomanip>
# include <iostream>

using namespace std;

# include "polynomial_root_bound.hpp"

//****************************************************************************80

double polynomial_root_bound ( int n, complex <double> c[] )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_root_bound() bounds the roots of a polynomial.
//
//  Discussion:
//
//    The method is due to Cauchy, and the bound is called the Cauchy bound.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 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:
//
//    int n: the degree of the polynomial.
//
//    complex <double> c(0:n): the coefficients of the polynomial.
//    c(n) multiplies x^n and c(0) is the constant term (multiplying x^0).
//
//  Output:
//
//    real b: a bound on the norm of all roots of the polynomial.
//
{
  int d;
  int e;
  double fx;
  int i;
  int it;
  int job;
  double *q;
  double qpos;
  double value;
  double x;
  double xneg;
  double xpos;
//
//  Determine degree d of the polynomial.
//  If the constant term is zero, the problem reduces to finding roots of p(z)/z.
//
  d = n;
  e = 0;

  while ( true )
  {
//
//  Polynomial could be identically zero.
//
    if ( d + 1 <= 0 )
    {
      cout << "polynomial_root_bound(): Polynomial is identically 0.\n";
      value = 0.0;
      return value;
    }

    if ( c[e] != 0.0 )
    {
      break;
    }
//
//  Drop the first coefficient, which is zero and reduce the degree.
//
    d = d - 1;
    e = e + 1;
  }
//
//  Define associated real polynomial q(x):
//    q(x) = |cn| x^(n) - |cn-1| x^(n-1) - |cn-2| x^(n-2) - ... - |c0|
//
  q = new double[d+1];

  for ( i = 0; i <= d; i++ )
  {
    if ( i < d )
    {
      q[i] = - abs ( c[i+e] );
    }
    else
    {
      q[i] = abs ( c[i+e] );
    }
  }
//
//  Determine change of sign interval for Q(x).
//
  xneg = 0.0;

  xpos = 1.0;
  qpos = polyval ( d, q, xpos );
  while ( qpos <= 0.0 )
  {
    xpos = xpos * 2.0;
    qpos = polyval ( d, q, xpos );
  }
//
//  Use bisection to find X such that Q(X)=0.
//
  job = 0;
  it = 0;
  while ( true )
  {
    bisection_rc ( xneg, xpos, x, fx, job );
    fx = polyval ( d, q, x );
    it = it + 1;
    if ( 100 < it )
    {
      cout << "\n";
      cout << "polynomial_root_bound(): Warning!\n";
      cout << "  Too many bisection steps.\n";
      break;
    }
    if ( fabs ( fx ) < 1.0E-10 )
    {
      cout << "Convergence after " << it << " iterations\n";
      break;
    }
  }

  value = 0.5 * ( xneg + xpos );

  delete [] q;

  return value;
}
//****************************************************************************80

void bisection_rc ( double &a, double &b, double &x, double fx, int &job )

//****************************************************************************80
//
//  Purpose:
//
//    bisection_rc() seeks a zero of f(x) in a change of sign interval.
//
//  Discussion:
//
//    The bisection method is used.
//
//    This routine uses reverse communication, so that the function is always
//    evaluated in the calling program.
//
//    On the first call, the user sets JOB = 0, and the values of A and B.
//    Thereafter, the user checks the returned value of JOB and follows 
//    directions.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 December 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double &A, &B, the endpoints of the change of sign interval.
//
//    double FX, the function value at the point X returned on previous call.
//    On the first call, FX is ignored.
//
//    int &JOB, a communication flag.
//    The user sets JOB to 0 before the first call.
//
//    Output, double BISECTION_RC, a point X at which the function is to 
//    be evaluated.
//
//  Output:
//
//    double &A, &B, the updated endpoints of the change of  sign interval.
//
//    double &X, a point at which the function is to be evaluated for next call.
//
//    int &JOB, an updated communication flag.
//
{
  static double fa;
  static double fb;
  static int state;

  if ( job == 0 )
  {
    fa = 0.0;
    fb = 0.0;
    state = 1;
    x = a;
    job = 1;
  }
//
//  Accept value of f(a).
//  Request value of f(b).
//
  else if ( state == 1 )
  {
    if ( fx == 0.0 )
    {
      b = a;
      state = 3;
      return;
    }
    fa = fx;
    x = b;
    state = 2;
  }
//
//  Accept value of f(b).
//  Request value of f at midpoint.
//
  else if ( state == 2 )
  {
    if ( fx == 0.0 )
    {
      a = b;
      state = 3;
      return;
    }
    fb = fx;

    if ( 0.0 < fa * fb )
    {
      cerr << "\n";
      cerr << "bisection_rc(): Fatal error!\n";
      cerr << "  F(A) and F(B) have the same sign.\n";
      exit ( 1 );
    }

    x = ( a + b ) / 2.0;
    state = 3;
  }
//
//  Use sign of function value to bisect current interval.
//  Determine new midpoint and request function value.
//
  else
  {
    if ( 0.0 < fa * fx )
    {
      a = x;
      fa = fx;
    }
    else
    {
      b = x;
      fb = fx;
    }
    x = ( a + b ) / 2.0;
    state = 3;
  }

  return;
}
//****************************************************************************80

double polyval ( int m, double c[], double x )

//****************************************************************************80
//
//  Purpose:
//
//    polyval() evaluates a polynomial using a naive method.
//
//  Discussion:
//
//    The polynomial 
//
//      p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m
//
//    is to be evaluated at the value X.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    22 August 2017
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the degree of the polynomial.
//
//    double C[M+1], the coefficients of the polynomial.
//    A[0] is the constant term.
//
//    double X, the point at which the polynomial is to be evaluated.
//
//  Output:
//
//    double POLYVAL, the value of the polynomial at X.
//
{
  int i;
  double value;
  double xi;

  value = c[0];
  xi = 1.0;

  for ( i = 1; i <= m; i++ )
  {
    xi = xi * x;
    value = value + c[i] * xi;
  }

  return value;
}