subroutine bernstein_to_legendre01 ( n, bcoef, lcoef )

c*********************************************************************72
c
cc bernstein_to_legendre01() converts from Bernstein to Legendre01 form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision bcoef(0:n): the Bernstein coefficients of the polynomial.
c 
c  Output:
c
c    double precision lcoef(0:n): the Legendre01 coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision bcoef(0:n)
      double precision lcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )

      call bernstein_to_legendre01_matrix ( n, A )
c
c  Apply the linear transformation.
c
      lcoef = matmul ( A, bcoef )

      deallocate ( A )

      return
      end
      subroutine bernstein_to_legendre01_matrix ( n, A )

c*********************************************************************72
c
cc bernstein_to_legendre01_matrix() returns the Bernstein-to-Legendre matrix.
c
c  Discussion:
c
c    The Legendre polynomials are often defined on [-1,+1], while the
c    Bernstein polynomials are defined on [0,1].  For this function,
c    the Legendre polynomials have been shifted to share the [0,1]
c    interval of definition.
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 maximum degree of the polynomials.
c
c  Output:
c
c    double precision A(0:N,0:N), the matrix.
c
      implicit none

      integer n

      double precision A(0:n,0:n)
      integer i
      integer i4_choose
      integer j
      integer k

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

      do i = 0, n
        do j = 0, n
          do k = 0, i
            A(i,j) = A(i,j) 
     &        + ( -1.0D+00 ) ** ( i + k ) * i4_choose ( i, k ) ** 2 
     &        / i4_choose ( n + i, j + k )
          end do
          A(i,j) = A(i,j) * i4_choose ( n, j ) 
     &      * dble ( 2 * i + 1 ) / dble ( n + i + 1 )
        end do
      end do

      return
      end
      subroutine bernstein_to_monomial ( n, bcoef, mcoef )

c*********************************************************************72
c
cc bernstein_to_monomial() converts from Bernstein to monomial form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 April 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision bcoef(0:n): the Bernstein coefficients of the polynomial.
c 
c  Output:
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision bcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )

      call bernstein_to_monomial_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      mcoef = matmul ( A, bcoef )

      deallocate ( A )

      return
      end
      subroutine bernstein_to_monomial_matrix ( n, A )

c*********************************************************************72
c
cc bernstein_to_monomial_matrix() converts Bernstein to monomial form.
c
c  Discussion:
c
c    The Bernstein matrix of order N is an NxN matrix A which can be used to
c    transform a vector of power basis coefficients C representing a polynomial 
c    P(X) to a corresponding Bernstein basis coefficient vector B:
c
c      B = A * C
c
c    The N power basis vectors are ordered as (1,X,X^2,...X^(N-1)) and the N 
c    Bernstein basis vectors as ((1-X)^(N-1), X*(1_X)^(N-2),...,X^(N-1)).
c
c  Example:
c
c    N = 5
c
c    1    -4     6    -4     1
c    0     4   -12    12    -4
c    0     0     6   -12     6
c    0     0     0     4    -4
c    0     0     0     0     1
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 i0
      integer i4_choose
      integer j0
      integer n0

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

      n0 = n - 1

      do j0 = 0, n0
        do i0 = 0, j0
          A(i0+1,j0+1) = ( -1 ) ** ( j0 - i0 ) 
     &      * i4_choose ( n0 - i0, j0 - i0 ) 
     &      * i4_choose ( n0, i0 )
        end do
      end do

      return
      end
      subroutine chebyshev_to_monomial ( n, ccoef, mcoef )

c*********************************************************************72
c
cc chebyshev_to_monomial() converts from Chebyshev to monomial form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    18 February 2024
c
c  Author:
c
c    Original Fortran77 version by Fred Krogh.
c    This version by John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision ccoef(0:n): the Chebyshev coefficients of the polynomial.
c 
c  Output:
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c
      implicit none

      integer n

      double precision ccoef(0:n)
      integer i
      integer j
      double precision mcoef(0:n)
      double precision tp
c
      if ( n < 0 ) then
        return
      end if

      tp = 1.0D+00

      mcoef(0:n) = ccoef(0:n)

      do j = 0, n - 2
        do i = n - 2, j, -1
          mcoef(i) = mcoef(i) - mcoef(i+2)
        end do
        mcoef(j+1) = 0.5D+00 * mcoef(j+1)
        mcoef(j) = tp * mcoef(j)
        tp = 2.0D+00 * tp
      end do

      mcoef(n) = tp * mcoef(n)
      if ( 0 < n ) then
        mcoef(n-1) = tp * mcoef(n-1)
      end if

      return
      end
      subroutine chebyshev_to_monomial_matrix ( n, A )

