function [stable] = stabilitycheck(A);

N = length(A)-1; % Order of A(z)
stable = 1;      % stable unless shown otherwise
A = A(:);        % make sure it's a column vector
for i=N:-1:1
  rci=A(i+1);
  if abs(rci) >= 1
    stable=0;
    return;
  end
  A = (A(1:i) - rci * A(i+1:-1:2))/(1-rci^2);
% disp(sprintf('A[%d]=',i)); A(1:i)'
end

% TSC - test function stabilitycheck, comparing against
%       pole radius computation

N = 200; % polynomial order
M = 20; % number of random polynomials to generate
disp('Random polynomial test');
nunstp = 0; % count of unstable A polynomials
sctime = 0; % total time in stabilitycheck()
rftime = 0; % total time computing pole radii
for pol=1:M
  A = [1; rand(N,1)]'; % random polynomial
  tic; 
  stable = stabilitycheck(A); 
  et=toc;  % Typ. 0.02 sec Octave/Linux, 2.8GHz Pentium
  sctime = sctime + et;
  % Now do it the old fashioned way
  tic; 
  poles = roots(A); % system poles
  pr = abs(poles);  % pole radii
  unst = (pr >= 1); % bit vector
  nunst = sum(unst);% number of unstable poles
  et=toc; % Typ. 2.9 sec Octave/Linux, 2.8GHz Pentium
  rftime = rftime + et;
  if stable, nunstp = nunstp + 1; end
  if (stable & nunst>0) | (~stable & nunst==0)
    error('*** stabilitycheck() and poles DISAGREE ***'); 
  end
end
disp(sprintf(...
     ['Out of %d random polynomials of order %d,',...
      ' %d were unstable'], M,N,nunstp));