# include <cmath>
# include <cstdlib>
# include <iostream>
# include <iomanip>

using namespace std;

# include "polynomial_conversion.hpp"

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

double *bernstein_to_legendre01 ( int n, double bcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    bernstein_to_legendre01() converts from Bernstein to Legendre01 form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double bcoef(1:n+1): the Bernstein coefficients of the polynomial.
//
//  Output:
//
//    double lcoef(1:n+1): the Legendre01 coefficients of the polynomial.
//
{
  double *A;
  double *lcoef;

  if ( n < 0 )
  {
    return NULL;
  }
/*
  Get the matrix.
*/
  A = bernstein_to_legendre01_matrix ( n );
/*
  Apply the transformation.
*/
  lcoef = r8mat_mv_new ( n + 1, n + 1, A, bcoef );

  delete [] A;

  return lcoef;
}
//****************************************************************************80

double *bernstein_to_legendre01_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    bernstein_to_legendre01_matrix() returns the Bernstein-to-Legendre01 matrix.
//
//  Discussion:
//
//    The Legendre polynomials are often defined on [-1,+1], while the
//    Bernstein polynomials are defined on [0,1].  For this function,
//    the Legendre polynomials have been shifted to share the [0,1]
//    interval of definition.
//
//  Example:
//
//    bernstein_to_legendre01_matrix ( 5 ):
//
//    0.1667    0.1667    0.1667    0.1667    0.1667    0.1667
//   -0.3571   -0.2143   -0.0714    0.0714    0.2143    0.3571
//    0.2976   -0.0595   -0.2381   -0.2381   -0.0595    0.2976
//   -0.1389    0.1944    0.1111   -0.1111   -0.1944    0.1389
//    0.0357   -0.1071    0.0714    0.0714   -0.1071    0.0357
//   -0.0040    0.0198   -0.0397    0.0397   -0.0198    0.0040
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the maximum degree of the polynomials.
//
//  Output:
//
//    double A[(N+1)*(N+1)]: the matrix.
//
{
  double *A;
  int i;
  int j;
  int k;

  A = new double[(n+1)*(n+1)];

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

  for ( i = 0; i <= n; i++ )
  {
    for ( j = 0; j <= n; j++ )
    {
      for ( k = 0; k <= i; k++ )
      {
        A[i+j*(n+1)] = A[i+j*(n+1)]
          + r8_mop ( i + k ) * pow ( r8_choose ( i, k ), 2 ) 
          / r8_choose ( n + i, j + k );
      }
      A[i+j*(n+1)] = A[i+j*(n+1)] * r8_choose ( n, j ) 
        * ( double ) ( 2 * i + 1 ) / ( double ) ( n + i + 1 );
    }
  }

  return A;
}
//****************************************************************************80

double *bernstein_to_monomial ( int n, double bcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    bernstein_to_monomial() converts from Bernstein to monomial form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double bcoef(1:n+1): the Bernstein coefficients of the polynomial.
//
//  Output:
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
{
  double *A;
  double *mcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = bernstein_to_monomial_matrix ( n + 1 );
//
//  Apply the transformation.
//
  mcoef = r8mat_mv_new ( n + 1, n + 1, A, bcoef );

  delete [] A;

  return mcoef;
}
//****************************************************************************80

double *bernstein_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    bernstein_to_monomial_matrix() returns the Bernstein-to-monomial matrix.
//
//  Discussion:
//
//    The N power basis vectors are ordered as (1,X,X^2,...X^(N-1)) and the N 
//    Bernstein basis vectors as ((1-X)^(N-1), X*(1_X)^(N-2),...,X^(N-1)).
//
//    1    -4     6    -4     1
//    0     4   -12    12    -4
//    0     0     6   -12     6
//    0     0     0     4    -4
//    0     0     0     0     1
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the order of the matrix.
//
//  Output:
//
//    double bernstein_to_monomial_matrix[n*n], the matrix.
//
{
  double *A;
  int i;
  int j;

  A = new double[n*n];

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

  for ( j = 0; j <= n - 1; j++ )
  {
    for ( i = 0; i <= j; i++ )
    {
      A[i+j*n] = r8_mop ( j - i ) * r8_choose ( n - 1 - i, j - i ) 
        * r8_choose ( n - 1, i );
    }
  }

  return A;
}
//****************************************************************************80

double *chebyshev_to_monomial ( int n, double ccoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    chebyshev_to_monomial() converts from Chebyshev to monomial form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 February 2024
//
//  Author:
//
//    Original Fortran77 version by Fred Krogh.
//    This version by John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double ccoef(0:n): the Chebyshev coefficients of the polynomial.
//
//  Output:
//
//    double mcoef(0:n): the monomial coefficients of the polynomial.
//
{
  int i;
  int j;
  double *mcoef;
  double tp;

  if ( n < 0 )
  {
    return NULL;
  }

  mcoef = new double[n+1];

  for ( i = 0; i <= n; i++ )
  {
    mcoef[i] = ccoef[i];
  }

  tp = 1.0;

  for ( j = 0; j <= n - 2; j++ )
  {
    for ( i = n - 2; j <= i; i-- )
    {
      mcoef[i] = mcoef[i] - mcoef[i+2];
    }
    mcoef[j+1] = 0.5 * mcoef[j+1];
    mcoef[j] = tp * mcoef[j];
    tp = 2.0 * tp;
  }

  mcoef[n] = tp * mcoef[n];
  if ( 0 < n )
  {
    mcoef[n-1] = tp * mcoef[n-1];
  }

  return mcoef;
}
//****************************************************************************80

double *chebyshev_to_monomial_matrix ( int n )

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

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    i = 0;
    A[i+j*n] = - A[i+(j-2)*n];
    for ( i = 1; i < n; i++ )
    {
      A[i+j*n] = 2.0 * A[i-1+(j-1)*n] - A[i+(j-2)*n];
    }
  }
  return A;
}
//****************************************************************************80

double *gegenbauer_to_monomial ( int n, double alpha, double gcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_to_monomial() converts from Gegenbauer to monomial form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 April 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double alpha: the Gegenbauer parameter.
//
//    double gcoef(1:n+1): the Gegenbauer coefficients of the polynomial.
//
//  Output:
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
{
  double *A;
  double *mcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = gegenbauer_to_monomial_matrix ( n + 1, alpha );
//
//  Apply the transformation.
//
  mcoef = r8mat_mv_new ( n + 1, n + 1, A, gcoef );

  delete [] A;

  return mcoef;
}
//****************************************************************************80

double *gegenbauer_to_monomial_matrix ( int n, double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the order of A.
//
//    double alpha: the parameter.
//
//  Output:
//
//    double A[N,N]: the matrix.
//
{
  double *A;
  double c1;
  double c2;
  int i;
  int j;
  int nn;

  A = ( double * ) malloc ( n * n * sizeof ( double ) );

  if ( n <= 0 )
  {
    return A;
  }

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 2.0 * alpha;
//
//  Perform convex sum.
//  Translating "(n+1) C(n+1) = 2 (n+alpha) x C(n) - ( n + 2 alpha - 1 ) C(n-1)"
//  drove me nuts, between indexing at 1 rather than 0, and dealing with
//  the interpretation of "n+1", because I now face the rare "off by 2" error!
//
  for ( j = 2; j < n; j++ )
  {
    nn = j - 1;
    c1 = ( 2.0 * nn + 2.0 * alpha       ) / ( nn + 1 );
    c2 = (     - nn - 2.0 * alpha + 1.0 ) / ( nn + 1 );
    for ( i = 1; i <= j; i++ )
    {
      A[i+j*n] = c1 * A[i-1+(j-1)*n];
    }
    for ( i = 0; i <= j - 2; i++ )
    {
      A[i+j*n] = A[i+j*n] + c2 * A[i+(j-2)*n];
    }
  }

  return A;
}
//****************************************************************************80

double *hermite_to_monomial ( int n, double hcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    hermite_to_monomial() converts from Hermite to monomial form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double hcoef(1:n+1): the Hermite coefficients of the polynomial.
//
//  Output:
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
{
  double *A;
  double *mcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = hermite_to_monomial_matrix ( n + 1 );
//
//  Apply the transformation.
//
  mcoef = r8mat_mv_new ( n + 1, n + 1, A, hcoef );

  delete [] A;

  return mcoef;
}
//****************************************************************************80

double *hermite_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    hermite_to_monomial_matrix() converts from Hermite to monomial form.
//
//  Example:
//
//    N = 11
//
//      1   .   -2      .      12     .   -120     .   1680      .  -30240
//      .   2    .    -12       .   120      . -1680      .  30240       .
//      .   .    4      .     -48     .    720     . -13440      .  302400
//      .   .    .      8       .  -160      .  3360      . -80640       .
//      .   .    .      .      16     .   -480     .  13440      . -403200
//      .   .    .      .       .    32      . -1344      .  48384       .
//      .   .    .      .       .     .     64     .  -3584      .  161280
//      .   .    .      .       .     .      .   128      .  -9216       .
//      .   .    .      .       .     .      .     .    256      .  -23040
//      .   .    .      .       .     .      .     .      .    512       .
//      .   .    .      .       .     .      .     .      .      .    1024
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double A[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 2.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == 0 )
      {
        A[i+j*n] =                      
          - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n];
      }
      else
      {
        A[i+j*n] = 2.0 * A[i-1+(j-1)*n] 
          - 2.0 * ( double ) ( j - 1 ) * A[i+(j-2)*n];
      }
    }
  }

  return A;
}
//****************************************************************************80