c*********************************************************************72
c
cc chebyshev_to_monomial_matrix() converts from Chebyshev to monomial form.
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 ( n, alpha, gcoef, mcoef )

c*********************************************************************72
c
cc gegenbauer_to_monomial() converts from Gegenbauer to monomial form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision alpha: the Gegenbauer parameter.
c
c    double precision gcoef(0:n): the Gegenbauer coefficients of the polynomial.
c 
c  Output:
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision alpha
      double precision gcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )

      call gegenbauer_to_monomial_matrix ( n + 1, alpha, A )
c
c  Apply the linear transformation.
c
      mcoef = matmul ( A, gcoef )

      deallocate ( A )

      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    16 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 hermite_to_monomial ( n, hcoef, mcoef )

c*********************************************************************72
c
cc hermite_to_monomial() converts from Hermite to monomial form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision hcoef(0:n): the Hermite coefficients of the polynomial.
c 
c  Output:
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision hcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )

      call hermite_to_monomial_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      mcoef = matmul ( A, hcoef )

      deallocate ( A )

      return
      end
      subroutine hermite_to_monomial_matrix ( n, a )

c*********************************************************************72
c
cc hermite_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  Parameters:
c
c    Input, integer N, the order of the matrix.
c
c    Output, 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
      function i4_choose ( n, k )

c*********************************************************************72
c
cc i4_choose() computes the binomial coefficient C(N,K).
c
c  Discussion:
c
c    The value is calculated in such a way as to avoid overflow and
c    roundoff.  The calculation is done in integer arithmetic.
c
c    The formula used is:
c
c      C(N,K) = Nc / ( Kc * (N-K)c )
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 October 2014
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    ML Wolfson, HV Wright,
c    Algorithm 160:
c    Combinatorial of M Things Taken N at a Time,
c    Communications of the ACM,
c    Volume 6, Number 4, April 1963, page 161.
c
c  Input:
c
c    integer N, K, are the values of N and K.
c
c  Output:
c
c    integer I4_CHOOSE, the number of combinations of N
c    things taken K at a time.
c
      implicit none

      integer i
      integer i4_choose
      integer k
      integer mn
      integer mx
      integer n
      integer value

      mn = min ( k, n - k )
      mx = max ( k, n - k )

      if ( mn .lt. 0 ) then

        value = 0

      else if ( mn .eq. 0 ) then

        value = 1

      else

        value = mx + 1

        do i = 2, mn
          value = ( value * ( mx + i ) ) / i
        end do

      end if

      i4_choose = value

      return
      end
      subroutine laguerre_to_monomial ( n, lcoef, mcoef )

c*********************************************************************72
c
cc laguerre_to_monomial() converts from Laguerre to monomial form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision lcoef(0:n): the Laguerre coefficients of the polynomial.
c 
c  Output:
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision lcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )

      call laguerre_to_monomial_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      mcoef = matmul ( A, lcoef )

      deallocate ( A )

      return
      end
      subroutine laguerre_to_monomial_matrix ( n, a )

c*********************************************************************72
c
cc laguerre_to_monomial_matrix() converts from Laguerre to monomial basis.
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 ( n, lcoef, mcoef )

c*********************************************************************72
c
cc legendre_to_monomial() converts from Legendre to monomial form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision lcoef(0:n): the Legendre coefficients of the polynomial.
c 
c  Output:
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision lcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )

      call legendre_to_monomial_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      mcoef = matmul ( A, lcoef )

      deallocate ( A )

      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 legendre01_to_bernstein ( n, lcoef, bcoef )

c*********************************************************************72
c
cc legendre01_to_bernstein() converts from Legendre01 to Bernstein form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision lcoef(1:n+1): the Legendre01 coefficients of the polynomial.
c 
c  Output:
c
c    double precision bcoef(1:n+1): the Bernstein coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision bcoef(0:n)
      double precision lcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )
      call legendre01_to_bernstein_matrix ( n, A )
c
c  Apply the linear transformation.
c
      bcoef = matmul ( A, lcoef )

      deallocate ( A )

      return
      end
      subroutine legendre01_to_bernstein_matrix ( n, a )

