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

using namespace std;

# include "companion_matrix.hpp"
# include "eigs.hpp"
# include <gsl/gsl_poly.h>

int main ( );
void companion_chebyshev_test ( );
void companion_gegenbauer_test ( );
void companion_hermite_test ( );
void companion_laguerre_test ( );
void companion_legendre_test ( );
void companion_monomial_test ( );
double *chebyshev_to_monomial_matrix ( int n );
double *gegenbauer_to_monomial_matrix ( int n, double alpha );
double *hermite_to_monomial_matrix ( int n );
double *laguerre_to_monomial_matrix ( int n );
double *legendre_to_monomial_matrix ( int n );
void polynomial_chebyshev_print ( int d, double p[], string label );
void polynomial_gegenbauer_print ( int d, double p[], string label );
void polynomial_hermite_print ( int d, double p[], string label );
void polynomial_laguerre_print ( int d, double p[], string label );
void polynomial_legendre_print ( int d, double p[], string label );
void polynomial_monomial_print ( int d, double p[], string label );
void c8vec_print ( int n, complex <double> a[], string title );
double *r8mat_mv_new ( int m, int n, double a[], double x[] );
void r8mat_print ( int m, int n, double a[], string title );
void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title );
complex <double> *r8poly_roots ( int n, double p[] );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_matrix_test() tests companion_matrix().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "companion_matrix_test():\n";
  cout << "  C++ version \n";
  cout << "  companion_matrix() computes the companion matrix\n";
  cout << "  of a polynomial, in various bases.\n";

  companion_chebyshev_test ( );
  companion_gegenbauer_test ( );
  companion_hermite_test ( );
  companion_laguerre_test ( );
  companion_legendre_test ( );
  companion_monomial_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "companion_matrix_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void companion_chebyshev_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_chebyshev_test() tests companion_chebyshev().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *A;
  double *C;
  int n = 5;
  double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 };
  double *q;
  complex <double> *r;
  complex <double> *s;

  cout << "\n";
  cout << "companion_chebyshev_test():\n";
  cout << "  companion_chebyshev() computes the companion matrix\n";
  cout << "  of a polynomial p(x) in the Chebyshev basis.\n";

  polynomial_chebyshev_print ( n, p, "  p(x)" );
//
//  Compute monomial form so we can get roots directly.
//
  A = chebyshev_to_monomial_matrix ( n + 1 );
  q = r8mat_mv_new ( n + 1, n + 1, A, p );
  polynomial_monomial_print ( n, q, "  Monomial q(x) =" );

  r = r8poly_roots ( n, q );
  c8vec_print ( n, r, "  Roots of q(x):" );

  C = companion_chebyshev ( n, p );
  r8mat_print ( n, n, C, "  Chebyshev companion matrix C(p):" );

  s = eigs ( n, C );
  c8vec_print ( n, s, "  Eigenvalues of C(p):" );

  delete [] A;
  delete [] C;
  delete [] q;
  delete [] r;
  delete [] s;

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

void companion_gegenbauer_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_gegenbauer_test() tests companion_gegenbauer().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *A;
  double alpha;
  double *C;
  int i;
  int n = 5;
  double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 };
  double *q;
  complex <double> *r;
  complex <double> *s;

  cout << "\n";
  cout << "companion_gegenbauer_test():\n";
  cout << "  companion_gegenbauer() computes the companion matrix\n";
  cout << "  of a polynomial p(x) in the Gegenbauer basis.\n";

  polynomial_chebyshev_print ( n, p, "  p(x)" );