double *laguerre_to_monomial ( int n, double lcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    laguerre_to_monomial() converts from Laguerre to monomial form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double lcoef(1:n+1): the Laguerre coefficients of the polynomial.
//
//  Output:
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
{
  double *A;
  double *mcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = laguerre_to_monomial_matrix ( n + 1 );
//
//  Apply the transformation.
//
  mcoef = r8mat_mv_new ( n + 1, n + 1, A, lcoef );

  delete [] A;

  return mcoef;
}
//****************************************************************************80

double *laguerre_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    laguerre_to_monomial_matrix() converts from Laguerre to monomial form.
//
//  Example:
//
//    N = 8 (each column must be divided by the factor below it.)
//
//      1      1      2      6     24    120    720   5040
//      .     -1     -4    -18    -96   -600  -4320 -35280
//      .      .      1      9     72    600   5400  52920
//      .      .      .      1    -16   -200  -2400 -29400
//      .      .      .      .      1     25    450   7350
//      .      .      .      .      .     -1    -36   -882
//      .      .      .      .      .      .      1     49
//      .      .      .      .      .      .      .     -1
//
//     /1     /1     /2     /6    /24   /120   /720  /5040
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    US Department of Commerce, 1964.
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double A[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[0+1*n] =   1.0;
  A[1+1*n] = - 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    A[0+j*n] = ( ( double ) ( 2 * j - 1 ) * A[0+(j-1)*n]
               + ( double ) (   - j + 1 ) * A[0+(j-2)*n] )
               / ( double ) (     j     );
    for ( i = 1; i < n; i++ )
    {
      A[i+j*n] = ( ( double ) ( 2 * j - 1 ) * A[i+(j-1)*n] 
                 - ( double ) (         1 ) * A[i-1+(j-1)*n]
                 + ( double ) (   - j + 1 ) * A[i+(j-2)*n] ) 
                 / ( double ) (     j     );
    }
  }
  return A;
}
//****************************************************************************80