c*********************************************************************72
c
cc legendre01_to_bernstein_matrix() returns the Legendre01-to-Bernstein matrix.
c
c  Discussion:
c
c    The Legendre polynomials are often defined on [-1,+1], while the
c    Bernstein polynomials are defined on [0,1].  For this function,
c    the Legendre polynomials have been shifted to share the [0,1]
c    interval of definition.
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 maximum degree of the polynomials.
c
c  Output:
c
c    double precision A(0:N,0:N), the matrix.
c
      implicit none

      integer n

      double precision A(0:n,0:n)
      integer i
      integer i4_choose
      integer j
      integer k

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

      do i = 0, n
        do j = 0, n
          do k = max ( 0, i + j - n ), min ( i, j )
            A(i,j) = A(i,j) 
     &        + ( -1.0D+00 ) ** ( j + k ) * i4_choose ( j, k ) ** 2 
     &        * i4_choose ( n - j, i - k )
          end do
          A(i,j) = A(i,j) / i4_choose ( n, i )
        end do
      end do

      return
      end
      subroutine monomial_to_bernstein ( n, mcoef, bcoef )

c*********************************************************************72
c
cc monomial_to_bernstein() converts from monomial to Bernstein form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision mcoef(1:n+1): the monomial coefficients of the polynomial.
c 
c  Output:
c
c    real bcoef(1:n+1): the Bernstein coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision bcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )
      call monomial_to_bernstein_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      bcoef = matmul ( A, mcoef )

      deallocate ( A )

      return
      end
      subroutine monomial_to_bernstein_matrix ( n, A )

c*********************************************************************72
c
cc monomial_to_bernstein_matrix() converts monomial to Bernsteim form.
c
c  Discussion:
c
c    The inverse Bernstein matrix of order N is an NxN matrix A which can 
c    be used to transform a vector of Bernstein basis coefficients B
c    representing a polynomial P(X) to a corresponding power basis 
c    coefficient vector C:
c
c      C = A * B
c
c    The N power basis vectors are ordered as (1,X,X^2,...X^(N-1)) and the N 
c    Bernstein basis vectors as ((1-X)^(N-1), X*(1-X)^(N-2),...,X^(N-1)).
c
c  Example:
c
c    N = 5
c
c   1.0000    1.0000    1.0000    1.0000    1.0000
c        0    0.2500    0.5000    0.7500    1.0000
c        0         0    0.1667    0.5000    1.0000
c        0         0         0    0.2500    1.0000
c        0         0         0         0    1.0000
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 inverse matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      integer i0
      integer i4_choose
      integer j0
      integer n0

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

      n0 = n - 1

      do j0 = 0, n0
        do i0 = 0, j0
          A(i0+1,j0+1) = dble ( i4_choose ( j0, i0 ) ) 
     &                 / dble ( i4_choose ( n0, i0 ) )
        end do
      end do

      return
      end
      subroutine monomial_to_chebyshev ( n, mcoef, ccoef )

c*********************************************************************72
c
cc monomial_to_chebyshev() converts from monomial to Chebyshev form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    18 February 2024
c
c  Author:
c
c    Original Fortran77 version by Fred Krogh.
c    This version by John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision mcoef(0:n): the monomial coefficients of the polynomial.
c 
c  Output:
c
c    double precision ccoef(0:n): the Chebyshev coefficients of the polynomial.
c
      implicit none

      integer n

      double precision ccoef(0:n)
      integer i
      integer j
      double precision mcoef(0:n)
      double precision tp

      if ( n < 0 ) then
        return
      end if

      ccoef(0:n) = mcoef(0:n)

      tp = 0.5D+00 ** ( n - 1 )
      ccoef(n) = tp * ccoef(n)

      if ( n == 0 ) then
        return
      end if

      ccoef(n-1) = tp * ccoef(n-1)
      do j = n - 2, 0, -1
        tp = 2.0D+00 * tp
        ccoef(j) = tp * ccoef(j)
        ccoef(j+1) = 2.0D+00 * ccoef(j+1)
        do i = j, n - 2
          ccoef(i) = ccoef(i) + ccoef(i+2)
        end do
      end do

      return
      end
      subroutine monomial_to_chebyshev_matrix ( n, A )

c*********************************************************************72
c
cc monomial_to_chebyshev_matrix() converts from monomial to Chebyshev T form.
c
c  Example:
c
c    N = 11
c    Each column must be divided by the divisor below it.
c
c      1   .   1   .   3   .  10   .   35    .  126
c      .   1   .   3   .  10   .  35    .  126    .
c      .   .   1   .   4   .  15   .   56    .  210
c      .   .   .   1   .   5   .  21    .   84    .
c      .   .   .   .   1   .   6   .   28    .  120
c      .   .   .   .   .   1   .   7    .   36    .
c      .   .   .   .   .   .   1   .    8    .   45
c      .   .   .   .   .   .   .   1    .    9    .
c      .   .   .   .   .   .   .   .    1    .   10
c      .   .   .   .   .   .   .   .    .    1    .
c      .   .   .   .   .   .   .   .    .    .    1  
c     /1  /1  /2  /4  /8 /16 /32 /64 /128 /256 /512 
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+1,j-1)   / 2.0D+00
          else if ( i .eq. 2 ) then
            A(i,j) = ( 2.0D+00 * A(i-1,j-1) + A(i+1,j-1) ) / 2.0D+00
          else if ( i .lt. n ) then
            A(i,j) = (           A(i-1,j-1) + A(i+1,j-1) ) / 2.0D+00
          else
            A(i,j) =             A(i-1,j-1)                / 2.0D+00
          end if
        end do
      end do

      return
      end
      subroutine monomial_to_gegenbauer ( n, alpha, mcoef, gcoef )