//
//  Try several parameter values.
//
  for ( i = 1; i <= 2; i++ )
  {
    alpha = i / 2.0;

    cout << "\n";
    cout << "  Gegenbauer parameter alpha will be " << alpha << "\n";
//
//   Compute monomial form so we can get roots directly.
//
    A = gegenbauer_to_monomial_matrix ( n + 1, alpha );
    q = r8mat_mv_new ( n + 1, n + 1, A, p );
    polynomial_monomial_print ( n, q, "  Monomial q(x) =" );

    r = r8poly_roots ( n, q );
    c8vec_print ( n, r, "  Roots of q(x):" );

    C = companion_gegenbauer ( n, p, alpha );
    r8mat_print ( n, n, C, "  Gegenbauer companion matrix C(p):" );

    s = eigs ( n, C );
    c8vec_print ( n, s, "  Eigenvalues of C(p):" );

    delete [] A;
    delete [] C;
    delete [] q;
    delete [] r;
    delete [] s;
  }

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

void companion_hermite_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_hermite_test() tests companion_hermite().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *A;
  double *C;
  int n = 5;
  double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 };
  double *q;
  complex <double> *r;
  complex <double> *s;

  cout << "\n";
  cout << "companion_hermite_test():\n";
  cout << "  companion_hermite() computes the companion matrix\n";
  cout << "  of a polynomial p(x) in the Hermite basis.\n";

  polynomial_hermite_print ( n, p, "  p(x)" );
//
//  Compute monomial form so we can get roots directly.
//
  A = hermite_to_monomial_matrix ( n + 1 );
  q = r8mat_mv_new ( n + 1, n + 1, A, p );
  polynomial_monomial_print ( n, q, "  Monomial q(x) =" );

  r = r8poly_roots ( n, q );
  c8vec_print ( n, r, "  Roots of q(x):" );

  C = companion_hermite ( n, p );
  r8mat_print ( n, n, C, "  Hermite companion matrix C(p):" );

  s = eigs ( n, C );
  c8vec_print ( n, s, "  Eigenvalues of C(p):" );

  delete [] A;
  delete [] C;
  delete [] q;
  delete [] r;
  delete [] s;

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

void companion_laguerre_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_laguerre_test() tests companion_laguerre().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *A;
  double *C;
  int n = 5;
  double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 };
  double *q;
  complex <double> *r;
  complex <double> *s;

  cout << "\n";
  cout << "companion_laguerre_test():\n";
  cout << "  companion_laguerre() computes the companion matrix\n";
  cout << "  of a polynomial p(x) in the Laguerre basis.\n";

  polynomial_laguerre_print ( n, p, "  p(x)" );
//
//  Compute monomial form so we can get roots directly.
//
  A = laguerre_to_monomial_matrix ( n + 1 );
  q = r8mat_mv_new ( n + 1, n + 1, A, p );
  polynomial_monomial_print ( n, q, "  Monomial q(x) =" );

  r = r8poly_roots ( n, q );
  c8vec_print ( n, r, "  Roots of q(x):" );

  C = companion_laguerre ( n, p );
  r8mat_print ( n, n, C, "  Laguerre companion matrix C(p):" );

  s = eigs ( n, C );
  c8vec_print ( n, s, "  Eigenvalues of C(p):" );

  delete [] A;
  delete [] C;
  delete [] q;
  delete [] r;
  delete [] s;

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

void companion_legendre_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_legendre_test() tests companion_legendre().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *A;
  double *C;
  int n = 5;
  double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 };
  double *q;
  complex <double> *r;
  complex <double> *s;

  cout << "\n";
  cout << "companion_legendre_test():\n";
  cout << "  companion_legendre() computes the companion matrix\n";
  cout << "  of a polynomial p(x) in the Legendre basis.\n";

  polynomial_legendre_print ( n, p, "  p(x)" );
//
//  Compute monomial form so we can get roots directly.
//
  A = legendre_to_monomial_matrix ( n + 1 );
  q = r8mat_mv_new ( n + 1, n + 1, A, p );
  polynomial_monomial_print ( n, q, "  Monomial q(x) =" );

  r = r8poly_roots ( n, q );
  c8vec_print ( n, r, "  Roots of q(x):" );

  C = companion_legendre ( n, p );
  r8mat_print ( n, n, C, "  Legendre companion matrix C(p):" );

  s = eigs ( n, C );
  c8vec_print ( n, s, "  Eigenvalues of C(p):" );

  delete [] A;
  delete [] C;
  delete [] q;
  delete [] r;
  delete [] s;

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

