program main

c*********************************************************************72
c
cc companion_matrix_test() tests companion_matrix().
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
      call timestamp ( )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_matrix_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) 
     &  '  companion_matrix() computes the companion matrix'
      write ( *, '(a)' ) '  of a polynomial, in various bases.'

      call companion_chebyshev_test ( )
      call companion_gegenbauer_test ( )
      call companion_hermite_test ( )
      call companion_laguerre_test ( )
      call companion_legendre_test ( )
      call companion_monomial_test ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_matrix_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ''
      call timestamp ( )

      stop ( 0 )
      end
      subroutine companion_chebyshev_test ( )

c*********************************************************************72
c
cc companion_chebyshev_test() tests companion_chebyshev().
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
      implicit none

      integer, parameter :: n = 5

      double precision A(n+1,n+1)
      double precision C(n,n)
      double precision p(n+1)
      double precision q(n+1)
      double complex r(n)
      double complex r1(n)

      data p /
     &  6.0,
     &  5.0,
     &  4.0,
     &  3.0,
     &  2.0,
     &  1.0 /

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_chebyshev_test():'
      write ( *, '(a)' ) 
     &  '  companion_chebyshev() computes the companion matrix'
      write ( *, '(a)' ) 
     &  '  of a polynomial p(x) in the Chebyshev basis.'

      call polynomial_chebyshev_print ( n, p, '  p(x)' )
c
c  Compute monomial form so we can get roots directly.
c
      call chebyshev_to_monomial_matrix ( n + 1, A )
      q = matmul ( A, p )
      call polynomial_monomial_print ( n, q, '  Monomial q(x)' )

      call r8poly_roots ( n, q, r1 )
      call c8vec_print ( n, r1, '  roots of p(x)' )

      call companion_chebyshev ( n, p, C )
      call r8mat_print ( n, n, C, '  Chebyshev companion matrix C(p):' )

      call eigs ( n, C, r )
      call c8vec_print ( n, r, '  Eigenvalues of C(p)' )

      return
      end
      subroutine companion_gegenbauer_test ( )

c*********************************************************************72
c
cc companion_gegenbauer_test() tests companion_gegenbauer().
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
      implicit none

      integer, parameter :: n = 5

      double precision A(n+1,n+1)
      double precision alpha
      double precision C(n,n)
      integer i
      double precision p(n+1)
      double precision q(n+1)
      double complex r(n)
      double complex r1(n)

      data p /
     &  6.0,
     &  5.0,
     &  4.0,
     &  3.0,
     &  2.0,
     &  1.0 /

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_gegenbauer_test():'
      write ( *, '(a)' ) 
     &  '  companion_gegenbauer() computes the companion matrix'
      write ( *, '(a)' ) 
     &  '  of a polynomial p(x) in the Gegenbauer basis.'

      call polynomial_gegenbauer_print ( n, p, '  p(x)' )
c
c  Try several parameter values.
c
      do i = 1, 2

        if ( i == 1 ) then
          alpha = 0.5
        else if ( i == 2 ) then
          alpha = 1.0
        end if

        write ( *, '(a)' ) ''
        write ( *, '(a,g14.6)' ) '  Parameter alpha will be', alpha
c
c  Compute monomial form so we can get roots directly.
c
        call gegenbauer_to_monomial_matrix ( n + 1, alpha, A )
        q = matmul ( A, p )
        call polynomial_monomial_print ( n, q, '  Monomial q(x)' )

        call r8poly_roots ( n, q, r1 )
        call c8vec_print ( n, r1, '  roots of p(x)' )

        call companion_gegenbauer ( n, p, alpha, C )
        call r8mat_print ( n, n, C, 
     &    '  Gegenbauer companion matrix C(p)' )

        call eigs ( n, C, r )
        call c8vec_print ( n, r, '  Eigenvalues of C(p)' )

      end do

      return
      end
      subroutine companion_hermite_test ( )

