function [v,u,y,x]=gpfx(b,a,test)
%  [v,u,y,x]=gpfx(b,a,test)
% gross partial fraction expansion: V/U + Y/X = B/A,
% UX=A, where the roots of U satisfy
% the conditions in the string 'test'.
% examples: test='abs(z)<1', test='real(z)<0'
% Polynomials do NOT use the MATLAB convention!
% rather, [b1 b2 b2] represents b1+b2*z^-1 +b2*z^-2
% Thus conv(u,x)=a, and conv(v,x)+conv(y,u)=b 
% (each term extended to the right). The leading
% coefficients of u and x are unity.
% The leading coefficient of y is zero.

	na=length(a);nb=length(b);u=1;x=1;r=roots(a);
	for k=[1:na-1]
		z=r(k);
		if eval(test)
			u=conv(u,[1 -z]);
		else
			x=conv(x,[1 -z]);
		end
	end
	u=real(u);x=real(x);
	mu=length(u);mx=length(x);
	
	if na >= nb
		bb=[b,zeros(1,na-nb)];
		uu=u;
		M=zeros(na,na);
	else
		bb=b;
		mu=nb+1-mx;
		uu=[u,zeros(1,mu-length(u))];
		M=zeros(nb,nb);
	end		
	for k=1:mu
		M(k,[k:k+mx-1])=x;
	end
	for k=2:mx
		M(k+mu-1,[k:k+mu-1])=uu;
	end
	s=bb*inv(M);
	v=s(1:mu);
	y=[0,s(mu+1:mu+mx-1)];