function c2 = r8poly_power ( d1, c1, p )

%*****************************************************************************80
%
%% r8poly_power() computes a positive integer power of a polynomial.
%
%  Discussion:
%
%    The power sum form of a polynomial is:
%
%      p(x) = c(0) + c(1) * x + ... + c(n-1) * x^(n-1) + c(n) * x^(n)
%
%    Because MATLAB doesn't allow 0 indexing, we use an adjusted formula:
%
%      p(x) = c(1) + c(2) * x + ... + c(n) * x^(n-1) + c(n+1) * x^(n)
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    21 October 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer D1: the degree of the polynomial.
%
%    real C1(1:D1+1): the polynomial to be raised to the power.
%
%    integer P: the nonnegative power to which the polynomial is raised.
%
%  Output:
%
%    real C2(P*D1+1): the power of the polynomial.
%

%
%  Zero out C2.
%
  d2 = p * d1;
  c2(1:d2+1) = 0.0;
%
%  Search for the first nonzero element.
%
  nonzer = -1;

  for i = 0 : d1
    if ( c1(i+1) ~= 0.0 )
      nonzer = i;
      break
    end
  end

  if ( nonzer == -1 )
    return;
  end

  c2(0+1) = c1(nonzer+1)^p;

  for i = 1 : p*(d1-nonzer)

    if ( i + nonzer <= d1 )
      c2(i+1) = i * p * c2(1) * c1(i+nonzer+1);
    else
      c2(i+1) = 0.0;
    end

    for j = 1 : i - 1

      if ( j + nonzer <= d1 )
        c2(i+1) = c2(i+1) - ( i - j ) * c1(j+nonzer+1) * c2(i-j+1);
      end

      if ( i - j + nonzer <= d1 )
        c2(i+1) = c2(i+1) + ( i - j ) * p * c2(j+1) * c1(i-j+nonzer+1);
      end

    end

    c2(i+1) = c2(i+1) / ( i * c1(nonzer+1) );

  end
%
%  Shift up.
%
  for i = p*nonzer : p*d1
    c2(i+1) = c2(i-p*nonzer+1);
  end

  for i = 0 : p * nonzer-1
    c2(i+1) = 0.0;
  end

  return
end