double *legendre_to_monomial ( int n, double lcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    legendre_to_monomial() converts from Legendre to monomial form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double lcoef(1:n+1): the Legendre coefficients of the polynomial.
//
//  Output:
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
{
  double *A;
  double *mcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = legendre_to_monomial_matrix ( n + 1 );
//
//  Apply the transformation.
//
  mcoef = r8mat_mv_new ( n + 1, n + 1, A, lcoef );

  delete [] A;

  return mcoef;
}
//****************************************************************************80

double *legendre_to_monomial_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    legendre_to_monomial_matrix() converts from Legendre to monomial form.
//
//  Example:
//
//    N = 11 (each column must be divided by factor at bottom)
//
//     1    .    -1     .      3     .     -5      .      35     .   -63
//     .    1     .    -3      .    15      .    -25       .   315     .
//     .    .     3     .    -30     .    105      .   -1260     .  3465
//     .    .     .     5      .   -70      .    315       . -4620     .
//     .    .     .     .     35     .   -315      .    6930     .-30030
//     .    .     .     .      .    63      .   -693       . 18018     .
//     .    .     .     .      .     .    231      .  -12012     . 90090
//     .    .     .     .      .     .      .    429       .-25740     .
//     .    .     .     .      .     .      .      .    6435     -109395
//     .    .     .     .      .     .      .      .       . 12155     .
//     .    .     .     .      .     .      .      .       .     . 46189
//
//    /1   /1    /2    /2     /8    /8    /16    /16    /128  /128  /256
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double A[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 3; j <= n; j++ )
  {
    for ( i = 1; i <= n; i++ )
    {
      if ( i == 1 )
      {
        A[i-1+(j-1)*n] = - ( j - 2 ) * A[i-1+(j-3)*n] 
                         / ( j - 1 );
      }
      else
      {
        A[i-1+(j-1)*n] = ( ( 2 * j - 3 ) * A[i-2+(j-2)*n] 
                         + (   - j + 2 ) * A[i-1+(j-3)*n] ) 
                         / (     j - 1 );
      }
    }
  }

  return A;
}
//****************************************************************************80

