function roots=BairstowRoots(pln,epsilon)
% BaristowRoots - Nalazenje nula polinoma Berstovljevom metodom
%
% Argumenti funkcije su :
%   pln - vektor vrste koji sadzi koeficijente polinoma
%   epsilon - zeljena tacnost, ako se ne navede podrazumeva se MATLABOV eps
%
% Funkcija vraca:
%   roots - vektor nula polinoma

% Autor : Filip Maric
% Verzija : 1.0
% Datum : 30.3.2003.

if nargin<2
    epsilon=eps;
end


% Izracunavamo dimenzije vektora pln
[m,n]=size(pln);

% Proveravamo da li se radi o vektoru
if m>1
    error('Polinom se zadaje vektorom!')
end

% Eliminisemo vodece nule polinoma
nenula = find(pln);
if isempty(nenula),
    % Svi elementi su nule
    return
end
pln = pln(nenula(1):nenula(end));


% Polinom prevodimo u oblik x^m+a1x^(m-1)+...+am
if (pln(1)~=1)
    pln=pln./pln(1);
end

% Racunamo vektor koeficijenata ai
a=pln(2:end);
m=length(a);

% Promenjiva u kojoj cemo smestati nadjene nule polinoma
roots=[];

% Dok god je nas polinom stepena >2 pokusavamo da izdvojimo kvadratni faktor
% tj. x^m+a1x^(m-1)+...+am=(x^2+px+q)(x^(m-2)+b1x^(m-3)+...+b_(m-2))

while m>2
    % Pretpostavljamo da je p=q=0 i onda iterativno, Newton-Rapshonovom metodom
    % vrsimo popravke
    p=0;
    q=0;
    while(1)
        
        % Izracunavamo vektor koeficijenata b
        b(1)=a(1)-p;
        b(2)=a(2)-p*b(1)-q;
        for k=3:m
            b(k)=a(k)-p*b(k-1)-q*b(k-2);
        end

        % Izracunavamo izvode koeficijenata bk po p i q
        % Njih izracunavamo po rekurentnoj formuli :
        %     ck=bk-pc_(k-1)-qc_(k-2), k=1,...,m-1  c_(-1)=0, c_1=0
        % Primetimo, da posto su nam potrebni samo izvodi c_(m-2) i c_(m-3)
        % Nema potrebe da pamtimo ceo niz c
        
        cpp=1;
        cp=b(1)-p;
        for k=2:m-2
            tmp=cp;
            cp=b(k)-p*cp-q*cpp;
            cpp=tmp;
        end
        
        % Resavamo sistem jednacina po popravkama
         c_prim=-(p*cp+q*cpp);
         D=cp^2-c_prim*cpp;
         Dp=b(m-1)*cp-b(m)*cpp;
         Dq=-b(m-1)*c_prim+b(m)*cp;

         deltap=Dp/D;
         deltaq=Dq/D;
        
       
        
        % Popravke dodajemo na pretpostavljano resenje
        p=p+deltap;
        q=q+deltaq;
    
        % Ukoliko je postignuta zeljena tacnost zavrsavamo sa primenom
        % Newton-Rapshsonove metode
        if (abs(deltap)<epsilon & abs(deltaq)<epsilon)
            break;
        end
    end
    
    % Nalazimo korene polinoma x^2+px+q
    discr=sqrt(p^2-4*q);
    x0=(-p+discr)/2;
    x1=(-p-discr)/2;
    
    % Dodajemo korene vektoru roots
    roots=[x0; x1; roots];
    
    % Prelazimo na trazenje nula nadjenog kolicnika
    %    x^(m-2)+b1x^(m-3)+...+b_(m-2)
    m=m-2;

    % Izracunavamo kolicnik polinoma a i (x^2+px+q)
    a(1)=a(1)-p;
    a(2)=a(2)-p*a(1)-q;
    for k=3:m
      a(k)=a(k)-p*a(k-1)-q*a(k-2);
    end
    
end

% U ovom slucaju nam je ostao polinom stepena <=2
% x^m+a1x^(m-1)+a2x^(m-2)
if m==2
    % Ukoliko se radi o kvadratnom polinomu, nalazimo njegove nule
    discr=sqrt(a(1)^2-4*a(2));
    x0=(-a(1)+discr)/2;
    x1=(-a(1)-discr)/2;
    roots=[x0; x1; roots];
else 
    if m==1 
        %Ukoliko se radi o linarnom polinomu nalazimo njegovu jedinu nulu
        x0=-a(1);
        roots=[x0; roots];
    else
        %Ukoliko je polinom 0 stepena svodi se na 1=0
        warning('Polinom nema nula')
    end
end