/* -*- MaTX -*-
 *
 *  NAME
 *      resolvent() - Resolvent matrix (s*I - A)^(-1) of A
 *
 *  SYNOPSIS
 *      {N,ch} = resolvent(A)
 *         PoMatrix N;
 *         Polynomial ch;
 *         Matrix A;
 *  
 *  DESCRIPTION
 *      resolvent(A) returns the numerator N and denominator ch, which
 *      is the characteristic polynomial, of the resolvent matrix
 *
 *           N(s)           -1
 *          ----- = (s*I - A)
 *          ch(s)
 *
 *      of matrix A.
 *
 *  SEE ALSO
 *      faddeev, makepoly, and charpoly
 */

Func List resolvent(A)
	Matrix	A;
{
	Integer	i, n;
	Polynomial s,ch;
	Matrix	Gamma,AA;
	PoMatrix N;

	n = Rows(A);
	AA = (A,:);

	{Gamma, ch} = faddeev(A);

	s = Polynomial("s");
	N = Z(n, n, s);
	if (version() == 4) {
		for (i = 0; i < n; i++) {
			N = N + Gamma(0, i, AA) * s^i;
		}
	} else {
		for (i = 0; i < n; i++) {
			N = N + Gamma(1, i+1, AA) * s^i;
		}
	}

	return {N,ch};
}