double *legendre01_to_bernstein ( int n, double lcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    legendre01_to_bernsteinl() converts from Legendre01 to Bernstein form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double lcoef(1:n+1): the Legendre01 coefficients of the polynomial.
//
//  Output:
//
//    double bcoef(1:n+1): the Bernstein coefficients of the polynomial.
//
{
  double *A;
  double *bcoef;

  if ( n < 0 )
  {
    return NULL;
  }
/*
  Get the matrix.
*/
  A = legendre01_to_bernstein_matrix ( n );
/*
  Apply the transformation.
*/
  bcoef = r8mat_mv_new ( n + 1, n + 1, A, lcoef );

  delete [] A;

  return bcoef;
}
//****************************************************************************80

double *legendre01_to_bernstein_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    legendre01_to_bernstein_matrix() returns the Legendre-to-Bernstein matrix.
//
//  Discussion:
//
//    The Legendre polynomials are often defined on [-1,+1], while the
//    Bernstein polynomials are defined on [0,1].  For this function,
//    the Legendre polynomials have been shifted to share the [0,1]
//    interval of definition.
//
//  Example:
//
//    legendre01_to_bernstein_matrix(5):
//
//    1.0000   -1.0000    1.0000   -1.0000    1.0000   -1.0000
//    1.0000   -0.6000   -0.2000    1.4000   -3.0000    5.0000
//    1.0000   -0.2000   -0.8000    0.8000    2.0000  -10.0000
//    1.0000    0.2000   -0.8000   -0.8000    2.0000   10.0000
//    1.0000    0.6000   -0.2000   -1.4000   -3.0000   -5.0000
//    1.0000    1.0000    1.0000    1.0000    1.0000    1.0000
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the maximum degree of the polynomials.
//
//  Output:
//
//    double A[(N+1)*(N+1)], the matrix.
//
{
  double *A;
  int i;
  int j;
  int k;

  A = new double[(n+1)*(n+1)];

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

  for ( i = 0; i <= n; i++ )
  {
    for ( j = 0; j <= n; j++ )
    {
      for ( k = max ( 0, i + j - n ); k <= min ( i, j ); k++ )
      {
        A[i+j*(n+1)] = A[i+j*(n+1)]
          + r8_mop ( j + k ) * pow ( r8_choose ( j, k ), 2 )
          * r8_choose ( n - j, i - k );
      }
      A[i+j*(n+1)] = A[i+j*(n+1)] / r8_choose ( n, i );
    }
  }

  return A;
}
//****************************************************************************80

double *monomial_to_bernstein ( int n, double mcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_bernstein() converts from monomial to Bernstein form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
//  Output:
//
//    double bcoef(1:n+1): the Bernstein coefficients of the polynomial.
//
{
  double *A;
  double *bcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = monomial_to_bernstein_matrix ( n + 1 );
//
//  Apply the transformation.
//
  bcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef );

  delete [] A;

  return bcoef;
}
//****************************************************************************80

double *monomial_to_bernstein_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_bernstein_matrix() converts monomial to Bernstein form.
//
//  Discussion:
//
//    N = 5
//
//   1.0000    1.0000    1.0000    1.0000    1.0000
//        0    0.2500    0.5000    0.7500    1.0000
//        0         0    0.1667    0.5000    1.0000
//        0         0         0    0.2500    1.0000
//        0         0         0         0    1.0000
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the maximum degree of the polynomials.
//
//  Output:
//
//    double A: the Power-to-Bernstein matrix.
//
{
  double *A;
  int i;
  int j;

  A = new double[n*n];

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

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i <= j; i++ )
    {
      A[i+j*n] = r8_choose ( j, i ) / r8_choose ( n - 1, i );
    }
  }

  return A;
}
//****************************************************************************80