c*********************************************************************72
c
cc monomial_to_gegenbauer() converts from monomial to Gegenbauer form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision alpha: the Gegenbauer parameter.
c
c    double precision mcoef(1:n+1): the monomial coefficients of the polynomial.
c 
c  Output:
c
c    double precision gcoef(1:n+1): the Gegenbauer coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision alpha
      double precision gcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )
      call monomial_to_gegenbauer_matrix ( n + 1, alpha, A )
c
c  Apply the linear transformation.
c
      gcoef = matmul ( A, mcoef )

      deallocate ( A )

      return
      end
      subroutine monomial_to_gegenbauer_matrix ( n, alpha, A )

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

      integer n

      double precision A(0:n-1,0:n-1)
      double precision alpha
      double precision bot
      integer i
      integer ilo
      integer j
      double precision r8_factorial
      double precision top

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

      if ( n <= 0 ) then
        return
      end if

      do j = 0, n - 1

        ilo = mod ( j, 2 )

        do i = ilo, j, 2

          top = ( i + alpha ) * r8_factorial ( j ) * gamma ( alpha )

          bot = 2.0D+00 ** j * r8_factorial ( ( j - i ) / 2 ) 
     &      * gamma ( ( j + i ) / 2 + alpha + 1.0D+00 )

          A(i,j) = top / bot

        end do
      end do

      return
      end
      subroutine monomial_to_hermite ( n, mcoef, hcoef )

c*********************************************************************72
c
cc monomial_to_hermite() converts from monomial to Hermite form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision mcoef(1:n+1): the monomial coefficients of the polynomial.
c 
c  Output:
c
c    double precision hcoef(1:n+1): the Hermite coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision hcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )
      call monomial_to_hermite_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      hcoef = matmul ( A, mcoef )

      deallocate ( A )

      return
      end
      subroutine monomial_to_hermite_matrix ( n, a )

c*********************************************************************72
c
cc monomial_to_hermite_matrix() converts from monomial to Hermite basis.
c
c  Example:
c
c    N = 11 (Each column must be divided by the factor below it).
c
c    1     .     2     .    12     .   120     .  1680     . 30240
c    .     1     .     6     .    60     .   840     . 15120     .
c    .     .     1     .    12     .   180     .  3360     . 75600
c    .     .     .     1     .    20     .   420     . 10080     .
c    .     .     .     .     1     .    30     .   840     . 25200
c    .     .     .     .     .     1     .    42     .  1512     .
c    .     .     .     .     .     .     1     .    56     .  2520
c    .     .     .     .     .     .     .     1     .    72     .
c    .     .     .     .     .     .     .     .     1     .    90
c    .     .     .     .     .     .     .     .     .     1     .
c    .     .     .     .     .     .     .     .     .     .     1
c
c   /1    /2    /4    /8   /16   /32   /64  /128  /256  /512 /1024
c
c  Properties:
c
c    A is generally not symmetric: A' /= A.
c
c    A is nonnegative.
c
c    A is lower triangular.
c
c    det ( A ) = 1 / 2^((N*(N-1))/2)
c
c    LAMBDA(I) = 1 / 2^(I-1)
c
c    A is reducible.
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) = 0.5D+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 ( j - 2 ) * A(i,j-2)              ) / 2.0D+00
          else
            A(i,j) = 
     &        ( dble ( j - 2 ) * A(i,j-2) + A(i-1,j-1) ) / 2.0D+00
          end if
        end do
      end do

      return
      end
      subroutine monomial_to_laguerre ( n, mcoef, lcoef )

c*********************************************************************72
c
cc monomial_to_laguerre() converts from monomial to Laguerre form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision mcoef(1:n+1): the monomial coefficients of the polynomial.
c 
c  Output:
c
c    double precision lcoef(1:n+1): the Laguerre coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision lcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )
      call monomial_to_laguerre_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      lcoef = matmul ( A, mcoef )

      deallocate ( A )

      return
      end
      subroutine monomial_to_laguerre_matrix ( n, a )