void companion_monomial_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    companion_monomial_test() tests companion_monomial().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
{
  double *C;
  int n = 5;
  double p[6] = { 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 };
  complex <double> *r;
  complex <double> *s;

  cout << "\n";
  cout << "companion_monomial_test():\n";
  cout << "  companion_monomial() computes the companion matrix\n";
  cout << "  of a polynomial p(x) in the monomial basis.\n";

  polynomial_monomial_print ( n, p, "  p(x) =" );
//
//  Get roots directly.
//
  r = r8poly_roots ( n, p );
  c8vec_print ( n, r, "  Roots of q(x):" );

  C = companion_monomial ( n, p );
  r8mat_print ( n, n, C, "  Monomial companion matrix C(p):" );

  s = eigs ( n, C );
  c8vec_print ( n, s, "  Eigenvalues of C(p):" );

  delete [] C;
  delete [] r;
  delete [] s;

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

double *chebyshev_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    chebyshev_to_monomial_matrix() returns the chebyshev_to_monomial matrix.
//
//  Discussion
//
//    This is the Chebyshev T matrix, associated with the Chebyshev
//    "T" polynomials, or Chebyshev polynomials of the first kind.
//
//  Example:
//
//    N = 11
//
//    1  .  -1    .    1    .   -1    .     1    .    -1
//    .  1   .   -3    .    5    .   -7     .    9     .
//    .  .   2    .   -8    .   18    .   -32    .    50
//    .  .   .    4    .  -20    .   56     . -120     .
//    .  .   .    .    8    .  -48    .   160    .  -400
//    .  .   .    .    .   16    . -112     .  432     .
//    .  .   .    .    .    .   32    .  -256    .  1120
//    .  .   .    .    .    .    .   64     . -576     .
//    .  .   .    .    .    .    .    .   128    . -1280
//    .  .   .    .    .    .    .    .     .  256     .
//    .  .   .    .    .    .    .    .     .    .   512
//
//  Properties:
//
//    A is generally not symmetric: A' /= A.
//
//    A is integral, therefore det ( A ) is integral, and 
//    det ( A ) * inverse ( A ) is integral.
//
//    A is reducible.
//
//    A is lower triangular.
//
//    Each row of A sums to 1.
//
//    det ( A ) = 2^( (N-1) * (N-2) / 2 )
//
//    A is not normal: A' * A /= A * A'.
//
//    For I = 1:
//      LAMBDA(1) = 1
//    For 1 < I
//      LAMBDA(I) = 2^(I-2)
//
//    The family of matrices is nested as a function of N.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double CHEBY_T[N*N]: the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

  for ( i = 0; i < n; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      A[i+j*n] = 0.0;
    }
  }

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    i = 0;
    A[i+j*n] = - A[i+(j-2)*n];
    for ( i = 1; i < n; i++ )
    {
      A[i+j*n] = 2.0 * A[i-1+(j-1)*n] - A[i+(j-2)*n];
    }
  }
  return A;
}
//****************************************************************************80

