      program rootf
*
* This is a polynomial root finder program.
* The polynomials must have real coefficients.
* The method applied is a modification of BAIRSTOW's 
* i.e. division by a quadratic polynomial x**2 +p1*x +p2
* determining p1 & p2 to be such that the remainder is vanishing.
* Then a deflation is applied by synthetic division to reduce
* the polynomial degree by two, and the procedure is repeated.
* A check is always made if the degree is reduced to either
* one or two to use special finishing procedure
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      common/toroot/a(0:100),b(0:100),kco
      COMPLEX ROOT(100),disc
      parameter (N=2,M=2)
      dimension icode(4),xp(2),xll(2),xrl(2),ixat(2)
      character line*78
      data icode/1,3*0/
      do 100 i=1,78
100   line(I:I)='='
* Read in, tolerance, polynomial degree, polynomial coefficients
      write(*,*) 'Enter a termination tolerance: '
      read(*,*) eps
      write(*,*) 'Enter the degree of the polynomial:'
      read(*,*)kco
      write(*,*) 'Enter the polynomial coefficients:( a(0),a(1),...)'
      read(*,*)(a(i),i=0,kco)
* Type the above input for reference
      write(*,'(A)')line
      write(*,*) '    RUNNING WITH THE FOLLOWING INPUT  '
      write(*,*) 'TOLERANCE: ',eps
      write(*,*) 'DEGREE:    ',kco
      write(*,*)'        COEFFICIENTS  a(i), i=0,1,... '
      write(*,*)(a(i),i=0,kco)
      write(*,'(A)')line
      nor = 0
30    continue
* Check if the degree is one, and if so solve a linear equation.
      if (kco.eq.1) then
        write(*,*)' The roots of the polynomial are:'
        nor = nor + 1
        root(nor)= cmplx(-a(0)/a(1),0)
        do 1 i=1,nor
           write(*,19)i,root(i)
1       continue
19      format(2x,'Root # ',i4,t17,'(',g14.7,',  ',g14.7,')')
        stop
      endif
* Check if the degree is two, and if so solve a quadratic equation.
      if (kco .eq. 2) then
        write(*,*)' The roots of the polynomial are:'
        nor = nor + 1
        disc = cmplx(a(1)**2 - 4*a(2)*a(0),0.d0)
        root(nor) = (-a(1)+csqrt(disc))/2./a(2)
        nor = nor + 1
        root(nor) = (-a(1)-csqrt(disc))/2./a(2)
        do 2 i=1,nor
           write(*,19)i,root(i)
2       continue
        stop
      endif
* Initialize randomly in (-1,1)
      p1 = 2*ranm()-1
      p2 = 2*ranm()-1
      xp(1) = p1
      xp(2) = p2
* Minimize the remainder down to zero !!!
      call OPTIMA(N,M,XP,XV,XLL,XRL,IXAT,ICODE,
     >            'in.dat','out.dat',grms,NF,NG,NHE,NJA)
      if (sqrt(xv) .le. eps) then
        p1 = xp(1)
        p2 = xp(2)
        disc = cmplx( p1**2-4*p2, 0.d0)
        nor = nor + 1
        root(nor)=(-p1 + csqrt(disc))/2
        nor = nor + 1
        root(nor)=(-p1 - csqrt(disc))/2
        call syndiv(a,kco,p1,p2,b,R,Q)
        kco = kco - 2
        do 3 i = 0,kco
          a(i) = b(i)
3       continue
      endif
      go to 30
      end
*
      subroutine syndiv(a,n,p1,p2,b,R,Q)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
* Given a Polynomial: a(n)*x**n + ...+a(1)*x +a(0)
* and a quadratic: x**2 + p1*x + p2
* calculate the coefficients of the quotient polynomial: b(k),k=0,1,...,n-2
* and the remainder terms R and Q.
      dimension a(0:*),b(0:*)

      b(n-2)=a(n)
      b(n-3)=a(n-1)-p1*b(n-2)
      do 1 k = n-2,2,-1
        b(k-2) = a(k)-p1*b(k-1)-p2*b(k)
1     continue
      R = a(1) - p1*b(0)-p2*b(1)
      Q = a(0) - p2*b(0)
      end
*
      subroutine synder(a,n,p1,p2,b,R,Q,b1,b2,R1,R2,Q1,Q2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
* Given a Polynomial: a(n)*x**n + ...+a(1)*x +a(0)
* and a quadratic: x**2 + p1*x + p2
* Calculate the coefficients of the quotient polynomial: b(k),k=0,1,...,n-2
* and the remainder terms R and Q.
* Calculate the derivatives: b1(k),R1,Q1 of b(k),R and Q with respect to p1, 
* and the derivatives: b2(k),R2,Q2 of b(k),R and Q with respect to p2.
      dimension a(0:*),b(0:*),b1(0:*),b2(0:*)
*
      b(n-2)=a(n)
      b(n-3)=a(n-1) - p1*b(n-2)
      do 1 k = n-2,2,-1
        b(k-2) = a(k) - p1*b(k-1) - p2*b(k)
1     continue
      R = a(1) - p1*b(0)-p2*b(1)
      Q = a(0) - p2*b(0)
*
      b1(n-2)= 0
      b1(n-3)= -a(n)
      do 2 k = n-2,2,-1
        b1(k-2) = -b(k-1) - p1*b1(k-1) - p2*b1(k)
2     continue
      R1 = -b(0) - p1*b1(0) - p2*b1(1)
      Q1 = -p2*b1(0)
*
      b2(n-2)= 0
      b2(n-3)= 0
      do 3 k = n-2,2,-1
        b2(k-2) = -p1*b2(k-1) - b(k) - p2*b2(k)
3     continue
      R2 = -p1*b2(0) - b(1) - p2*b2(1)
      Q2 = -b(0) - p2*b2(0)
      end
C  ---------------------------------------------------------------------
      SUBROUTINE SUBSUM ( M, N, X, F )
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N), F(M)
      common/toroot/a(0:100),b(0:100),kco
       p1 = x(1)
       p2 = x(2)
       call syndiv(a,kco,p1,p2,b,R,Q)
       f(1) = R
       f(2) = Q
      END
C  ---------------------------------------------------------------------
      SUBROUTINE GRANAL ( N, X, GRAD )
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N), GRAD(N)
      common/toroot/a(0:100),b(0:100),kco
      dimension b1(0:100),b2(0:100)
       p1 = x(1)
       p2 = x(2)
       call  synder(a,kco,p1,p2,b,R,Q,b1,b2,R1,R2,Q1,Q2)
       grad(1) = 2*(R*R1 + Q*Q1)
       grad(2) = 2*(R*R2 + Q*Q2)
      end
C  ---------------------------------------------------------------------
      SUBROUTINE JANAL ( M, N, X, FJ, LD )
C  ---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N), FJ(LD,N)
      common/toroot/a(0:100),b(0:100),kco
      dimension b1(0:100),b2(0:100)
       p1 = x(1)
       p2 = x(2)
       call  synder(a,kco,p1,p2,b,R,Q,b1,b2,R1,R2,Q1,Q2)
       FJ(1,1) = R1
       FJ(1,2) = R2
       FJ(2,1) = Q1
       FJ(2,2) = Q2
      END
