!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  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.                                                     !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  compiled with gfortran -fdefault-real-8 Bairstow_Root_Finding.f95  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  ABW 2012 - 2015  !
!!!!!!!!!!!!!!!!!!!!!
program bairstow
  implicit none
  
  integer,parameter :: ns = 6  !  Order of the starting polynomial
  integer,parameter :: itmax = 200  !  Max number of iterations
  real(kind=8),parameter :: tol = 1d-9  !  Tolerance for zero remainder
  real(kind=8) :: u,v,detJ,s1,s2,du,dv
  real(kind=8),allocatable :: a(:),b(:),c(:)
  integer :: j,k,n
  
  allocate(a(0:ns), b(0:ns), c(0:ns))
  
  !  Initialize the starting polynomial:
  !  a(n)*x^n + a(n-1)*x^(n-1) +...+ a(1)*x + a(0)
  a(0) = -7.0
  a(1) = 6.0
  a(2) = 5.0
  a(3) = 4.0
  a(4) = 3.0
  a(5) = 2.0
  a(6) = 1.0
  
  write(*,*) 'Input polynomial:'
  write(*,101) a(6),'*x^6 +',a(5),'*x^5 +',a(4),'*x^4 +',a(3),&
       '*x^3 +',a(2),'*x^2 +',a(1),'*x -', -a(0)
  101 format (F4.1,A6,F4.1,A6,F4.1,A6,F4.1,A6,F4.1,A6,F4.1,A4,F4.1)
  
  !  Initialize the 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.
  do n = ns,2,-2
     
     b(n) = a(n)
     c(n) = 0.0
     c(n-1) = a(n)
     
     write(*,*)
     write(*,"(A26,I2)") 'Input polynomial of order',n
     write(*,*)
     write(*,"(A4,A14,A20)") 'it#','u','v'
     write(*,*)
     
     do j = 1,itmax
	!  Recursively determine the coefficients of the 
	!  new polynomial.  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-1) = a(n-1) + u*b(n)         
        do k = n-2,0,-1
           b(k) = a(k) + u*b(k+1) + v*b(k+2)
           c(k) = b(k+1) + u*c(k+1) + v*c(k+2)
        enddo
        
        detJ = c(0)*c(2) - c(1)*c(1)
        if (detJ == 0.0) then
           write(*,*) 'ERROR:  detJ = 0, division by zero'
           exit
        endif
        du = (c(1)*b(1) - c(2)*b(0))/detJ
        dv = (c(1)*b(0) - c(0)*b(1))/detJ
        u = u + du
        v = v + dv
        
        write(*,"(I3,F20.9,F20.9)") j,u,v
        
        if (max(abs(du),abs(dv)) < tol) then   
           exit
        elseif (max(abs(du),abs(dv)) > tol .and. j == itmax) then
           write(*,*)
           write(*,*) 'ERROR:  Tolerance not met, increase itmax or adjust (u,v).'
           exit
        endif
     enddo
     
     if ((u*u + 4.0*v) < 0.0 .and. n > 1) then
        write(*,*)
        write(*,*) 'Solution is complex with roots x1 and x2:'
        s1 = u/2.0
        s2 = (sqrt(abs(u*u + 4.0*v)))/2.0
	write(*,"(A6,F14.9,A7,F14.9,A2)") 'x1 = ',s1,' + i *(',s2,')'
	write(*,"(A6,F14.9,A7,F14.9,A2)") 'x2 = ',s1,' - i *(',s2,')'
     elseif ((u*u + 4.0*v) > 0.0 .and. n > 1) then
        write(*,*)
        write(*,*) 'Solution is real with roots x1 and x2:'
        s1 = (u + sqrt(u*u + 4.0*v))/2.0
        s2 = (u - sqrt(u*u + 4.0*v))/2.0
        write(*,"(A5,F14.9)") 'x1 =',s1
        write(*,"(A5,F14.9)") 'x2 =',s2
     endif
     
     !  If the order of the input polynomial is 3, then
     !  the order or the polynomial for the next step is
     !  1 (linear) and the solution is already known.
     if (n == 3) then  
        write(*,*)
        write(*,*) 'Input polynomial of order',n-2
        write(*,*)
        write(*,*) 'Solution is real with root x:'
        write(*,"(A4,F14.9)") 'x =',-b(2)/b(3)
        exit
     endif
     deallocate(a)
     allocate(a(0:ns-2))
     
     !  Reinitialize the input polynomial with order n-2
     do k = n,2,-1
        a(k-2) = b(k)
     enddo
     deallocate(b,c)
     allocate(b(0:ns-2),c(0:ns-2))
  enddo
  deallocate(a,b,c)
  
end program bairstow