c*********************************************************************72
c
cc companion_hermite_test() tests companion_hermite().
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
      implicit none

      integer, parameter :: n = 5

      double precision A(n+1,n+1)
      double precision C(n,n)
      double precision p(n+1)
      double precision q(n+1)
      double complex r(n)
      double complex r1(n)

      data p /
     &  6.0,
     &  5.0,
     &  4.0,
     &  3.0,
     &  2.0,
     &  1.0 /

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_hermite_test():'
      write ( *, '(a)' ) 
     &  '  companion_hermite() computes the companion matrix'
      write ( *, '(a)' ) 
     &  '  of a polynomial p(x) in the Hermite basis.'

      call polynomial_hermite_print ( n, p, '  Hermite p(x)' )
c
c  Compute monomial form so we can get roots directly.
c
      call h_to_monomial_matrix ( n + 1, A )
      q = matmul ( A, p )
      call polynomial_monomial_print ( n, q, '  Monomial q(x)' )

      call r8poly_roots ( n, q, r1 )
      call c8vec_print ( n, r1, '  roots of p(x)' )
     
      call companion_hermite ( n, p, C )
      call r8mat_print ( n, n, C, '  Hermite companion matrix C(p):' )

      call eigs ( n, C, r )
      call c8vec_print ( n, r, '  Eigenvalues of C(p)' )

      return
      end
      subroutine companion_laguerre_test ( )

c*********************************************************************72
c
cc companion_laguerre_test() tests companion_laguerre().
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
      implicit none

      integer, parameter :: n = 5

      double precision A(n+1,n+1)
      double precision C(n,n)
      double precision p(n+1)
      double precision q(n+1)
      double complex r(n)
      double complex r1(n)

      data p /
     &  6.0,
     &  5.0,
     &  4.0,
     &  3.0,
     &  2.0,
     &  1.0 /

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_laguerre_test():'
      write ( *, '(a)' ) 
     &  '  companion_laguerre() computes the companion matrix'
      write ( *, '(a)' ) 
     &  '  of a polynomial p(x) in the Laguerre basis.'

      call polynomial_laguerre_print ( n, p, '  Laguerre p(x)' )
c
c  Compute monomial form so we can get roots directly.
c
      call laguerre_to_monomial_matrix ( n + 1, A )
      q = matmul ( A, p )
      call polynomial_monomial_print ( n, q, '  Monomial q(x)' )

      call r8poly_roots ( n, q, r1 )
      call c8vec_print ( n, r1, '  roots of p(x)' )
     
      call companion_laguerre ( n, p, C )
      call r8mat_print ( n, n, C, '  Laguerre companion matrix C(p):' )

      call eigs ( n, C, r )
      call c8vec_print ( n, r, '  Eigenvalues of C(p)' )

      return
      end
      subroutine companion_legendre_test ( )

c*********************************************************************72
c
cc companion_legendre_test() tests companion_legendre().
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
      implicit none

      integer, parameter :: n = 5

      double precision A(n+1,n+1)
      double precision C(n,n)
      double precision p(n+1)
      double precision q(n+1)
      double complex r(n)
      double complex r1(n)

      data p /
     &  6.0,
     &  5.0,
     &  4.0,
     &  3.0,
     &  2.0,
     &  1.0 /

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_legendre_test():'
      write ( *, '(a)' ) 
     &  '  companion_legendre() computes the companion matrix'
      write ( *, '(a)' ) 
     &  '  of a polynomial p(x) in the Legendre basis.'

      call polynomial_legendre_print ( n, p, '  Legendre p(x)' )
c
c  Compute monomial form so we can get roots directly.
c
      call legendre_to_monomial_matrix ( n + 1, A )
      q = matmul ( A, p )
      call polynomial_monomial_print ( n, q, '  Monomial q(x)' )

      call r8poly_roots ( n, q, r1 )
      call c8vec_print ( n, r1, '  roots of p(x)' )
     
      call companion_legendre ( n, p, C )
      call r8mat_print ( n, n, C, '  Legendre companion matrix C(p):' )

      call eigs ( n, C, r )
      call c8vec_print ( n, r, '  Eigenvalues of C(p)' )

      return
      end
      subroutine companion_monomial_test ( )