double *monomial_to_chebyshev ( int n, double mcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_chebyshev() converts from monomial to Chebyshev form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    19 February 2024
//
//  Author:
//
//    Original Fortran77 version by Fred Krogh.
//    This version by John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double mcoef(0:n): the monomial coefficients of the polynomial.
//
//  Output:
//
//    double ccoef(0:n): the Chebyshev coefficients of the polynomial.
//
{
  double *ccoef;
  int i;
  int j;
  double tp;

  if ( n < 0 )
  {
    return NULL;
  }

  ccoef = new double[n+1];

  for ( i = 0; i <= n; i++ )
  {
    ccoef[i] = mcoef[i];
  }

  tp = pow ( 0.5, n - 1 );
  ccoef[n] = tp * ccoef[n];

  if ( n == 0 )
  {
    return ccoef;
  }

  ccoef[n-1] = tp * ccoef[n-1];
  for ( j = n - 2; 0 <= j; j-- )
  {
    tp = 2.0 * tp;
    ccoef[j] = tp * ccoef[j];
    ccoef[j+1] = 2.0 * ccoef[j+1];
    for ( i = j; i <= n - 2; i++ )
    {
      ccoef[i] = ccoef[i] + ccoef[i+2];
    }
  }

  return ccoef;
}
//****************************************************************************80

double *monomial_to_chebyshev_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_chebyshev_matrix() maps from monomial to Chebyshev T form.
//
//  Example:
//
//    N = 11
//    Each column must be divided by the divisor below it.
//
//      1   .   1   .   3   .  10   .   35    .  126
//      .   1   .   3   .  10   .  35    .  126    .
//      .   .   1   .   4   .  15   .   56    .  210
//      .   .   .   1   .   5   .  21    .   84    .
//      .   .   .   .   1   .   6   .   28    .  120
//      .   .   .   .   .   1   .   7    .   36    .
//      .   .   .   .   .   .   1   .    8    .   45
//      .   .   .   .   .   .   .   1    .    9    .
//      .   .   .   .   .   .   .   .    1    .   10
//      .   .   .   .   .   .   .   .    .    1    .
//      .   .   .   .   .   .   .   .    .    .    1  
//     /1  /1  /2  /4  /8 /16 /32 /64 /128 /256 /512 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N: the order of the matrix.
//
//  Output:
//
//    double CHEBY_T_INVERSE[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == 0 )
      {
        A[i+j*n] =                          A[i+1+(j-1)*n]   / 2.0;
      }
      else if ( i == 1 )
      {
        A[i+j*n] = ( 2.0 * A[i-1+(j-1)*n] + A[i+1+(j-1)*n] ) / 2.0;
      }
      else if ( i < n - 1 )
      {
        A[i+j*n] = (       A[i-1+(j-1)*n] + A[i+1+(j-1)*n] ) / 2.0;
      }
      else
      {
        A[i+j*n] =         A[i-1+(j-1)*n]                    / 2.0;
      }
    }
  }

  return A;
}
//****************************************************************************80

