function isStable = idRecursiveEstimationCheckBistritzCondition(p)
% idRecursiveEstimationCheckBistritzCondition Test if all roots of the
% polynomial with coefficients p are in the unit disk. p must be a
% row vector of real numbers. p(1) must be positive.
%
%   isStable = idRecursiveEstimationCheckBistritzCondition(p)
%
%   Algorithm: Yuval Bistritz, Zero Location With Respect To The Unit
%   Circle Of Discrete-Time Linear System Polynomials, 1984 IEEE

%   Copyright 2013-2015 The MathWorks, Inc.

%#codegen

%% Check simple, necessary conditions
if isscalar(p)
    % Static gain
    isStable = true();
    return;
end

% Necessary conditions with the p(1)>1 assumption (pg. 1134)
if p(1)<=abs(p(end)) || sum(p,2)<=0
    isStable=false();
    return;
end

%% Initizalization

% pr: reciprocal polynomial of p
pr = p(end:-1:1);

% Eq. (12a)
T_n = p + pr; 
if T_n(1)==0 % sum(p,2)<=0 implies sum(T_n,2)>0 necessary&sufficient condition
    isStable = false();
    return;
end

% T_nm1=deconv(Tp,[1 -1]) or Tp/(z-1) (Eq. (12b))
nt = uint16(numel(p));
T_nm1 = zeros(1,nt-1,'like',p);
T_nm1(1) = p(1)-pr(1);
for kk=2:nt-1
    T_nm1(kk) = p(kk)-pr(kk)+T_nm1(kk-1);
end

% Stability:
% Necessary condition Eq. (14)
% Necessary and sufficient condition Eq. (31)
if T_nm1(1)==0 || sum(T_nm1,2)<=0; 
    isStable = false();
    return;
end

%% Recursion step
for kk=nt-1:-1:2
    delta = T_n(1)/T_nm1(1); % Eq. (25)
    % delta>0 is a necessary stability condition but skip cheking it
    for mm=1:kk-1 
        T_n(mm) = T_nm1(mm);
        T_nm1(mm) = delta*(T_nm1(mm)+T_nm1(mm+1))-T_n(mm+1); % Eq. (26)
    end
    T_n(kk) = T_nm1(kk);
    % Stability:
    % Necessary condition Eq. (14)
    % Necessary and sufficient condition Eq. (31)
    sigmaI = T_nm1(1);
    if kk>=2
        for mm=2:kk-1
           sigmaI = sigmaI +  T_nm1(mm);
        end
    end
    if T_nm1(1)==0 || sigmaI<=0; % sign change
        isStable = false();
        return;
    end
end
isStable = true();
end