c*********************************************************************72
c
cc companion_monomial_test() tests companion_monomial().
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
      implicit none

      integer, parameter :: n = 5

      double precision C(n,n)
      double precision p(n+1)
      double complex r(n)
      double complex r1(n)

      data p /
     &  6.0,
     &  5.0,
     &  4.0,
     &  3.0,
     &  2.0,
     &  1.0 /

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'companion_monomial_test():'
      write ( *, '(a)' ) 
     &  '  companion_monomial() computes the companion matrix'
      write ( *, '(a)' ) 
     &  '  of a polynomial p(x) in the monomial basis.'

      call polynomial_monomial_print ( n, p, '  p(x)' )
c
c  Get roots directly.
c
      call r8poly_roots ( n, p, r1 )
      call c8vec_print ( n, r1, '  roots of p(x)' )

      call companion_monomial ( n, p, C )
      call r8mat_print ( n, n, C, '  Monomial companion matrix C(p):' )

      call eigs ( n, C, r )
      call c8vec_print ( n, r, '  Eigenvalues of C(p)' )

      return
      end
      subroutine chebyshev_to_monomial_matrix ( n, A )

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

      integer n

      double precision A(n,n)
      integer i
      integer j

      if ( n .le. 0 ) then
        return
      end if

      do j = 1, n
        do i = 1, n
          A(i,j) = 0.0D+00
        end do
      end do

      A(1,1) = 1.0D+00

      if ( n .eq. 1 ) then
        return
      end if

      A(2,2) = 1.0D+00

      if ( n .eq. 2 ) then
        return
      end if

      do j = 3, n
        do i = 1, n
          if ( i .eq. 1 ) then
            A(i,j) = - A(i,j-2)
          else
            A(i,j) = 2.0D+00 * A(i-1,j-1) - A(i,j-2)
          end if
        end do
      end do

      return
      end
      subroutine gegenbauer_to_monomial_matrix ( n, alpha, A )

c*********************************************************************72
c
cc gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    05 April 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N: the order of A.
c
c    real alpha: the parameter.
c
c  Output:
c
c    double precision A(N,N): the matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      double precision alpha
      double precision c1
      double precision c2
      integer i
      integer j
      integer nn

      if ( n <= 0 ) then
        return
      end if

      do j = 1, n
        do i = 1, n
          A(i,j) = 0.0D+00
        end do
      end do

      A(1,1) = 1.0D+00

      if ( n == 1 ) then
        return
      end if

      A(2,2) = 2.0D+00 * alpha
c
c  Perform convex sum.
c  Translating "(n+1) C(n+1) = 2 (n+alpha) x C(n) - ( n + 2 alpha - 1 ) C(n-1)"
c  drove me nuts, between indexing at 1 rather than 0, and dealing with
c  the interpretation of "n+1", because I now face the rare "off by 2" errorc
c
      do j = 3, n
        nn = j - 2
        c1 = ( 2.0 * nn + 2.0 * alpha       ) / dble ( nn + 1 )
        c2 = (     - nn - 2.0 * alpha + 1.0 ) / dble ( nn + 1 )
        A(2:j,j)   =              c1 * A(1:j-1,j-1)
        A(1:j-2,j) = A(1:j-2,j) + c2 * A(1:j-2,j-2)
      end do

      return
      end
      subroutine h_to_monomial_matrix ( n, a )