double *monomial_to_gegenbauer ( int n, double alpha, double mcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_gegenbauer() converts from monomial to Gegenbauer form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 April 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double alpha: the Gegenbauer parameter.
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
//  Output:
//
//    double gcoef(1:n+1): the Gegenbauer coefficients of the polynomial.
//
{
  double *A;
  double *gcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = monomial_to_gegenbauer_matrix ( n + 1, alpha );
//
//  Apply the transformation.
//
  gcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef );

  delete [] A;

  return gcoef;
}
//****************************************************************************80

double *monomial_to_gegenbauer_matrix ( int n, double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_gegenbauer_matrix(): monomial to Gegenbauer conversion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the order of A.
//
//    double alpha: the parameter.
//
//  Output:
//
//    double A[N,N]: the matrix.
//
{
  double *A;
  double bot;
  int i;
  int ilo;
  int j;
  double top;

  A = new double[n*n];

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

  if ( n <= 0 )
  {
    return A;
  }

  for ( j = 0; j < n; j++ )
  {
    ilo = ( j % 2 );
    for ( i = ilo; i <= j; i = i + 2 )
    {
      top = ( i + alpha ) * r8_factorial ( j ) * tgamma ( alpha );
      bot = pow ( 2.0, j ) * r8_factorial ( ( j - i ) / 2 ) 
        * tgamma ( ( j + i ) / 2 + alpha + 1.0 );
      A[i+j*n] = top / bot;
    }
  }

  return A;
}
//****************************************************************************80

double *monomial_to_hermite ( int n, double mcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_hermite() converts from monomial to Hermite form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
//  Output:
//
//    double hcoef(1:n+1): the Hermite coefficients of the polynomial.
//
{
  double *A;
  double *hcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = monomial_to_hermite_matrix ( n + 1 );
//
//  Apply the transformation.
//
  hcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef );

  delete [] A;

  return hcoef;
}
//****************************************************************************80

double *monomial_to_hermite_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_hermite_matrix() converts from monomial to Hermite form.
//
//  Example:
//
//    N = 11 (each column must be divided by the factor below it)
//
//    1     .     2     .    12     .   120     .  1680     . 30240
//    .     1     .     6     .    60     .   840     . 15120     .
//    .     .     1     .    12     .   180     .  3360     . 75600
//    .     .     .     1     .    20     .   420     . 10080     .
//    .     .     .     .     1     .    30     .   840     . 25200
//    .     .     .     .     .     1     .    42     .  1512     .
//    .     .     .     .     .     .     1     .    56     .  2520
//    .     .     .     .     .     .     .     1     .    72     .
//    .     .     .     .     .     .     .     .     1     .    90
//    .     .     .     .     .     .     .     .     .     1     .
//    .     .     .     .     .     .     .     .     .     .     1
//
//   /1    /2    /4    /8   /16   /32   /64  /128  /256  /512 /1024
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    21 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double A[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 0.5;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == 0 )
      {
        A[i+j*n] = ( ( double ) ( j - 1 ) * A[i+(j-2)*n]                  ) / 2.0;
      }
      else
      {
        A[i+j*n] = ( ( double ) ( j - 1 ) * A[i+(j-2)*n] + A[i-1+(j-1)*n] ) / 2.0;
      }
    }
  }
  return A;
}
//****************************************************************************80

