# v_rootstab

## PURPOSE

V_ROOTSTAB determines the number of polynomial roots outside, inside and on the unit circle [NO,NI,NC]=v_rootstab(P)

## SYNOPSIS

function [no,ni,nc]=v_rootstab(p)

## DESCRIPTION

    V_ROOTSTAB determines the  number of polynomial roots outside, inside and on the unit circle [NO,NI,NC]=v_rootstab(P)

      Inputs:  p   Polynomial with real or complex coefficients

     Outputs: no   Number of roots outside the unit circle
              ni   Number of roots inside the unit circle
              nc   Number of roots lying on the unit circle

     This routine uses the algorithm given in [1]. Rounding errors may cause roots that lie on the unit circle to migrate to either inside or outside.

     Refs:
       [1] Messaoud Benidir. On the root distribution of general polynomials with respect to the unit circle.
           Signal Processing, 53 (1): 75–82, August 1996. ISSN 0165-1684. doi: 10.1016/0165-1684(96)00077-1.
           [ Note error in equation (52) which should read k^ = +k - Q(0)/(conj(c0)*(1/c - c)) ]

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 function [no,ni,nc]=v_rootstab(p)
    0002 %V_ROOTSTAB determines the  number of polynomial roots outside, inside and on the unit circle [NO,NI,NC]=v_rootstab(P)
    0003 %
    0004 %  Inputs:  p   Polynomial with real or complex coefficients
    0005 %
    0006 % Outputs: no   Number of roots outside the unit circle
    0007 %          ni   Number of roots inside the unit circle
    0008 %          nc   Number of roots lying on the unit circle
    0009 %
    0010 % This routine uses the algorithm given in [1]. Rounding errors may cause roots that lie on the unit circle to migrate to either inside or outside.
    0011 %
    0012 % Refs:
    0013 %   [1] Messaoud Benidir. On the root distribution of general polynomials with respect to the unit circle.
    0014 %       Signal Processing, 53 (1): 75–82, August 1996. ISSN 0165-1684. doi: 10.1016/0165-1684(96)00077-1.
    0015 %       [ Note error in equation (52) which should read k^ = +k - Q(0)/(conj(c0)*(1/c - c)) ]
    0016
    0017 %      Copyright (C) Mike Brookes 2025
    0018 %
    0019 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0020 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0021 %
    0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0023 %   This program is free software; you can redistribute it and/or modify
    0024 %   it under the terms of the GNU Lesser General Public License as published by
    0025 %   the Free Software Foundation; either version 3 of the License, or
    0026 %   (at your option) any later version.
    0027 %
    0028 %   This program is distributed in the hope that it will be useful,
    0029 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0030 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0031 %   GNU Lesser General Public License for more details.
    0032 %
    0033 %   You can obtain a copy of the GNU Lesser General Public License from
    0034 %   https://www.gnu.org/licenses/ .
    0035 %
    0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0037
    0038 no=0;                                       % initialize count of roots outside unit circle
    0039 nc=0;                                       % initialize count of roots on the unit circle
    0040 if all(p==0)                                % p is all zero
    0041     ni=0;
    0042 else
    0043     p=p(:).';                               % force to be a row vector
    0044     p=p(find(p~=0,1):end);                  % trim leading zeros
    0045     np=length(p);                           % initial 1+order(p)
    0046     np0=np;                                 % save initial order for calculating ni later
    0047     npd=0;                                  % initialize saved order for calculating nc
    0048     while length(p)>1
    0049         p=p/sqrt(p*p');                     % normalize p each loop
    0050         pf=conj(p(np:-1:1));                % flipped version of p
    0051         k=-p(np)/pf(np);
    0052         q=p+k*pf;                           % null out the constant coefficient
    0053         q(np)=0;                            % force exact zero in case of rounding errors
    0054         if all(q==0)
    0055             p=p(1:end-1).*(np-1:-1:1);      % take derivative
    0056             if npd==0
    0057                 npd=np;                     % save current 1+order(p)
    0058                 nod=no;                     % save current count of roots outside unit circle
    0059             end
    0060         elseif q(1)==0                      % if |k|=1 and q~=0
    0061             q=q(1:find(q~=0,1,'last'));     % trim trailing zeros from q
    0062             dr=-q(end)/(pf(end)*k);         % direction we will move normalized k in
    0063             if abs(real(dr))>abs(imag(dr))
    0064                 cf=abs(real(dr));           % increase or decrease real part by 0.5
    0065             else
    0066                 cf=0.25*abs(imag(dr));      % increase or decrease imag part by 2
    0067             end
    0068             c=sqrt(1+cf^2)-cf;             % choose c in range (0,1)
    0069             p=p/c+k*c*pf+[zeros(1,np-length(q)) q];
    0070         elseif abs(k)>1
    0071             q=q(1:find(q~=0,1,'last'));     % trim trailing zeros from q
    0072             p=conj(q(end:-1:1));            % flip q
    0073             no=no+np-length(p);             % increement count of zeros outside unit circle
    0074         else
    0075             p=q(1:find(q~=0,1,'last'));     % trim trailing zeros from q and take as p
    0076         end
    0077         np=length(p);                       % 1+order(p)
    0078     end
    0079     if npd>0
    0080         nc=npd-1-2*(no-nod);
    0081     end
    0082     ni=np0-1-nc-no;
    0083 end