double *gegenbauer_to_monomial_matrix ( int n, double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    05 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N: the order of A.
//
//    double alpha: the parameter.
//
//  Output:
//
//    double A(N,N): the matrix.
//
{
  double *A;
  double c1;
  double c2;
  int i;
  int j;
  int nn;

  A = new double[n*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      A[i+j*n] = 0.0;
    }
  }

  if ( n <= 0 )
  {
    return A;
  }

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 2.0 * alpha;
//
//  Perform convex sum.
//  Translating "(n+1) C(n+1) = 2 (n+alpha) x C(n) - ( n + 2 alpha - 1 ) C(n-1)"
//  drove me nuts, between indexing at 1 rather than 0, and dealing with
//  the interpretation of "n+1", because I now face the rare "off by 2" error!
//
  for ( j = 2; j < n; j++ )
  {
    nn = j - 1;
    c1 = ( 2.0 * nn + 2.0 * alpha       ) / ( nn + 1 );
    c2 = (     - nn - 2.0 * alpha + 1.0 ) / ( nn + 1 );
    for ( i = 1; i <= j; i++ )
    {
      A[i+j*n] = c1 * A[i-1+(j-1)*n];
    }
    for ( i = 0; i <= j - 2; i++ )
    {
      A[i+j*n] = A[i+j*n] + c2 * A[i+(j-2)*n];
    }
  }

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

double *hermite_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    hermite_to_monomial_matrix() returns the HERMITE matrix.
//
//  Example:
//
//    N = 11 (Each column must be divided by the factor below it).
//
//    1     .     2     .    12     .   120     .  1680     . 30240
//    .     1     .     6     .    60     .   840     . 15120     .
//    .     .     1     .    12     .   180     .  3360     . 75600
//    .     .     .     1     .    20     .   420     . 10080     .
//    .     .     .     .     1     .    30     .   840     . 25200
//    .     .     .     .     .     1     .    42     .  1512     .
//    .     .     .     .     .     .     1     .    56     .  2520
//    .     .     .     .     .     .     .     1     .    72     .
//    .     .     .     .     .     .     .     .     1     .    90
//    .     .     .     .     .     .     .     .     .     1     .
//    .     .     .     .     .     .     .     .     .     .     1
//
//   /1    /2    /4    /8   /16   /32   /64  /128  /256  /512 /1024
//
//  Properties:
//
//    A is generally not symmetric: A' /= A.
//
//    A is lower triangular.
//
//    det ( A ) = 2^((N*(N-1))/2)
//
//    LAMBDA(I) = 2^(I-1).
//
//    A is integral, therefore det ( A ) is integral, and 
//    det ( A ) * inverse ( A ) is integral.
//
//    A is reducible.
//
//    Successive diagonals are zero, positive, zero, negative.
//
//    A is generally not normal: A' * A /= A * A'.
//
//    The family of matrices is nested as a function of N.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double HERMITE[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      A[i+j*n] = 0.0;
    }
  }

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 2.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == 0 )
      {
        A[i+j*n] = - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n];
      }
      else
      {
        A[i+j*n] = - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n] 
          + 2.0 * A[i-1+(j-1)*n];
      }
    }
  }

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

double *laguerre_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    laguerre_to_monomial_matrix() returns the LAGUERRE matrix.
//
//  Example:
//
//    N = 8 (each column must be divided by the factor below it.)
//
//      1      1      2      6     24    120    720   5040
//      .     -1     -4    -18    -96   -600  -4320 -35280
//      .      .      1      9     72    600   5400  52920
//      .      .      .      1    -16   -200  -2400 -29400
//      .      .      .      .      1     25    450   7350
//      .      .      .      .      .     -1    -36   -882
//      .      .      .      .      .      .      1     49
//      .      .      .      .      .      .      .     -1
//
//     /1     /1     /2     /6    /24   /120   /720  /5040
//
//  Properties:
//
//    A is generally not symmetric: A' /= A.
//
//    A is lower triangular.
//
//    The columns of A alternate in sign.
//
//    A(I,1) = 1
//    A(I,I) = (-1)^(I-1) / (I-1).
//
//    LAMBDA(I) = (-1)^(I-1) / (I-1).
//
//    det ( A ) = product ( 1 <= I <= N ) (-1)^(I-1) / (I-1)
//
//    The I-th row contains the coefficients of the Laguerre
//    polynomial of degree I-1.
//
//    The family of matrices is nested as a function of N.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    US Department of Commerce, 1964.
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double LAGUERRE[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      A[i+j*n] = 0.0;
    }
  }

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[0+1*n] =   1.0;
  A[1+1*n] = - 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    A[0+j*n] = ( ( double ) ( 2 * j - 1 ) * A[0+(j-1)*n]
               + ( double ) (   - j + 1 ) * A[0+(j-2)*n] )
               / ( double ) (     j     );
    for ( i = 1; i < n; i++ )
    {
      A[i+j*n] = ( ( double ) ( 2 * j - 1 ) * A[i+(j-1)*n] 
                 - ( double ) (         1 ) * A[i-1+(j-1)*n]
                 + ( double ) (   - j + 1 ) * A[i+(j-2)*n] ) 
                 / ( double ) (     j     );
    }
  }
  return A;
}
//****************************************************************************80

