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

using namespace std;

# include "polynomial_root_bound.hpp"

int main ( );
void polynomial1_test ( );
void polynomial2_test ( );
complex <double> *c8vec_normal_01_new ( int n );
complex <double> *roots_to_c8poly ( int n, complex <double> r[] );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_root_bound_test() tests polynomial_root_bound().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 December 2023
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "polynomial_root_bound_test():\n";
  cout << "  C++ version\n";
  cout << "  Test polynomial_root_bound()\n";

  polynomial1_test ( );
  polynomial2_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "polynomial_root_bound_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void polynomial1_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial1_test() deals with a particular polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 December 2023
//
//  Author:
//
//    John Burkardt
//
{
  double b;
  complex <double> c[6] = { 
    complex <double> (  0.0, 23.0 ),
    complex <double> (  0.0,  0.0 ),
    complex <double> (  2.0,  0.0 ),
    complex <double> (  0.0,  0.0 ),
    complex <double> (  0.0,  0.0 ),
    complex <double> ( 12.0,  0.0 ) };
  int i;
  int n = 5;
  complex <double> r[5] = {
    complex <double> ( -1.10401598, - 0.33686293 ),
    complex <double> ( -0.66152059, + 0.89701442 ),
    complex <double> (  0.02570877, - 1.13896251 ),
    complex <double> (  0.67740946, + 0.94586564 ),
    complex <double> (  1.06241834, - 0.36705461 ) };

  cout << "\n";
  cout << "polynomial1_test():\n";
  cout << "  Bound the roots of a specific polynomial:\n";
  cout << "    12z^5 + 2z^2 + 23i.\n";
 
  cout << "\n";
  cout << "  Polynomial coefficients are:\n";
  for ( i = 0; i <= n; i++ )
  {
    cout << "  ( " 
         << real ( c[i] )
         << ", "
         << imag ( c[i] ) 
         << " )\n";
  }

  b = polynomial_root_bound ( n, c );

  cout << "\n";
  cout << "  Root magnitude bound is " << b << "\n";

  cout << "\n";
  cout << "  Polynomial roots and norms are:\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  ( " 
         << real ( r[i] )
         << ", "
         << imag ( r[i] ) 
         << " )  "
         << abs ( r[i] )
         << "\n";
  }

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

void polynomial2_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial2_test() deals with a random polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 December 2023
//
//  Author:
//
//    John Burkardt
//
{
  double b;
  complex <double> *c;
  int i;
  int n = 5;
  complex <double> *r;

  cout << "\n";
  cout << "polynomial2_test():\n";
  cout << "  Bound the roots of a random polynomial:\n";

  r = c8vec_normal_01_new ( n );
  cout << "\n";
  cout << "  Polynomial roots and norms are:\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  ( " 
         << real ( r[i] )
         << ", "
         << imag ( r[i] ) 
         << " )  "
         << abs ( r[i] )
         << "\n";
  }

  c = roots_to_c8poly ( n, r );

  cout << "\n";
  cout << "  Polynomial coefficients are:\n";
  for ( i = 0; i <= n; i++ )
  {
    cout << "  ( " 
         << real ( c[i] )
         << ", "
         << imag ( c[i] ) 
         << " )\n";
  }

  b = polynomial_root_bound ( n, c );

  cout << "\n";
  cout << "  Root magnitude bound is " << b << "\n";

  delete [] c;
  delete [] r;

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

complex <double> *c8vec_normal_01_new ( int n )

/******************************************************************************/
/*
  Purpose:

    c8vec_normal_01_new() returns a unit pseudonormal C8VEC.

  Discussion:

    The standard normal probability distribution function (PDF) has
    mean 0 and standard deviation 1.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    07 December 2023

  Author:

    John Burkardt

  Input:

    int N, the number of values desired.

  Output:

    complex <double> C8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF.
*/
{
  int i;
  const double r8_pi = 3.141592653589793;
  double v1;
  double v2;
  complex <double> *x;

  x = new complex <double>[n];

  for ( i = 0; i < n; i++ )
  {
    v1 = drand48 ( );
    v2 = drand48 ( );

    x[i] = sqrt ( - 2.0 * log ( v1 ) ) * complex <double>
      ( cos ( 2.0 * r8_pi * v2 ), sin ( 2.0 * r8_pi * v2 ) );
  }

  return x;
}
//*****************************************************************************/

complex <double> *roots_to_c8poly ( int n, complex <double> r[] )

//*****************************************************************************/
//
//  Purpose:
//
//    roots_to_c8poly() converts polynomial roots to polynomial coefficients.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 December 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of roots.
//
//    complex <double> R[N], the roots.
//
//  Output:
//
//    complex <double> ROOTS_TO_C8POLY[N+1], the coefficients of the polynomial.
//
{
  complex <double> *c;
  int i;
  int j;

  c = new complex <double>[n+1];
//
//  Initialize C to (0, 0, ..., 0, 1).
//  Essentially, we are setting up a divided difference table.
//
  for ( i = 0; i < n; i++ )
  {
    c[i] = 0.0;
  }
  c[n] = 1.0;
//
//  Convert to standard polynomial form by shifting the abscissas
//  of the divided difference table to 0.
//
  for ( j = 1; j <= n; j++ )
  {
    for ( i = 1; i <= n + 1 - j; i++ )
    {
      c[n-i] = c[n-i] - r[n+1-i-j] * c[n-i+1];
    }
  }
  return c;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    19 March 2018
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}