double *monomial_to_laguerre ( int n, double mcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_laguerre() converts from monomial to Laguerre form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
//  Output:
//
//    double lcoef(1:n+1): the Laguerre coefficients of the polynomial.
//
{
  double *A;
  double *lcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = monomial_to_laguerre_matrix ( n + 1 );
//
//  Apply the transformation.
//
  lcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef );

  delete [] A;

  return lcoef;
}
//****************************************************************************80

double *monomial_to_laguerre_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_laguerre_matrix() converts from monomial to Laguerre form.
//
//  Example:
//
//    N = 9
//
//        1       1       2       6      24     120     720    5040    40320
//        .      -1      -4     -18     -96    -600   -4320  -35280  -322560
//        .       .       2      18     144    1200   10800  105840  1128960
//        .       .       .      -6     -96   -1200  -14400 -176400 -2257920
//        .       .       .       .      24     600   10800  176400  2822400
//        .       .       .       .       .    -120   -4320 -105840 -2257920
//        .       .       .       .       .       .     720   35280  1128960
//        .       .       .       .       .       .       .   -5040  -322560
//        .       .       .       .       .       .       .       .    40320
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double A[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[0+1*n] =   1.0;
  A[1+1*n] = - 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j < n; j++ )
  {
    i = 0;
    A[i+j*n] = ( double ) ( j ) * A[i+(j-1)*n];
    for ( i = 1; i < n; i++ )
    {
      A[i+j*n] = ( double ) ( j ) * ( A[i+(j-1)*n] - A[i-1+(j-1)*n] );
    }
  }
  return A;
}
//****************************************************************************80

double *monomial_to_legendre ( int n, double mcoef[] )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_legendre() converts from monomial to Legendre form.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt.
//
//  Input:
//
//    int n: the order of the polynomial.
//
//    double mcoef(1:n+1): the monomial coefficients of the polynomial.
//
//  Output:
//
//    double lcoef(1:n+1): the Legendre coefficients of the polynomial.
//
{
  double *A;
  double *lcoef;

  if ( n < 0 )
  {
    return NULL;
  }
//
//  Get the matrix.
//
  A = monomial_to_legendre_matrix ( n + 1 );
//
//  Apply the transformation.
//
  lcoef = r8mat_mv_new ( n + 1, n + 1, A, mcoef );

  delete [] A;

  return lcoef;
}
//****************************************************************************80

double *monomial_to_legendre_matrix ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    monomial_to_legendre_matrix() converts from monomial to Legendre form.
//
//  Example:
//
//    N = 11 (each column must be divided by the underlying factor).
//
//       1     .      1     .      7     .    33    .   715    . 4199
//       .     1      .     3      .    27     .  143     . 3315    .
//       .     .      2     .     20     .   110    .  2600    .16150
//       .     .      .     2      .    28     .  182     . 4760    .
//       .     .      .     .      8     .    72    .  2160    .15504
//       .     .      .     .      .     8     .   88     . 2992    .
//       .     .      .     .      .     .    16    .   832    . 7904
//       .     .      .     .      .     .     .   16     .  960    .
//       .     .      .     .      .     .     .    .   128    . 2176
//       .     .      .     .      .     .     .    .     .  128    .
//       .     .      .     .      .     .     .    .     .    .  256
//
//      /1    /1     /3    /5    /35   /63  /231 /429 /6435/12155/46189  
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    20 February 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the order of the matrix.
//
//  Output:
//
//    double A[N*N], the matrix.
//
{
  double *A;
  int i;
  int j;

  if ( n <= 0 )
  {
    return NULL;
  }

  A = new double[n*n];

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 1.0;

  if ( n == 2 )
  {
    return A;
  }

  for ( j = 2; j <= n - 1; j++ )
  {
    for ( i = 0; i <= n - 1; i++ )
    {
      if ( i == 0 )
      {
        A[i+j*n] = (     i + 1 ) * A[i+1+(j-1)*n] / ( 2 * i + 3 );
      }
      else if ( i < n - 1 )
      {
        A[i+j*n] = (     i     ) * A[i-1+(j-1)*n] / ( 2 * i - 1 )
                 + (     i + 1 ) * A[i+1+(j-1)*n] / ( 2 * i + 3 );
      }
      else
      {
        A[i+j*n] = (     i     ) * A[i-1+(j-1)*n] / ( 2 * i - 1 );
      }
    }
  }

  return A;
}
//****************************************************************************80