double *legendre_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    legendre_to_monomial_matrix() returns the LEGENDRE matrix.
//
//  Example:
//
//    N = 11 (each column must be divided by factor at bottom)
//
//     1    .    -1     .      3     .     -5      .      35     .   -63
//     .    1     .    -3      .    15      .    -25       .   315     .
//     .    .     3     .    -30     .    105      .   -1260     .  3465
//     .    .     .     5      .   -70      .    315       . -4620     .
//     .    .     .     .     35     .   -315      .    6930     .-30030
//     .    .     .     .      .    63      .   -693       . 18018     .
//     .    .     .     .      .     .    231      .  -12012     . 90090
//     .    .     .     .      .     .      .    429       .-25740     .
//     .    .     .     .      .     .      .      .    6435     -109395
//     .    .     .     .      .     .      .      .       . 12155     .
//     .    .     .     .      .     .      .      .       .     . 46189
//
//    /1   /1    /2    /2     /8    /8    /16    /16    /128  /128  /256
//
//  Properties:
//
//    A is generally not symmetric: A' /= A.
//
//    A is lower triangular.
//
//    The elements of each row sum to 1.
//
//    Because it has a constant row sum of 1,
//    A has an eigenvalue of 1, and
//    a (right) eigenvector of ( 1, 1, 1, ..., 1 ).
//
//    A is reducible.
//
//    The diagonals form a pattern of zero, positive, zero, negative.
//
//    The family of matrices is nested as a function of N.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double LEGENDRE_MATRIX[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      A[i+j*n] = 0.0;
    }
  }

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 3; j <= n; j++ )
  {
    for ( i = 1; i <= n; i++ )
    {
      if ( i == 1 )
      {
        A[i-1+(j-1)*n] = - ( j - 2 ) * A[i-1+(j-3)*n] 
                         / ( j - 1 );
      }
      else
      {
        A[i-1+(j-1)*n] = ( ( 2 * j - 3 ) * A[i-2+(j-2)*n] 
                         + (   - j + 2 ) * A[i-1+(j-3)*n] ) 
                         / (     j - 1 );
      }
    }
  }

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

