/*  -*- MaTX -*-
 *
 *  NAME
 *      deconv() - Deconvolution
 *
 *  SYNOPSIS
 *      {q,r} = deconv(b,a)
 *         Matrix q,r;
 *         Matrix b,a;
 *
 *  DESCRIPTION
 *      deconv(b,a) produces the deconvolution of b by a vector a. 
 *
 *      The result is returned in q and the remainder in r such that 
 *
 *        b = conv(q,a) + r
 *
 *  SEE ALSO
 *      conv
 */

Func List deconv(b, a)
	Matrix b, a;
{
	Matrix q, r;
	Integer na, nb;

	nb = length(b);
	na = length(a);

	if (na > nb) {
		q = Z(1);
		r = b;
		return {q,r};
	}
	if (a(1) == 0.0) {
		error("deconv(): First coefficient must be nonzero.\n");
	}

	{q} = filter(b, a, [[1], Z(1,nb-na)]);

	if (Rows(b) != 1) {
		q = makecolv(q);
	}

	r = b - Matrix(conv(q, a));
	return {q,r};
}