c*********************************************************************72
c
cc h_to_monomial_matrix() returns the HERMITE matrix.
c
c  Example:
c
c    N = 11
c
c      1   .   -2      .      12     .   -120     .   1680      .  -30240
c      .   2    .    -12       .   120      . -1680      .  30240       .
c      .   .    4      .     -48     .    720     . -13440      .  302400
c      .   .    .      8       .  -160      .  3360      . -80640       .
c      .   .    .      .      16     .   -480     .  13440      . -403200
c      .   .    .      .       .    32      . -1344      .  48384       .
c      .   .    .      .       .     .     64     .  -3584      .  161280
c      .   .    .      .       .     .      .   128      .  -9216       .
c      .   .    .      .       .     .      .     .    256      .  -23040
c      .   .    .      .       .     .      .     .      .    512       .
c      .   .    .      .       .     .      .     .      .      .    1024
c
c  Properties:
c
c    A is generally not symmetric: A' /= A.
c
c    A is lower triangular.
c
c    det ( A ) = 2^((N*(N-1))/2)
c
c    LAMBDA(I) = 2^(I-1).
c
c    A is integral, therefore det ( A ) is integral, and 
c    det ( A ) * inverse ( A ) is integral.
c
c    A is reducible.
c
c    Successive diagonals are zero, positive, zero, negative.
c
c    A is generally not normal: A' * A /= A * A'.
c
c    The family of matrices is nested as a function of N.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    04 April 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the order of the matrix.
c
c  Output:
c
c    double precision A(N,N), the matrix.
c
      implicit none

      integer n

      double precision a(n,n)
      integer i
      integer j

      if ( n .le. 0 ) then
        return
      end if

      do j = 1, n
        do i = 1, n
          a(i,j) = 0.0D+00
        end do
      end do

      a(1,1) = 1.0D+00

      if ( n .eq. 1 ) then
        return
      end if

      a(2,2) = 2.0D+00

      if ( n .eq. 2 ) then
        return
      end if

      do j = 3, n
        do i = 1, n
          if ( i == 1 ) then
            A(i,j) =                      - 2.0D+00 
     &        * dble ( j - 2 ) * A(i,j-2)
          else
            A(i,j) = 2.0D+00 * a(i-1,j-1) - 2.0D+00 
     &        * dble ( j - 2 ) * A(i,j-2)
          end if
        end do
      end do

      return
      end
      subroutine laguerre_to_monomial_matrix ( n, a )

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

      integer n

      double precision a(n,n)
      integer i
      integer j

      if ( n <= 0 ) then
        return
      end if

      do j = 1, n
        do i = 1, n
          a(i,j) = 0.0D+00
        end do
      end do

      a(1,1) = 1.0D+00

      if ( n .eq. 1 ) then
        return
      end if

      a(1,2) = 1.0D+00
      a(2,2) = -1.0D+00

      if ( n .eq. 2 ) then
        return
      end if

      do j = 3, n
        do i = 1, n
          if ( i == 1 ) then
            A(i,j) = ( dble ( 2 * j - 3 ) * A(i,j-1) 
     &               + dble (   - j + 2 ) * A(i,j-2) ) 
     &               / dble (     j - 1 )
          else
            A(i,j) = ( dble ( 2 * j - 3 ) * A(i,j-1) 
     &               - dble (         1 ) * A(i-1,j-1) 
     &               + dble (   - j + 2 ) * A(i,j-2) ) 
     &               / dble (     j - 1 )
          end if
        end do
      end do

      return
      end
      subroutine legendre_to_monomial_matrix ( n, A )

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

      integer n

      double precision A(n,n)
      integer i
      integer j

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

      if ( n <= 0 ) then
        return
      end if

      A(1,1) = 1.0D+00

      if ( n == 1 ) then
        return
      end if

      A(2,2) = 1.0D+00

      if ( n == 2 ) then
        return
      end if

      do j = 3, n
        do i = 1, n
          if ( i == 1 ) then
            A(i,j) = - ( j - 2 ) * A(i,j-2) 
     &               / ( j - 1 )
          else
            A(i,j) = ( ( 2 * j - 3 ) * A(i-1,j-1) 
     &               + (   - j + 2 ) * A(i,j-2) ) 
     &               / (     j - 1 )
          end if
        end do
      end do

      return
      end
      subroutine polynomial_chebyshev_print ( d, c, label )