c*********************************************************************72
c
cc monomial_to_laguerre_matrix() converts from monomial to Laguerre basis.
c
c  Example:
c
c    N = 9
c
c        1       1       2       6      24     120     720    5040    40320
c        .      -1      -4     -18     -96    -600   -4320  -35280  -322560
c        .       .       2      18     144    1200   10800  105840  1128960
c        .       .       .      -6     -96   -1200  -14400 -176400 -2257920
c        .       .       .       .      24     600   10800  176400  2822400
c        .       .       .       .       .    -120   -4320 -105840 -2257920
c        .       .       .       .       .       .     720   35280  1128960
c        .       .       .       .       .       .       .   -5040  -322560
c        .       .       .       .       .       .       .       .    40320
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) = (I-1)c
c    A(I,I) = (I-1)c * ( -1 )^(N+1)
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(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 ( j - 1 ) * ( A(i,j-1)              )
          else
            A(i,j) = dble ( j - 1 ) * ( A(i,j-1) - A(i-1,j-1) )
          end if
        end do
      end do

      return
      end
      subroutine monomial_to_legendre ( n, mcoef, lcoef )

c*********************************************************************72
c
cc monomial_to_legendre() converts from monomial to Legendre form.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    26 February 2024
c
c  Author:
c
c    John Burkardt.
c
c  Input:
c
c    integer n: the order of the polynomial.
c
c    double precision mcoef(1:n+1): the monomial coefficients of the polynomial.
c 
c  Output:
c
c    double precision lcoef(1:n+1): the Legendre coefficients of the polynomial.
c
      implicit none

      integer n

      double precision, allocatable :: A(:,:)
      double precision lcoef(0:n)
      double precision mcoef(0:n)

      if ( n < 0 ) then
        return
      end if
c
c  Get the matrix.
c
      allocate ( A(0:n,0:n) )
      call monomial_to_legendre_matrix ( n + 1, A )
c
c  Apply the linear transformation.
c
      lcoef = matmul ( A, mcoef )

      deallocate ( A )

      return
      end
      subroutine monomial_to_legendre_matrix ( n, a )

c*********************************************************************72
c
cc monomial_to_legendre_matrix() returns the inverse of the LEGENDRE matrix.
c
c  Example:
c
c    N = 11
c
c       1    .     .    .     .    .    .   .    .   .   .
c       .    1     .    .     .    .    .   .    .   .   .
c       1    .     2    .     .    .    .   .    .   .   . /     3
c       .    3     .    2     .    .    .   .    .   .   . /     5
c       7    .    20    .     8    .    .   .    .   .   . /    35
c       .   27     .   28     .    8    .   .    .   .   . /    63
c      33    .   110    .    72    .   16   .    .   .   . /   231
c       .  143     .  182     .   88    .  16    .   .   . /   429
c     715    .  2600    .  2160    .  832   .  128   .   . /  6435
c       . 3315     . 4760     . 2992    . 960    . 128   . / 12155
c    4199    . 16150    . 15504    . 7904   . 2176   . 256 / 46189
c
c  Properties:
c
c    A is nonnegative.
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 lower triangular.
c
c    A is reducible.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    20 June 2011
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 i = 3, n
        do j = 1, n
          if ( j .eq. 1 ) then

            a(i,j) = dble ( j ) * a(i-1,j+1) 
     &             / dble ( 2 * j + 1 )

          else if ( j .lt. n ) then

            a(i,j) = dble (     j - 1 ) * a(i-1,j-1) 
     &             / dble ( 2 * j - 3 ) 
     &             + dble (     j     ) * a(i-1,j+1) 
     &             / dble ( 2 * j + 1 )

          else

            a(i,j) = dble ( j - 1 ) * a(i-1,j-1) 
     &             / dble ( 2 * j - 3 )

          end if
        end do
      end do

      return
      end
      function r8_factorial ( n )

c*********************************************************************72
c
cc r8_factorial() computes the factorial of N.
c
c  Discussion:
c
c    factorial ( N ) = product ( 1 <= I <= N ) I
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    07 June 2008
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the argument of the factorial function.
c    If N is less than 1, the function value is returned as 1.
c
c  Output:
c
c    double precision R8_FACTORIAL, the factorial of N.
c
      implicit none

      integer i
      integer n
      double precision r8_factorial

      r8_factorial = 1.0D+00

      do i = 1, n
        r8_factorial = r8_factorial * dble ( i )
      end do

      return
      end