function companion_matrix_test ( )

%*****************************************************************************80
%
%% companion_matrix_test() tests companion_matrix().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 March 2024
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../companion_matrix' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_matrix_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  companion_matrix() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial, in various bases.\n' );

  companion_chebyshev_test ( );
  companion_gegenbauer_test ( );
  companion_hermite_test ( );
  companion_laguerre_test ( );
  companion_legendre_test ( );
  companion_monomial_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_matrix_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../companion_matrix' );

  return
end
function companion_chebyshev_test ( )

%*****************************************************************************80
%
%% companion_chebyshev_test() tests companion_chebyshev().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_chebyshev_test():\n' );
  fprintf ( 1, '  companion_chebyshev() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial p(x) in the Chebyshev basis.\n' );

  p = [ 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 ];

  polynomial_chebyshev_print ( p, '  p(x)' );
%
%  Compute monomial form so we can get roots directly.
%
  n = length ( p );
  A = chebyshev_to_monomial_matrix ( n );
  q = A * p';
  polynomial_monomial_print ( q, '  Monomial q(x)' );
  r1 = roots ( q(end:-1:1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Roots of q(x):\n' );
  disp ( r1 );

  A = companion_chebyshev ( p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Chebyshev companion matrix A(p):\n' );

  disp ( A );

  r = eigs ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Eigenvalues of A(p):\n' );
  disp ( r );

  return
end
function companion_gegenbauer_test ( )

%*****************************************************************************80
%
%% companion_gegenbauer_test() tests companion_gegenbauer().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_gegenbauer_test():\n' );
  fprintf ( 1, '  companion_gegenbauer() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial p(x) in the Gegenbauer basis.\n' );

  p = [ 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 ];

  polynomial_gegenbauer_print ( p, '  p(x)' );
%
%  Try several parameter values.
%
  for alpha = [ 0.5, 1.0 ]

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Gegenbauer parameter alpha will be %g\n', alpha );
%
%  Compute monomial form so we can get roots directly.
%
    n = length ( p );
    A = gegenbauer_to_monomial_matrix ( n, alpha );
    q = A * p';
    polynomial_monomial_print ( q, '  Monomial q(x)' );
    r1 = roots ( q(end:-1:1) );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Roots of q(x):\n' );
    disp ( r1 );

    A = companion_gegenbauer ( p, alpha );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Gegenbauer companion matrix A(p), alpha = %g:\n', alpha );

    disp ( A );

    r = eigs ( A );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Eigenvalues of A(p):\n' );
    disp ( r );

  end

  return
end
function companion_hermite_test ( )

%*****************************************************************************80
%
%% companion_hermite_test() tests companion_hermite().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_hermite_test():\n' );
  fprintf ( 1, '  companion_hermite() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial p(x) in the Hermite basis.\n' );

  p = [ 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 ];

  polynomial_hermite_print ( p, '  Hermite p(x)' );
%
%  Compute monomial form so we can get roots directly.
%
  n = length ( p );
  A = h_to_monomial_matrix ( n );
  q = A * p';
  polynomial_monomial_print ( q, '  Monomial q(x)' );
  r1 = roots ( q(end:-1:1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Roots of q(x):\n' );
  disp ( r1 );
 
  A = companion_hermite ( p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Hermite companion matrix A(p):\n' );

  disp ( A );

  r = eigs ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Eigenvalues of A(p):\n' );
  disp ( r );

  return
end
function companion_laguerre_test ( )

%*****************************************************************************80
%
%% companion_laguerre_test() tests companion_laguerre().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_laguerre_test():\n' );
  fprintf ( 1, '  companion_laguerre() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial p(x) in the Laguerre basis.\n' );

  p = [ 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 ];

  polynomial_laguerre_print ( p, '  Laguerre p(x)' );
%
%  Compute monomial form so we can get roots directly.
%
  n = length ( p );
  A = laguerre_to_monomial_matrix ( n );
  q = A * p';
  polynomial_monomial_print ( q, '  Monomial q(x)' );
  r1 = roots ( q(end:-1:1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Roots of q(x):\n' );
  disp ( r1 );
 
  A = companion_laguerre ( p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Laguerre companion matrix A(p):\n' );

  disp ( A );

  r = eigs ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Eigenvalues of A(p):\n' );
  disp ( r );

  return
end
function companion_legendre_test ( )

%*****************************************************************************80
%
%% companion_legendre_test() tests companion_legendre().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    30 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_legendre_test():\n' );
  fprintf ( 1, '  companion_legendre() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial p(x) in the Legendre basis.\n' );

  p = [ 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 ];

  polynomial_legendre_print ( p, '  Legendre p(x)' );
%
%  Compute monomial form so we can get roots directly.
%
  n = length ( p );
  A = legendre_to_monomial_matrix ( n );
  q = A * p';
  polynomial_monomial_print ( q, '  Monomial q(x)' );
  r1 = roots ( q(end:-1:1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Roots of q(x):\n' );
  disp ( r1 );
 
  A = companion_legendre ( p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Legendre companion matrix A(p):\n' );

  disp ( A );

  r = eigs ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Eigenvalues of A(p):\n' );
  disp ( r );

  return
end
function companion_monomial_test ( )

%*****************************************************************************80
%
%% companion_monomial_test() tests companion_monomial().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 March 2024
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'companion_monomial_test():\n' );
  fprintf ( 1, '  companion_monomial() computes the companion matrix\n' );
  fprintf ( 1, '  of a polynomial p(x) in the monomial basis.\n' );

  p = [ 6.0, 5.0, 4.0, 3.0, 2.0, 1.0 ];

  polynomial_monomial_print ( p, '  p(x)' );
%
%  Get roots directly.
%
  r1 = roots ( p(end:-1:1) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Roots of p(x):\n' );
  disp ( r1 );

  A = companion_monomial ( p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Monomial companion matrix A(p):\n' );

  disp ( A );

  r = eigs ( A );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Eigenvalues of A(p):\n' );
  disp ( r );

  return
end
function A = chebyshev_to_monomial_matrix ( n )

%*****************************************************************************80
%
%% chebyshev_to_monomial_matrix(): convert Chebyshev polynomial to monomial form.
%
%  Discussion:
%
%    1     0     -1      0       1     0     -1    0
%    0     1      0     -3       0     5      0   -7
%    0     0      2      0      -8     0     18    0
%    0     0      0      4       0   -20      0   56
%    0     0      0      0       8     0    -48    0
%    0     0      0      0       0    16      0 -112
%    0     0      0      0       0     0     32    0
%    0     0      0      0       0     0      0   64
%
%  Recursion:
%
%    T(0,X) = 1,
%    T(1,X) = X,
%    T(N,X) = 2 * X * T(N-1,X) - T(N-2,X)
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    31 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Input:
%
%    integer N, the order of the matrix.
%
%  Output:
%
%    real A(n,n), the matrix.
%
  A = zeros ( n, n );

  A(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  A(2,2) = 1.0;
 
  for i = 2 : n - 1
    A(1,i+1)     =                  - A(1,i-1);
    A(2:i-1,i+1) = 2.0 * A(1:i-2,i) - A(2:i-1,i-1);
    A(i,i+1) =     2.0 * A(i-1,i);
    A(i+1,i+1) =   2.0 * A(i,i);
  end
 
  return
end
function A = gegenbauer_to_monomial_matrix ( n, alpha )

%*****************************************************************************80
%
%% gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 April 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the order of A.
%
%    real alpha: the parameter.
%
%  Output:
%
%    real A(N,N): the matrix.
%
  A = zeros ( n, n );

  A(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  A(2,2) = 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 = 3 : n
    nn = j - 2; 
    c1 = ( 2 * nn + 2 * alpha     ) / ( nn + 1 );
    c2 = (   - nn - 2 * alpha + 1 ) / ( 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

  return
end
function A = h_to_monomial_matrix ( n )

%*****************************************************************************80
%
%% h_to_monomial_matrix(): physicist's Hermite to monomial conversion matrix.
%
%  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:
%
%    31 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the order of A.
%
%  Output:
%
%    real A(N,N), the matrix.
%
  A = zeros ( n, n );

  A(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  A(2,2) = 2.0;

  if ( n == 2 )
    return
  end

  for j = 3 : n
    for i = 1 : n
      if ( i == 1 )
        A(i,j) =                  - 2.0 * ( j - 2 ) * A(i,j-2);
      else
        A(i,j) = 2.0 * A(i-1,j-1) - 2.0 * ( j - 2 ) * A(i,j-2);
      end
    end
  end

  return
end
function A = laguerre_to_monomial_matrix ( n )

%*****************************************************************************80
%
%% 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:
%
%    integer N: the order of the matrix.
%
%  Output:
%
%    real A(N,N): the matrix.
%
  A = zeros ( n, n );

  A(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  A(1,2) = 1.0;
  A(2,2) = -1.0;

  for j = 3 : n
    for i = 1 : n
      if ( i == 1 )
        A(i,j) = (  ( 2 * j - 3 ) * A(i,j-1) ...
                 +  (   - j + 2 ) * A(i,j-2) ) ...
                 /  (     j - 1 );
      else
        A(i,j) = (  ( 2 * j - 3 ) * A(i,j-1) ...
                 -  (         1 ) * A(i-1,j-1) ...
                 +  (   - j + 2 ) * A(i,j-2) ) ...
                 /  (     j - 1 );
      end
    end
  end

  return
end
function A = legendre_to_monomial_matrix ( n )

%*****************************************************************************80
%
%% legendre_to_monomial_matrix(): Legendre coefficient conversion matrix.
%
%  Discussion:
%
%    If PL(x) is a linear combination of Legendre polynomials
%    with coefficients CL, then PM(x) is a linear combination of
%    monomials with coefficients CM = A * CL.
%    
%    The coefficients are ordered so the constant term is first.
%
%  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:
%
%    21 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N, the order of A.
%
%  Output:
%
%    real A(N,N), the matrix.
%
  A = zeros ( n, n );

  A(1,1) = 1.0;

  if ( n == 1 )
    return
  end

  A(2,2) = 1.0;

  for j = 3 : n
    for i = 1 : n
      if ( i == 1 )
        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
    end
  end

  return
end
function polynomial_chebyshev_print ( c, label )

%*****************************************************************************80
%
%% polynomial_chebyshev_print() prints a polynomial in the Chebyshev basis.
%
%  Discussion:
%
%    The form of a Chebyshev polynomial is:
%
%      p(x) = c(0)*T0(x) + c(1)*T(1)(x) + ... + c(n)*T(n)(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    29 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      fprintf ( 1, '%+14f * T%d(x)\n', c(i+1), i );

    end

  end

  return
end
function polynomial_gegenbauer_print ( c, label )

%*****************************************************************************80
%
%% polynomial_gegenbauer_print() prints a polynomial in the Gegenbauer basis.
%
%  Discussion:
%
%    The form of a Ggenbauer polynomial is:
%
%      p(x) = c(0)*C0(x) + c(1)*C1(x) + ... + c(n)*C(n)(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    30 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      fprintf ( 1, '%+14f * C%d(x)\n', c(i+1), i );

    end

  end

  return
end
function polynomial_hermite_print ( c, label )

%*****************************************************************************80
%
%% polynomial_hermite_print() prints a polynomial in the Hermite basis.
%
%  Discussion:
%
%    The form of a Hermite polynomial is:
%
%      p(x) = c(0)*H0(x) + c(1)*H(1)(x) + ... + c(n)*H(n)(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    29 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      fprintf ( 1, '%+14f * H%d(x)\n', c(i+1), i );

    end

  end

  return
end
function polynomial_laguerre_print ( c, label )

%*****************************************************************************80
%
%% polynomial_laguerre_print() prints a polynomial in the Laguerre basis.
%
%  Discussion:
%
%    The form of a Laguerre polynomial is:
%
%      p(x) = c(0)*L0(x) + c(1)*L(1)(x) + ... + c(n)*L(n)(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    31 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial coefficients.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      fprintf ( 1, '%+14f * L%d(x)\n', c(i+1), i );

    end

  end

  return
end
function polynomial_legendre_print ( c, label )

%*****************************************************************************80
%
%% polynomial_legendre_print() prints a polynomial in the Legendre basis.
%
%  Discussion:
%
%    The form of a Legendre polynomial is:
%
%      p(x) = c(0)*P0(x) + c(1)*P(1)(x) + ... + c(n)*P(n)(x)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    30 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial coefficients.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      fprintf ( 1, '%+14f * P%d(x)\n', c(i+1), i );

    end

  end

  return
end
function polynomial_monomial_print ( c, label )

%*****************************************************************************80
%
%% polynomial_monomial_print() prints a polynomial in monomial basis.
%
%  Discussion:
%
%    The monomial form of a polynomial is:
%
%      p(x) = c(0)*x^0 + c(1)*x + ... + c(n)*x^(n)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    23 October 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real C(1:D+1): the polynomial.
%
%    character LABEL: an optional title.
%
  d = length ( c ) - 1;

  fprintf ( 1, '\n' );
  if ( 0 < length ( label ) )
    fprintf ( 1, '%s = \n', label );
  end

  if ( d < 0 )
    fprintf ( 1, '  Zero polynomial\n' );
    return
  end

  if ( all ( c == 0 ) )
    fprintf ( 1, '  0\n' );
    return
  end

  for i = d : -1 : 0

    if ( c(i+1) ~= 0 )

      if ( 2 <= i )
        fprintf ( 1, '%+14f * x^%d\n', c(i+1), i );
      elseif ( i == 1 )
        fprintf ( 1, '%+14f * x\n', c(i+1) );
      elseif ( i == 0 )
        fprintf ( 1, '%+14f\n', c(i+1) );
      end

    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
