# include <cmath>
# include <cstdlib>

using namespace std;

# include "companion_matrix.hpp"

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

double *companion_chebyshev ( int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    companion_chebyshev() returns the Chebyshev basis companion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Boyd,
//    Solving Transcendental Equations,
//    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
//    Perturbation Series, and Oracles,
//    SIAM, 2014,
//    ISBN: 978-1-611973-51-8,
//    LC: QA:353.T7B69
//
//  Input:
//
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, in order of increasing degree.
//
//  Output:
//
//    real C(n,n): the companion matrix.
//
{
  double *C;
  int i;
  int j;

  C = new double[n*n];

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

  C[0+1*n] = 1.0;

  for ( i = 1; i < n - 1; i++ )
  {
    C[i+(i-1)*n] = 0.5;
    C[i+(i+1)*n] = 0.5;
  }

  for ( j = 0; j < n; j++ )
  {
    C[n-1+j*n] = - 0.5 * p[j] / p[n];
  }

  C[n-1+(n-2)*n] = C[n-1+(n-2)*n] + 0.5;

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

double *companion_gegenbauer ( int n, double p[], double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    companion_gegenbauer() returns the Gegenbauer basis companion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Boyd,
//    Solving Transcendental Equations,
//    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
//    Perturbation Series, and Oracles,
//    SIAM, 2014,
//    ISBN: 978-1-611973-51-8,
//    LC: QA:353.T7B69
//
//  Input:
//
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, in order of increasing degree.
//
//    real alpha: the Gegenbauer parameter.
//
//  Output:
//
//    real C(n,n): the companion matrix.
//
{
  double *C;
  int i;
  int j;

  C = new double[n*n];

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

  C[0+1*n] = 0.5 / alpha;

  for ( i = 1; i < n - 1; i++ )
  {
    C[i+(i-1)*n] = 0.5 * ( i + 2 * alpha - 1 ) / ( i + alpha );
    C[i+(i+1)*n] = 0.5 * ( i             + 1 ) / ( i + alpha );
  }

  for ( j = 0; j < n; j++ )
  {
    C[n-1+j*n] = - 0.5 * p[j] * n / ( n - 1 + alpha ) / p[n];
  }

  C[n-1+(n-2)*n] = C[n-1+(n-2)*n] 
    + 0.5 * ( n - 2 + 2 * alpha ) / ( n - 1 + alpha );

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

double *companion_hermite ( int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    companion_hermite() returns the Hermite basis companion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Boyd,
//    Solving Transcendental Equations,
//    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
//    Perturbation Series, and Oracles,
//    SIAM, 2014,
//    ISBN: 978-1-611973-51-8,
//    LC: QA:353.T7B69
//
//  Input:
//
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, in order of increasing degree.
//
//  Output:
//
//    real C(n,n): the companion matrix.
//
{
  double *C;
  int i;
  int j;

  C = new double[n*n];

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

  for ( i = 0; i < n - 1; i++ )
  {
    C[i+(i+1)*n] = 0.5;
  }

  for ( i = 1; i < n; i++ )
  {
    C[i+(i-1)*n] = i;
  }

  for ( j = 0; j < n; j++ )
  {
    C[n-1+j*n] = C[n-1+j*n] - p[j] / 2.0 * p[n];
  }

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

double *companion_laguerre ( int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    companion_laguerre() returns the Laguerre basis companion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Boyd,
//    Solving Transcendental Equations,
//    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
//    Perturbation Series, and Oracles,
//    SIAM, 2014,
//    ISBN: 978-1-611973-51-8,
//    LC: QA:353.T7B69
//
//  Input:
// 
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, in order of increasing degree.
//
//  Output:
//
//    real C(n,n): the companion matrix.
//
{
  double *C;
  int i;
  int j;

  C = new double[n*n];

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

  for ( i = 0; i < n; i++ )
  {
    C[i+i*n] = 2 * i + 1;
  }

  for ( i = 1; i < n; i++ )
  {
    C[i+(i-1)*n] = - i;
  }

  for ( i = 0; i < n - 1; i++ )
  {
    C[i+(i+1)*n] = - i - 1;
  }

  for ( j = 0; j < n; j++ )
  {
    C[n-1+j*n] = C[n-1+j*n] + n * p[j] / p[n];
  }

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

double *companion_legendre ( int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    companion_legendre() returns the Legendre basis companion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Boyd,
//    Solving Transcendental Equations,
//    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
//    Perturbation Series, and Oracles,
//    SIAM, 2014,
//    ISBN: 978-1-611973-51-8,
//    LC: QA:353.T7B69
//
//  Input:
//
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, in order of increasing degree.
//
//  Output:
//
//    real C(n,n): the companion matrix.
//
{
  double *C;
  int i;
  int j;

  C = new double[n*n];

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

  C[0+1*n] = 1.0;

  for ( i = 1; i < n - 1; i++ )
  {
    C[i+(i-1)*n] = ( double ) ( i     ) / ( double ) ( 2 * i + 1 );
    C[i+(i+1)*n] = ( double ) ( i + 1 ) / ( double ) ( 2 * i + 1 );
  }

  for ( j = 0; j < n; j++ )
  {
    C[n-1+j*n] = - p[j] / p[n] * ( double ) ( n ) / ( double ) ( 2 * n - 1 );
  }

  C[n-1+(n-2)*n] = C[n-1+(n-2)*n] + ( double ) ( n - 1 ) 
    / ( double ) ( 2 * n - 1 );

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

double *companion_monomial ( int n, double p[] )

//****************************************************************************80
//
//  Purpose:
//
//    companion_monomial() returns the monomial basis companion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 June 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Boyd,
//    Solving Transcendental Equations,
//    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
//    Perturbation Series, and Oracles,
//    SIAM, 2014,
//    ISBN: 978-1-611973-51-8,
//    LC: QA:353.T7B69
//
//  Input:
//
//    integer n: the polynomial degree.
//
//    real p(n+1): the polynomial coefficients, in order of increasing degree.
//
//  Output:
//
//    real C(n,n): the companion matrix.
//
{
  double *C;
  int i;
  int j;

  C = new double[n*n];

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

  for ( i = 0; i < n - 1; i++ )
  {
    C[i+(i+1)*n] = 1.0;
  }

  for ( j = 0; j < n; j++ )
  {
    C[n-1+j*n] = - p[j] / p[n];
  }

  return C;
}