c*********************************************************************72
c
cc polynomial_chebyshev_print() prints a polynomial in the Chebyshev basis.
c
c  Discussion:
c
c    The form of a Chebyshev polynomial is:
c
c      p(x) = c(0)*T0(x) + c(1)*T(1)(x) + ... + c(n)*T(n)(x)
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  Input:
c
c    integer D: the polynomial degree.
c
c    real C(1:D+1): the polynomial coefficients, constant first.
c
c    character LABEL: an optional title.
c
      implicit none

      integer d

      double precision c(d+1)
      integer i
      character ( len = * ) label

      write ( *, '(a)' ) ''
      if ( 0 < len_trim ( label ) ) then
        write ( *, '(a)' ) trim ( label ) // ' = '
      end if

      if ( d < 0 ) then
        write ( *, '(a)' ) '  Zero polynomial'
        return
      end if

      if ( all ( c == 0 ) ) then
        write ( *, '(a)' ) '  0'
        return
      end if

      do i = d, 0, -1

        if ( c(i+1) /= 0.0 ) then

          write ( *, '(a,g14.6,a,i2,a)' ) 
     &      '  +', c(i+1), ' * T', i, '(x)'

        end if

      end do

      return
      end
      subroutine polynomial_gegenbauer_print ( d, c, label )

c*********************************************************************72
c
cc polynomial_gegenbauer_print() prints a polynomial in the Gegenbauer basis.
c
c  Discussion:
c
c    The form of a Gegenbauer polynomial is:
c
c      p(x) = c(0)*C0(x) + c(1)*C1(x) + ... + c(n)*C(n)(x)
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  Input:
c
c    integer D: the polynomial degree.
c
c    real C(1:D+1): the polynomial coefficients, constant first.
c
c    character LABEL: an optional title.
c
      implicit none

      integer d

      double precision c(d+1)
      integer i
      character ( len = * ) label

      write ( *, '(a)' ) ''
      if ( 0 < len_trim ( label ) ) then
        write ( *, '(a)' ) trim ( label ) // ' = '
      end if

      if ( d < 0 ) then
        write ( *, '(a)' ) '  Zero polynomial'
        return
      end if

      if ( all ( c == 0 ) ) then
        write ( *, '(a)' ) '  0'
        return
      end if

      do i = d, 0, -1

        if ( c(i+1) /= 0.0 ) then

          write ( *, '(a,g14.6,a,i2,a)' ) 
     &      '  +', c(i+1), ' * C', i, '(x)'

        end if

      end do

      return
      end
      subroutine polynomial_hermite_print ( d, c, label )

c*********************************************************************72
c
cc polynomial_hermite_print() prints a polynomial in the Hermite basis.
c
c  Discussion:
c
c    The form of a Hermite polynomial is:
c
c      p(x) = c(0)*H0(x) + c(1)*H(1)(x) + ... + c(n)*H(n)(x)
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  Input:
c
c    integer D: the polynomial degree.
c
c    real C(1:D+1): the polynomial coefficients, constant first.
c
c    character LABEL: an optional title.
c
      implicit none

      integer d

      double precision c(d+1)
      integer i
      character ( len = * ) label

      write ( *, '(a)' ) ''
      if ( 0 < len_trim ( label ) ) then
        write ( *, '(a)' ) trim ( label ) // ' = '
      end if

      if ( d < 0 ) then
        write ( *, '(a)' ) '  Zero polynomial'
        return
      end if

      if ( all ( c == 0 ) ) then
        write ( *, '(a)' ) '  0'
        return
      end if

      do i = d, 0, -1

        if ( c(i+1) /= 0.0 ) then

          write ( *, '(a,g14.6,a,i2,a)' ) 
     &      '  +', c(i+1), ' * H', i, '(x)'

        end if

      end do

      return
      end
      subroutine polynomial_laguerre_print ( d, c, label )