double r8_choose ( int n, int k )

//****************************************************************************80
//
//  Purpose:
//
//    r8_choose() computes the binomial coefficient C(N,K) as an R8.
//
//  Discussion:
//
//    The value is calculated in such a way as to avoid overflow and
//    roundoff.  The calculation is done in R8 arithmetic.
//
//    The formula used is:
//
//      C(N,K) = N! / ( K! * (N-K)! )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 July 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    ML Wolfson, HV Wright,
//    Algorithm 160:
//    Combinatorial of M Things Taken N at a Time,
//    Communications of the ACM,
//    Volume 6, Number 4, April 1963, page 161.
//
//  Input:
//
//    int N, K, the values of N and K.
//
//  Output:
//
//    double R8_CHOOSE, the number of combinations of N
//    things taken K at a time.
//
{
  int i;
  int mn;
  int mx;
  double value;

  mn = min ( k, n - k );

  if ( mn < 0 )
  {
    value = 0.0;
  }
  else if ( mn == 0 )
  {
    value = 1.0;
  }
  else
  {
    mx = max ( k, n - k );
    value = ( double ) ( mx + 1 );

    for ( i = 2; i <= mn; i++ )
    {
      value = ( value * ( double ) ( mx + i ) ) / ( double ) i;
    }
  }
  return value;
}
//****************************************************************************80

double r8_factorial ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    r8_factorial() computes the factorial of N.
//
//  Discussion:
//
//    factorial ( N ) = product ( 1 <= I <= N ) I
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 January 1999
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the argument of the factorial function.
//    If N is less than 1, the function value is returned as 1.
//
//  Output:
//
//    double R8_FACTORIAL, the factorial of N.
//
{
  int i;
  double value;

  value = 1.0;

  for ( i = 1; i <= n; i++ )
  {
    value = value * ( double ) ( i );
  }

  return value;
}
//****************************************************************************80

double r8_mop ( int i )

//****************************************************************************80
//
//  Purpose:
//
//    r8_mop() returns the I-th power of -1 as an R8 value.
//
//  Discussion:
//
//    An R8 is an double value.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 November 2007
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int I, the power of -1.
//
//  Output:
//
//    double R8_MOP, the I-th power of -1.
//
{
  double value;

  if ( ( i % 2 ) == 0 )
  {
    value = 1.0;
  }
  else
  {
    value = -1.0;
  }

  return value;
}
//****************************************************************************80

double *r8mat_mv_new ( int m, int n, double a[], double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    r8mat_mv_new() multiplies a matrix times a vector.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    For this routine, the result is returned as the function value.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 April 2007
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int M, N, the number of rows and columns of the matrix.
//
//    double A[M,N], the M by N matrix.
//
//    double X[N], the vector to be multiplied by A.
//
//  Output:
//
//    double R8MAT_MV_NEW[M], the product A*X.
//
{
  int i;
  int j;
  double *y;

  y = new double[m];

  for ( i = 0; i < m; i++ )
  {
    y[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      y[i] = y[i] + a[i+j*m] * x[j];
    }
  }

  return y;
}