//-------------------------------------------------------------------//
//  roots

//  Syntax: roots(P)

//  Description:

//  Find roots of polynomial P of degree n.
//
//  P is a vector (row or column) with n+1 components respresenting
//  polyniomial coefficients in decending powers.
// 
//  Return a column vetor containing n roots of P.

//  Example:
//  the following script finds the roots of 
//  x^3 + 3*x^2 + 5 = 0
// 
//  > p=[1,3,0,5]
//   p =
//          1          3          0          5
//  > r=roots(p)
//   r =
//                -3.43 + 0i
//             0.213 - 1.19i
//             0.213 + 1.19i

//  See Also: poly

//  Tzong-Shuoh Yang (tsyang@ce.berkeley.edu)  5/7/94
//
//-------------------------------------------------------------------//
roots = function(P)
{
   local(c,p,r,z,nz);

   if (P.n <= 1) {
      error("Argument of roots() is a scalar, not a polynomial.");
   else if (min(size(P)) != 1) {
      error("Argument of roots() must be a column or a row vector.");
   }}
   p = P[:]';
   // trim leading and trailing zeros
   z  = find(p != 0);
   nz = p.n - z[z.n];
   p  = p[z[1]:z[z.n]];
   // find companion matrix
   c = compan(p);
   // find roots and stuff with trailing zeros
   r = [zeros(nz,1);eign(c).val'];  
   // trim complex part if all real
   if (all(imag(r) == 0)) {
	  r = real(r);
   }
   return r;
};