c*********************************************************************72
c
cc polynomial_laguerre_print() prints a polynomial in the Laguerre basis.
c
c  Discussion:
c
c    The form of a Laguerre polynomial is:
c
c      p(x) = c(0)*L0(x) + c(1)*L(1)(x) + ... + c(n)*L(n)(x)
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  Input:
c
c    integer D: the polynomial degree.
c
c    real C(1:D+1): the polynomial coefficients, constant first.
c
c    character LABEL: an optional title.
c
      implicit none

      integer d

      double precision c(d+1)
      integer i
      character ( len = * ) label

      write ( *, '(a)' ) ''
      if ( 0 < len_trim ( label ) ) then
        write ( *, '(a)' ) trim ( label ) // ' = '
      end if

      if ( d < 0 ) then
        write ( *, '(a)' ) '  Zero polynomial'
        return
      end if

      if ( all ( c == 0 ) ) then
        write ( *, '(a)' ) '  0'
        return
      end if

      do i = d, 0, -1

        if ( c(i+1) /= 0.0 ) then

          write ( *, '(a,g14.6,a,i2,a)' ) 
     &      '  +', c(i+1), ' * L', i, '(x)'

        end if

      end do

      return
      end
      subroutine polynomial_legendre_print ( d, c, label )

c*********************************************************************72
c
cc polynomial_legendre_print() prints a polynomial in the Legendre basis.
c
c  Discussion:
c
c    The form of a Legendre polynomial is:
c
c      p(x) = c(0)*P0(x) + c(1)*P(1)(x) + ... + c(n)*P(n)(x)
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  Input:
c
c    integer D: the polynomial degree.
c
c    real C(1:D+1): the polynomial coefficients, constant first.
c
c    character LABEL: an optional title.
c
      implicit none

      integer d

      double precision c(d+1)
      integer i
      character ( len = * ) label

      write ( *, '(a)' ) ''
      if ( 0 < len_trim ( label ) ) then
        write ( *, '(a)' ) trim ( label ) // ' = '
      end if

      if ( d < 0 ) then
        write ( *, '(a)' ) '  Zero polynomial'
        return
      end if

      if ( all ( c == 0 ) ) then
        write ( *, '(a)' ) '  0'
        return
      end if

      do i = d, 0, -1

        if ( c(i+1) /= 0.0 ) then

          write ( *, '(a,g14.6,a,i2,a)' ) 
     &      '  +', c(i+1), ' * P', i, '(x)'

        end if

      end do

      return
      end
      subroutine polynomial_monomial_print ( d, c, label )

c*********************************************************************72
c
cc polynomial_monomial_print() prints a polynomial in monomial basis.
c
c  Discussion:
c
c    The monomial form of a polynomial is:
c
c      p(x) = c(0)*x^0 + c(1)*x + ... + c(n)*x^(n)
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  Input:
c
c    integer D: the polynomial degree.
c
c    real C(1:D+1): the polynomial coefficients, constant first.
c
c    character LABEL: an optional title.
c
      implicit none

      integer d

      double precision c(d+1)
      integer i
      character ( len = * ) label

      write ( *, '(a)' ) ''
      if ( 0 < len_trim ( label ) ) then
        write ( *, '(a)' ) trim ( label ) // ' = '
      end if

      if ( d < 0 ) then
        write ( *, '(a)' ) '  Zero polynomial'
        return
      end if

      if ( all ( c == 0 ) ) then
        write ( *, '(a)' ) '  0'
        return
      end if

      do i = d, 0, -1

        if ( c(i+1) /= 0.0 ) then

          if ( 2 <= i ) then
            write ( *, '(a,g14.6,a,i2)' ) '  +', c(i+1), ' * x^', i
          else if ( i == 1 ) then
            write ( *, '(a,g14.6,a)' ) '  +', c(i+1), ' * x'
          else if ( i == 0 ) then
            write ( *, '(a,g14.6)' ) '  +', c(i+1)
          end if

        end if

      end do

      return
      end
      subroutine c8vec_print ( n, a, title )

