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

# include "polynomial_root_bound.h"

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

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

int main ( )

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

    polynomial_root_bound_test() tests polynomial_root_bound().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 December 2023

  Author:

    John Burkardt
*/
{
  timestamp ( );
  printf ( "\n" );
  printf ( "polynomial_root_bound_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test polynomial_root_bound()\n" );

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

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

void polynomial1_test ( )

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

    polynomial1_test() deals with a particular polynomial.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 December 2023

  Author:

    John Burkardt
*/
{
  double b;
  double complex c[6] = { 23.0 * I, 0.0, 2.0, 0.0, 0.0, 12.0 };
  int i;
  int n = 5;
  double complex r[5] = {
    -1.10401598 - 0.33686293 * I,
    -0.66152059 + 0.89701442 * I,
     0.02570877 - 1.13896251 * I,
     0.67740946 + 0.94586564 * I,
     1.06241834 - 0.36705461 * I };

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

  b = polynomial_root_bound ( n, c );

  printf ( "\n" );
  printf ( "  Root magnitude bound is %g\n", b );

  printf ( "\n" );
  printf ( "  Polynomial roots and norms are:\n" );
  for ( i = 0; i < n; i++ )
  {
    printf ( "  ( %g, %g ), %g\n", 
      creal ( r[i] ), cimag ( r[i] ), cabs ( r[i] ) );
  }

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

void polynomial2_test ( )

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

    polynomial2_test() deals with a random polynomial.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 December 2023

  Author:

    John Burkardt
*/
{
  double b;
  double complex *c;
  int i;
  int n = 5;
  double complex *r;

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

  r = c8vec_normal_01_new ( n );
  printf ( "\n" );
  printf ( "  Polynomial roots and norms are:\n" );
  for ( i = 0; i < n; i++ )
  {
    printf ( "  ( %g, %g ), %g\n", 
      creal ( r[i] ), cimag ( r[i] ), cabs ( r[i] ) );
  }

  c = roots_to_c8poly ( n, r );

  printf ( "\n" );
  printf ( "  Polynomial coefficients are:\n" );
  for ( i = 0; i <= n; i++ )
  {
    printf ( "  ( %g, %g )\n", creal ( c[i] ), cimag ( c[i] ) );
  }

  b = polynomial_root_bound ( n, c );

  printf ( "\n" );
  printf ( "  Root magnitude bound is %g\n", b );

  free ( c );
  free ( r );

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

double complex *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:

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

  x = ( double complex * ) malloc ( n * sizeof ( double complex ) );

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

    x[i] = sqrt ( - 2.0 * log ( v1 ) ) *
      cos ( 2.0 * r8_pi * v2 ) + sin ( 2.0 * r8_pi * v2 ) * I;
  }

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

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

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

    roots_to_c8poly() converts polynomial roots to polynomial coefficients.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    07 December 2023

  Author:

    John Burkardt

  Input:

    int N, the number of roots.

    double complex R[N], the roots.

  Output:

    double complex ROOTS_TO_C8POLY[N+1], the coefficients of the polynomial.
*/
{
  double complex *c;
  int i;
  int j;

  c = ( double complex * ) malloc ( ( n + 1 ) * sizeof ( double complex ) );
/*
  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;
}
/******************************************************************************/

void timestamp ( )

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

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

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

  now = time ( NULL );
  tm = localtime ( &now );

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

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}