# include <complex.h>
# include <math.h>
# include <stdbool.h>
# include <stdlib.h>
# include <stdio.h>

# include "polynomial_root_bound.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    11 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.

    double complex 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 )
    {
      printf ( "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 = ( double * ) malloc ( ( d + 1 ) * sizeof ( double ) );

  for ( i = 0; i <= d; i++ )
  {
    if ( i < d )
    {
      q[i] = - cabs ( c[i+e] );
    }
    else
    {
      q[i] = cabs ( 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 )
    {
      printf ( "\n" );
      printf ( "polynomial_root_bound(): Warning!\n" );
      printf ( "  Too many bisection steps.\n" );
      break;
    }
    if ( fabs ( fx ) < 1.0E-10 )
    {
      printf ( "Convergence after %d iterations\n", it );
      break;
    }
  }

  value = 0.5 * ( xneg + xpos );

  free ( q );

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

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

/******************************************************************************/
/*
  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:

    11 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 the previous call.
    Not needed on the first call, with JOB = 0.

    int *JOB: a communication flag.
    The user sets JOB to 0 before the first call.

  Output:

    double *A, *B: new values for the change of sign interval.

    double *X: a new point at which the function is to be evaluated.

    int *JOB: a communication flag used by the program.
*/
{
  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 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "bisection_rc(): Fatal error\n" );
      fprintf ( stderr, "  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    11 December 2023

  Author:

    John Burkardt

  Input:

    int M, the degree of the polynomial.

    double C[0:M], the coefficients of the polynomial.
    C[0] is the constant term.

    double X, the point at which the polynomial is to be evaluated.

  Output:

    double R8POLY_VALUE, 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;
}