subroutine companion_chebyshev ( n, p, A )

c*********************************************************************72
c
cc companion_chebyshev() returns the Chebyshev basis companion matrix for a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    John Boyd,
c    Solving Transcendental Equations,
c    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
c    Perturbation Series, and Oracles,
c    SIAM, 2014,
c    ISBN: 978-1-611973-51-8,
c    LC: QA:353.T7B69
c
c  Input:
c
c    real p(n+1): the polynomial coefficients, in order of increasing degree.
c
c  Output:
c
c    real A(n,n): the companion matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      integer i
      integer j
      double precision p(n+1)

      A(1:n,1:n) = 0.0D+00

      A(1,2) = 1.0D+00

      do i = 2, n - 1
        A(i,i-1) = 0.5D+00
        A(i,i+1) = 0.5D+00
      end do

      do j = 1, n
        A(n,j) = - 0.5D+00 * p(j) / p(n+1)
      end do

      A(n,n-1) = A(n,n-1) + 0.5D+00

      return
      end
      subroutine companion_gegenbauer ( n, p, alpha, A )

c*********************************************************************72
c
cc companion_gegenbauer() returns the Gegenbauer basis companion matrix for a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    John Boyd,
c    Solving Transcendental Equations,
c    The gegenbauer Polynomial Proxy and other Numerical Rootfinders,
c    Perturbation Series, and Oracles,
c    SIAM, 2014,
c    ISBN: 978-1-611973-51-8,
c    LC: QA:353.T7B69
c
c  Input:
c
c    real p(n+1): the polynomial coefficients, in order of increasing degree.
c
c    real alpha: the Gegenbauer parameter.
c
c  Output:
c
c    real A(n,n): the companion matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      double precision alpha
      integer i
      integer j
      double precision p(n+1)

      A(1:n,1:n) = 0.0D+00

      A(1,2) = 0.5D+00 / alpha

      do i = 2, n - 1
        A(i,i-1) = 0.5 * ( i - 1.0 + 2.0 * alpha - 1.0 ) 
     &    / ( i - 1.0 + alpha )
        A(i,i+1) = 0.5 * ( i                     ) 
     &    / ( i - 1.0 + alpha )
      end do

      do j = 1, n
        A(n,j) = - 0.5D+00 * p(j) * n / ( n - 1.0 + alpha ) / p(n+1)
      end do

      A(n,n-1) = A(n,n-1) + 0.5D+00 
     &  * ( n - 2.0 + 2.0 * alpha ) / ( n - 1.0 + alpha )

      return
      end
      subroutine companion_hermite ( n, p, A )

c*********************************************************************72
c
cc companion_hermite() returns the Hermite basis companion matrix for a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    John Boyd,
c    Solving Transcendental Equations,
c    The legendre Polynomial Proxy and other Numerical Rootfinders,
c    Perturbation Series, and Oracles,
c    SIAM, 2014,
c    ISBN: 978-1-611973-51-8,
c    LC: QA:353.T7B69
c
c  Input:
c
c    real p(n+1): the polynomial coefficients, in order of increasing degree.
c
c  Output:
c
c    real A(n,n): the companion matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      integer i
      integer j
      double precision p(n+1)

      A(1:n,1:n) = 0.0D+00

      do i = 1, n - 1
        A(i,i+1) = 0.5
      end do

      do i = 2, n
        A(i,i-1) = i - 1
      end do

      do j = 1, n
        A(n,j) = A(n,j) - p(j) / 2.0 * p(n+1)
      end do

      return
      end
      subroutine companion_laguerre ( n, p, A )

c*********************************************************************72
c
cc companion_laguerre() returns the Laguerre basis companion matrix for a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    John Boyd,
c    Solving Transcendental Equations,
c    The gegenbauer Polynomial Proxy and other Numerical Rootfinders,
c    Perturbation Series, and Oracles,
c    SIAM, 2014,
c    ISBN: 978-1-611973-51-8,
c    LC: QA:353.T7B69
c
c  Input:
c
c    real p(n+1): the polynomial coefficients, in order of increasing degree.
c
c  Output:
c
c    real A(n,n): the companion matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      integer i
      integer j
      double precision p(n+1)

      A(1:n,1:n) = 0.0D+00

      do i = 1, n
        A(i,i) = 2 * i - 1
      end do

      do i = 2, n
        A(i,i-1) = - i + 1
      end do

      do i = 1, n - 1
        A(i,i+1) = - i
      end do

      do j = 1, n
        A(n,j) = A(n,j) + n * p(j) / p(n+1)
      end do

      return
      end
      subroutine companion_legendre ( n, p, A )

c*********************************************************************72
c
cc companion_legendre() returns the Legendre basis companion matrix for a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    John Boyd,
c    Solving Transcendental Equations,
c    The legendre Polynomial Proxy and other Numerical Rootfinders,
c    Perturbation Series, and Oracles,
c    SIAM, 2014,
c    ISBN: 978-1-611973-51-8,
c    LC: QA:353.T7B69
c
c  Input:
c
c    real p(n+1): the polynomial coefficients, in order of increasing degree.
c
c  Output:
c
c    real A(n,n): the companion matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      integer i
      integer j
      double precision p(n+1)

      A(1:n,1:n) = 0.0D+00

      A(1,2) = 1.0

      do i = 2, n - 1
        A(i,i-1) = dble ( i - 1 ) / dble ( 2 * i - 1 )
        A(i,i+1) = dble ( i ) / dble ( 2 * i - 1 )
      end do

      do j = 1, n
        A(n,j) = - p(j) / p(n+1) * dble ( n ) / dble ( 2 * n - 1 )
      end do

      A(n,n-1) = A(n,n-1) + dble ( n - 1 ) / dble ( 2 * n - 1 )

      return
      end
      subroutine companion_monomial ( n, p, A )

c*********************************************************************72
c
cc companion_monomial() returns the monomial basis companion matrix for a polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    John Boyd,
c    Solving Transcendental Equations,
c    The Chebyshev Polynomial Proxy and other Numerical Rootfinders,
c    Perturbation Series, and Oracles,
c    SIAM, 2014,
c    ISBN: 978-1-611973-51-8,
c    LC: QA:353.T7B69
c
c  Input:
c
c    real p(n+1): the polynomial coefficients, in order of increasing degree.
c
c  Output:
c
c    real A(n,n): the companion matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      integer i
      integer j
      double precision p(n+1)

      A(1:n,1:n) = 0.0D+00

      do i = 1, n - 1
        A(i,i+1) = 1.0D+00
      end do

      do j = 1, n
        A(n,j) = - p(j) / p(n+1)
      end do

      return
      end