c*********************************************************************72
c
cc c8vec_print() prints a C8VEC, with an optional title.
c
c  Discussion:
c
c    A C8VEC is a vector of C8's
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    07 December 2008
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of components of the vector.
c
c    double complex A(N), the vector to be printed.
c
c    character*(*) TITLE, a title.
c
      implicit none

      integer n

      double complex a(n)
      integer i
      character*(*) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )
      write ( *, '(a)' ) ' '
      do i = 1, n
        write ( *, '(2x,i8,a,1x,2g14.6)' ) i, ':', a(i)
      end do

      return
      end
      subroutine r8poly_roots ( n, p, r )

c*********************************************************************72
c
cc r8poly_roots() returns the roots of a polynomial.
c
c  Discussion:
c
c    The monomial form of a polynomial is:
c
c      p(x) = c(0)*x^0 + c(1)*x + ... + c(n)*x^(n)
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  Input:
c
c    integer n: the polynomial degree.
c
c    double precision p(n+1): the polynomial coefficients, constant first.
c
c  Output:
c
c    double complex r(n): the polynomial roots.
c
      implicit none

      logical fail
      integer i
      integer n
      double precision op(n+1)
      double precision p(n+1)
      double complex r(n)
      double precision zeroi(n)
      double precision zeror(n)
      
      op(1:n+1) = p(n+1:1:-1)

      call rpoly ( op, n, zeror, zeroi, fail )

      if ( fail ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'r8poly_roots(): Fatal error!'
        write ( *, '(a)' ) '  rpoly() from toms493() failed.'
        stop 1
      end if

      do i = 1, n
        r(i) = dcmplx ( zeror(i), zeroi(i) )
      end do
      
      return
      end
      subroutine r8mat_print ( m, n, a, title )

c*********************************************************************72
c
cc r8mat_print() prints an R8MAT.
c
c  Discussion:
c
c    An R8MAT is an array of R8's.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    20 May 2004
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer M, the number of rows in A.
c
c    integer N, the number of columns in A.
c
c    double precision A(M,N), the matrix.
c
c    character * ( * ) TITLE, a title.
c
      implicit none

      integer m
      integer n

      double precision a(m,n)
      character * ( * ) title

      call r8mat_print_some ( m, n, a, 1, 1, m, n, title )

      return
      end
      subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi,
     &  title )

c*********************************************************************72
c
cc r8mat_print_some() prints some of an R8MAT.
c
c  Discussion:
c
c    An R8MAT is an array of R8's.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer M, N, the number of rows and columns.
c
c    double precision A(M,N), an M by N matrix to be printed.
c
c    integer ILO, JLO, the first row and column to print.
c
c    integer IHI, JHI, the last row and column to print.
c
c    character * ( * ) TITLE, a title.
c
      implicit none

      integer incx
      parameter ( incx = 5 )
      integer m
      integer n

      double precision a(m,n)
      character * ( 14 ) ctemp(incx)
      integer i
      integer i2hi
      integer i2lo
      integer ihi
      integer ilo
      integer inc
      integer j
      integer j2
      integer j2hi
      integer j2lo
      integer jhi
      integer jlo
      character * ( * ) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )

      if ( m .le. 0 .or. n .le. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  (None)'
        return
      end if

      do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

        j2hi = j2lo + incx - 1
        j2hi = min ( j2hi, n )
        j2hi = min ( j2hi, jhi )

        inc = j2hi + 1 - j2lo

        write ( *, '(a)' ) ' '

        do j = j2lo, j2hi
          j2 = j + 1 - j2lo
          write ( ctemp(j2), '(i7,7x)') j
        end do

        write ( *, '(''  Col   '',5a14)' ) ( ctemp(j), j = 1, inc )
        write ( *, '(a)' ) '  Row'
        write ( *, '(a)' ) ' '

        i2lo = max ( ilo, 1 )
        i2hi = min ( ihi, m )

        do i = i2lo, i2hi

          do j2 = 1, inc

            j = j2lo - 1 + j2

            write ( ctemp(j2), '(g14.6)' ) a(i,j)

          end do

          write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc )

        end do

      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end