program main

c*********************************************************************72
c
cc toms493_test() tests toms493().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 December 2008
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'toms493_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test toms493().'

      call test01 ( )
      call test02 ( )
      call test03 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'toms493_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc test01() uses a polynomial suggested by Driessen and Hunt.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 December 2007
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    HB Driessen, EW Hunt,
c    Remark on Algorithm 429,
c    Communications of the ACM,
c    Volume 16, Number 9, page 579, September 1973.
c
      implicit none

      integer degree
      parameter ( degree = 4 )

      double precision c(degree+1)
      logical fail
      integer i
      double precision pvi
      double precision pvr
      double precision root_i(degree)
      double precision root_r(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test01():'
      write ( *, '(a)' ) '  p(x) = x^4 + 5.6562x^3 + 5.8854x^2'
      write ( *, '(a)' ) '             + 7.3646x   + 6.1354'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Approximate roots:'
      write ( *, '(a)' ) '    -1.001,'
      write ( *, '(a)' ) '    -4.7741, '
      write ( *, '(a)' ) '     0.0089 + 1.1457 i,'
      write ( *, '(a)' ) '     0.0089 - 1.1457 i.'

      c(1) = 1.0D+00
      c(2) = 5.6562D+00
      c(3) = 5.8854D+00
      c(4) = 7.3646D+00
      c(5) = 6.1354D+00

      call rpoly ( c, degree, root_r, root_i, fail )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 
     &  '      real(z)         imag(z)       real(f(z))      imag(f(z))'
      write ( *, '(a)' ) ''

      do i = 1, degree

        call polyev ( degree + 1, root_r(i), root_i(i), c, pvr, pvi )

        write ( *, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)' ) 
     &    root_r(i), root_i(i), pvr, pvi

      end do

      return
      end
      subroutine test02 ( )

c*********************************************************************72
c
cc test02() tests a polynomial with a single root.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    30 December 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer degree
      parameter ( degree = 4 )

      double precision c(degree+1)
      logical fail
      integer i
      double precision pvi
      double precision pvr
      double precision root_i(degree)
      double precision root_r(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02():'
      write ( *, '(a)' ) '  p(x) = x^4 - 8 x^3 + 24 x^2 - 32 x + 16'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Exact roots:'
      write ( *, '(a)' ) '    2 (multiplicity 4)'

      c(1) =   1.0D+00
      c(2) =  -8.0D+00
      c(3) =  24.0D+00
      c(4) = -32.0D+00
      c(5) =  16.0D+00

      call rpoly ( c, degree, root_r, root_i, fail )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 
     &  '      real(z)         imag(z)       real(f(z))      imag(f(z))'
      write ( *, '(a)' ) ''

      do i = 1, degree

        call polyev ( degree + 1, root_r(i), root_i(i), c, pvr, pvi )

        write ( *, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)' ) 
     &    root_r(i), root_i(i), pvr, pvi

      end do

      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
cc test03() tests a polynomial in the roots of unity.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer degree
      parameter ( degree = 5 )

      double precision c(degree+1)
      logical fail
      integer i
      double precision pvi
      double precision pvr
      double precision root_i(degree)
      double precision root_r(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test03():'
      write ( *, '(a)' ) '  p(x) = x^5 - 1'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Exact roots:'
      write ( *, '(a)' ) '    0.3090 + 0.9510 i'
      write ( *, '(a)' ) '   -0.8090 + 0.5877 i'
      write ( *, '(a)' ) '   -0.8090 - 0.5877 i'
      write ( *, '(a)' ) '    0.3090 - 0.9510 i'
      write ( *, '(a)' ) '    1'

      c(1) =  1.0D+00
      c(2) =  0.0D+00
      c(3) =  0.0D+00
      c(4) =  0.0D+00
      c(5) =  0.0D+00
      c(6) = -1.0D+00

      call rpoly ( c, degree, root_r, root_i, fail )

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 
     &  '      real(z)         imag(z)       real(f(z))      imag(f(z))'
      write ( *, '(a)' ) ''

      do i = 1, degree

        call polyev ( degree + 1, root_r(i), root_i(i), c, pvr, pvi )

        write ( *, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)' ) 
     &    root_r(i), root_i(i), pvr, pvi

      end do

      return
      end
      subroutine polyev ( nn, sr, si, c, pvr, pvi )

c*********************************************************************72
c
cc polyev() evaluates a real polynomial by the Horner recurrence.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    14 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Reference:
c
c    Michael Jenkins, Joseph Traub,
c    Algorithm 419: Zeros of a Complex Polynomial,
c    Communications of the Association for Computing Machinery,
c    Volume 15, Number 2, February 1972, pages 97-99.
c
c  Input:
c
c    integer nn: the polynomial order.
c
c    double precison sr, si: the real and imaginary parts of the
c    polynomial argument.
c
c    double precison c(nn): the polynomial coefficients.
c
c  Output:
c
c    double precision pvr, pvi: the real and imaginary parts of the
c    polynomial value.
c
      implicit none

      integer nn

      double precision c(nn)
      integer i
      double precision pvi
      double precision pvr
      double precision si
      double precision sr
      double precision t

      do i = 1, nn
        if ( i == 1 ) then
          pvr = c(1)
          pvi = 0.0
        else
          t =   pvr * sr - pvi * si + c(i)
          pvi = pvr * si + pvi * sr
          pvr = t
        end if
      end do

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 September 2005
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character ( len = 8 ) date
      character ( len = 10 ) time

      call date_and_time ( date, time )

      write ( *, '(a8,2x,a10)' ) date, time

      return
      end