%  Bairstow complex root finding method for polynomials of order n > 1
%=================================================================%
%  Consider a polynomial p(x) = a_n*x^n + a_n-1*x^n-1 +...+ a_0   % 
%  of aribtrary order, n.  A polynomial x^2 - u*x - v of order,   % 
%  2, is constructed using Newton-Raphson root finding that       % 
%  shares a pair of roots with the original polynomial.  p(x)     % 
%  is then divided by the constructed polynomial and produces a   % 
%  polynomial q(x) = b_n*x^n-2 + b_n-1*x^n-3 +...+ b_2 of order,  % 
%  n-2, and the Bairstow method is iterated until the remainder   % 
%  of the division is approximately zero.  q(x) is then used      % 
%  as the input polynomial.  This process is iterated until the   % 
%  order of the input polynomial is 1 or 2 and all roots of p(x)  % 
%  are found.                                                     %
%=================================================================%
%  ABW 2012 - 2015  %
%===================%
clc;

%  Iniitialize parameters
ns = 7;  %  Order of starting polynomial
it_max = 200; %  Maximum number of iterations
tol = 1.0e-9;  %  Tolerance for zero remainder

%  Initialize coefficient arrays
b = zeros(1,ns+1);
c = zeros(1,ns+1);

%  Initialize starting polynomial coefficients
%  a(n)*x^n + a(n-1)*x^(n-1) + ... + a(1)*x + a(0)
a = [-7.0,6.0,5.0,4.0,3.0,2.0,1.0,-1.0];

%  Display input polynomial
disp('Input polynomial:')
fprintf('%dx^7 + %dx^6 + %dx^5 + %dx^4 + %dx^3 + %dx^2 + %dx + %d \n \n'...
    ,a(8),a(7),a(6),a(5),a(4),a(3),a(2),a(1));

%  Initialize constructed polynomial
%  x^2 - u*x - v
u = 3.0;
v = -4.0;

%  Iteration starts with the order of the input polynomial and each
%  step is decreased by the order of the constructed polynomial, 2.
for n = ns:-2:2 
    
    b(n+1) = a(n+1);
    c(n+1) = 0.0;
    c(n) = a(n+1);
    
    fprintf('\nInput polynomial of order, %d \n \n',n)
    disp('it#           u                v')
    
    for i = 1:it_max
        %  Recursively determine the coefficients of the new polynomial.  
        %  The new polynomial should share two roots with the input 
        %  polynomial as long as the remainder (du,dv) is approximately 
        %  zero (du,dv < tol).
        b(n) = a(n) + u*b(n+1);
        for j = n-2:-1:0
            b(j+1) = a(j+1) + u*b(j+2) + v*b(j+3);
            c(j+1) = b(j+2) + u*c(j+2) + v*c(j+3);
        end
        detJ = c(1)*c(3) - c(2)*c(2);
        if detJ == 0
            disp('ERROR: detJ = 0, division by zero')
            break
        end
        du = (c(2)*b(2) - c(3)*b(1))/detJ;
        dv = (c(2)*b(1) - c(1)*b(2))/detJ;
        u = u + du;
        v = v + dv;
        
        fprintf('%d     %14.9f   %14.9f \n',i,u,v)
        
        if max(abs(du),abs(dv)) < tol
            break
        elseif max(abs(du),abs(dv)) > tol && i == it_max
            disp('ERROR:  Tolerance not met')
            disp('Increase it_max or adjust (u,v).')
            break
        end
    end
    
    if u*u + 4.0*v < 0.0 && n > 1
        fprintf('\nSolution is complex with roots x1 and x2: \n')
        s1 = u/2.0;
        s2 = sqrt(abs(u*u + 4.0*v))/2.0;
        fprintf('x1 = %14.9f + i*%14.9f \n',s1,s2)
        fprintf('x1 = %14.9f - i*%14.9f \n',s1,s2)
        fprintf('\n=========================================\n')
    elseif u*u + 4.0*v > 0.0 && n > 1
        fprintf('\nSolution is real with roots x1 and x2: \n')
        s1 = (u + sqrt(u*u + 4.0*v))/2.0;
        s2 = (u - sqrt(u*u + 4.0*v))/2.0;
        fprintf('x1 = %14.9f, x2 = %14.9f \n',s1,s2)
        fprintf('\n=========================================\n')
    end
    
    %  If the order of the input polynomial is 3, then the order of
    %  the polynomial for the next step is 1 (linear) and the
    %  solution is already known.
    if n == 3
        fprintf('\nInput polynomial of order, %d \n \n',n-2)
        fprintf('Solution is real with root x:\n')
        s3 = -b(3)/b(4);
        fprintf('x = %14.9f \n',s3)
        fprintf('\n=========================================\n')
        break
    end
    
    %  Reinitialize the input polynomial with order n-2
    a = zeros(1,n-1);
    for k = n:-1:2
        a(k-1) = b(k+1);
    end
    b = zeros(1,n-1);
    c = zeros(1,n-1);
end