/* characteristic polynomial - Dorny, P4.5, page 224 

   returns coefficients b of the polynomial: x^n + b1 x^(n-1) + ... + bn

   b1 = -trace(A)

   b2 = -(1/2)(b1 trace(A) + trace(A^2))
   ...

   bn = -(1/n)(bn-1 trace(A) + ... + b1 trace(A^n-1) + trace(A^n))

*/

vector charpoly( matrix A)
{
  int n = rows(A);  int i, c;  vector t = zeros(n,1), b = t;

  matrix B = eye(n,n);

  for( c = 1; c <= n; ++c)
  {
    B = A*B; /* A^c */

    t(c) = sum(diag(B)); /* trace(A^c) */

    b(c) = t(c);

    for( i = 1; i < c; ++i)
    {
      b(c) = b(c) + b(i)*t(c-i);
    }

    b(c) = -b(c)/c;
  }

  return b;
}