void polynomial_chebyshev_print ( int d, double p[], string label )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_chebyshev_print() prints a Chebyshev polynomial.
//
//  Discussion:
//
//    The power sum form is:
//
//      p(x) = a(0)*T0(x) + a(1) * T1(x) + ... + a(n-1) * T(n-1)(x) + a(n) * Tn(x)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer d: the degree of the polynomial.
//
//    real p(0:d), the polynomial coefficients.
//    p(0) is the constant term and
//    p(N) is the coefficient of X^N.
//
//    string label: a title.
//
{
  int i;
  double mag;
  char plus_minus;

  cout << "\n";
  cout << "  " << label << "\n";

  if ( d < 0 )
  {
    cout << "  0\n";
    return;
  }

  for ( i = d; 0 <= i; i-- )
  {
    if ( p[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( p[i] );

    if ( mag != 0.0 )
    {
      cout << "         " << plus_minus << mag << " * T" << i << "(x)\n";
    }
  }
  return;
}
//****************************************************************************80

void polynomial_gegenbauer_print ( int d, double p[], string label )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_gegenbauer_print() prints a Gegenbauer polynomial.
//
//  Discussion:
//
//    The power sum form is:
//
//      p(x) = a(0)*C0(x) + a(1) * C1(x) + ... + a(n-1) * C(n-1)(x) + a(n) * Cn(x)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer d: the degree of the polynomial.
//
//    real p(0:d), the polynomial coefficients.
//    p(0) is the constant term and
//    p(N) is the coefficient of X^N.
//
//    string label: a title.
//
{
  int i;
  double mag;
  char plus_minus;

  cout << "\n";
  cout << "  " << label << "\n";

  if ( d < 0 )
  {
    cout << "  0\n";
    return;
  }

  for ( i = d; 0 <= i; i-- )
  {
    if ( p[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( p[i] );

    if ( mag != 0.0 )
    {
      cout << "         " << plus_minus << mag << " * C" << i << "(x)\n";
    }
  }
  return;
}
//****************************************************************************80

void polynomial_hermite_print ( int d, double p[], string label )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_hermite_print() prints a Hermite polynomial.
//
//  Discussion:
//
//    The power sum form is:
//
//      p(x) = a(0)*H0(x) + a(1) * H1(x) + ... + a(n-1) * H(n-1)(x) + a(n) * Hn(x)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer d: the degree of the polynomial.
//
//    real p(0:d), the polynomial coefficients.
//    p(0) is the constant term and
//    p(N) is the coefficient of X^N.
//
//    string label: a title.
//
{
  int i;
  double mag;
  char plus_minus;

  cout << "\n";
  cout << "  " << label << "\n";

  if ( d < 0 )
  {
    cout << "  0\n";
    return;
  }

  for ( i = d; 0 <= i; i-- )
  {
    if ( p[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( p[i] );

    if ( mag != 0.0 )
    {
      cout << "         " << plus_minus << mag << " * H" << i << "(x)\n";
    }
  }
  return;
}
//****************************************************************************80

void polynomial_laguerre_print ( int d, double p[], string label )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_laguerre_print() prints a Laguerre polynomial.
//
//  Discussion:
//
//    The power sum form is:
//
//      p(x) = a(0)*L0(x) + a(1) * L1(x) + ... + a(n-1) * L(n-1)(x) + a(n) * Ln(x)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer d: the degree of the polynomial.
//
//    real p(0:d), the polynomial coefficients.
//    p(0) is the constant term and
//    p(N) is the coefficient of X^N.
//
//    string label: a title.
//
{
  int i;
  double mag;
  char plus_minus;

  cout << "\n";
  cout << "  " << label << "\n";

  if ( d < 0 )
  {
    cout << "  0\n";
    return;
  }

  for ( i = d; 0 <= i; i-- )
  {
    if ( p[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( p[i] );

    if ( mag != 0.0 )
    {
      cout << "         " << plus_minus << mag << " * L" << i << "(x)\n";
    }
  }
  return;
}
//****************************************************************************80

void polynomial_legendre_print ( int d, double p[], string label )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_legendre_print() prints a Legendre polynomial.
//
//  Discussion:
//
//    The power sum form is:
//
//      p(x) = a(0)*P0(x) + a(1) * P1(x) + ... + a(n-1) * P(n-1)(x) + a(n) * Pn(x)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer d: the degree of the polynomial.
//
//    real p(0:d), the polynomial coefficients.
//    p(0) is the constant term and
//    p(N) is the coefficient of X^N.
//
//    string label: a title.
//
{
  int i;
  double mag;
  char plus_minus;

  cout << "\n";
  cout << "  " << label << "\n";

  if ( d < 0 )
  {
    cout << "  0\n";
    return;
  }

  for ( i = d; 0 <= i; i-- )
  {
    if ( p[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( p[i] );

    if ( mag != 0.0 )
    {
      cout << "         " << plus_minus << mag << " * P" << i << "(x)\n";
    }
  }
  return;
}
//****************************************************************************80

void polynomial_monomial_print ( int d, double p[], string label )

//****************************************************************************80
//
//  Purpose:
//
//    polynomial_monomial_print() prints a polynomial.
//
//  Discussion:
//
//    The power sum form is:
//
//      p(x) = a(0) + a(1) * x + ... + a(n-1) * x^(n-1) + a(n) * x^(n)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer d: the degree of the polynomial.
//
//    real p(0:d), the polynomial coefficients.
//    p(0) is the constant term and
//    p(N) is the coefficient of X^N.
//
//    string label: a title.
//
{
  int i;
  double mag;
  char plus_minus;

  cout << "\n";
  cout << "  " << label << "\n";

  if ( d < 0 )
  {
    cout << "  0\n";
    return;
  }

  for ( i = d; 0 <= i; i-- )
  {
    if ( p[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( p[i] );

    if ( mag != 0.0 )
    {
      if ( 2 <= i )
      {
        cout << "         " << plus_minus << mag << " * x^" << i << "\n";
      }
      else if ( i == 1 )
      {
        cout << "         " << plus_minus << mag << " * x\n";
      }
      else if ( i == 0 )
      {
        cout << "         " << plus_minus << mag << "\n";
      }
    }
  }
  return;
}
//****************************************************************************80

void c8vec_print ( int n, complex <double> a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    c8vec_print() prints a C8VEC.
//
//  Discussion:
//
//    A C8VEC is a vector of complex <double> values.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of components of the vector.
//
//    complex <double> A[N], the vector to be printed.
//
//    string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";

  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8) << i
         << ": " << real ( a[i] )
         << "  " << imag ( a[i] ) << "\n";
  }

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

double *r8mat_mv_new ( int m, int n, double a[], double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_mv_new() multiplies a matrix times a vector.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    For this routine, the result is returned as the function value.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 April 2007
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the number of rows and columns of the matrix.
//
//    double A[M,N], the M by N matrix.
//
//    double X[N], the vector to be multiplied by A.
//
//  Output:
//
//    double R8MAT_MV_NEW[M], the product A*X.
//
{
  int i;
  int j;
  double *y;

  y = new double[m];

  for ( i = 0; i < m; i++ )
  {
    y[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      y[i] = y[i] + a[i+j*m] * x[j];
    }
  }

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

void r8mat_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_print() prints an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    Entry A(I,J) is stored as A[I+J*M]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the number of rows in A.
//
//    int N, the number of columns in A.
//
//    double A[M*N], the M by N matrix.
//
//    string TITLE, a title.
//
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

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

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_print_some() prints some of an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 June 2013
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, the number of rows of the matrix.
//    M must be positive.
//
//    int N, the number of columns of the matrix.
//    N must be positive.
//
//    double A[M*N], the matrix.
//
//    int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    string TITLE, a title.
//
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }
//
//  Print the columns of the matrix, in strips of 5.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }
    cout << "\n";
//
//  For each column J in the current range...
//
//  Write the header.
//
    cout << "  Col:    ";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(7) << j - 1 << "       ";
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "\n";
//
//  Determine the range of the rows in this strip.
//
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( ihi < m )
    {
      i2hi = ihi;
    }
    else
    {
      i2hi = m;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
//
//  Print out (up to) 5 entries in row I, that lie in the current strip.
//
      cout << setw(5) << i - 1 << ": ";
      for ( j = j2lo; j <= j2hi; j++ )
      {
        cout << setw(12) << a[i-1+(j-1)*m] << "  ";
      }
      cout << "\n";
    }
  }

  return;
# undef INCX
}
//****************************************************************************80

complex <double> *r8poly_roots ( int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8poly_roots() returns the roots of a polynomial.
//
//  Discussion:
//
//    The monomial form of a polynomial is:
//
//      p(x) = c(0)*x^0 + c(1)*x + ... + c(n)*x^(n)
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    11 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, constant first.
//
//  Output:
//
//    complex <double> r(n): the polynomial roots.
//
{
  int i;
  complex <double> *r;
  double *z;
//
//  Use GSL to compute roots.
//
  z = new double[2*n];
      
  gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc ( n + 1 );
  
  gsl_poly_complex_solve ( p, n + 1, w, z );

  gsl_poly_complex_workspace_free ( w );
//
//  Transfer pairs of real, imaginary parts to complex <double> vector.
//
  r = new complex <double>[n];
  for ( i = 0; i < n; i++ )
  {
    r[i] = complex <double> ( z[2*i], z[2*i+1] );
  }

  delete [] z;

  return r;
}
//****************************************************************************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
}