# include <math.h>
# include <stdlib.h>

# include "companion_matrix.h"

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

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

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

    companion_chebyshev() returns the Chebyshev basis companion matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 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 = ( double * ) malloc ( n * n * sizeof ( double ) );

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

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

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

    companion_gegenbauer() returns the Gegenbauer basis companion matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 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 = ( double * ) malloc ( n * n * sizeof ( double ) );

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

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

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

    companion_hermite() returns the Hermite basis companion matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    09 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 = ( double * ) malloc ( n * n * sizeof ( double ) );

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

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

/******************************************************************************/
/*
  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 = ( double * ) malloc ( n * n * sizeof ( double ) );

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

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

/******************************************************************************/
/*
  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 = ( double * ) malloc ( n * n * sizeof ( double ) );

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

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

/******************************************************************************/
/*
  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 = ( double * ) malloc ( n * n * sizeof ( double ) );

  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;
}