program main

c*********************************************************************72
c
cc toms419_test() tests toms419().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
      implicit none

      call timestamp ( )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'toms419_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test toms419().'

      call test01 ( )
      call test02 ( )
      call test03 ( )
      call test04 ( )
      call test05 ( )
c
c  Termination.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'toms419_test():'
      write ( *, '(a)' ) '  Normal end of execution.'

      call timestamp ( )

      stop
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc test01() tests the first polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
      implicit none

      integer, parameter :: degree = 10
      integer, parameter :: order = degree + 1

      logical fail
      double precision fzr(degree)
      double precision fzi(degree)
      integer i
      double precision pr(order)
      double precision pi(order)
      double precision zr(degree)
      double precision zi(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test01():'
      write ( *, '(a)' ) '  Polynomial with zeros 1,2,...,10.'
      write ( *, '(a,i2)' ) '  Degree ', degree
 
      pr(1) = 1.0
      pr(2) = -55.0
      pr(3) = 1320.0
      pr(4) = -18150.0
      pr(5) = 157773.0
      pr(6) = -902055.0
      pr(7) = 3416930.0
      pr(8) = -8409500.0
      pr(9) = 12753576.0
      pr(10) = -10628640.0
      pr(11) = 3628800.0
      do i = 1, order
        pi(i) = 0.0
      end do

      call print_coefficients ( order, pr, pi )
      call cpoly ( pr, pi, degree, zr, zi, fail )

      if ( fail ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  cpoly() failed on this example.'
      else
        call print_zeros ( degree, zr, zi )
        call compute_fz ( degree, zr, zi, pr, pi, fzr, fzi )
        call print_fz ( degree, fzr, fzi )
      end if

      return
      end
      subroutine test02 ( )

c*********************************************************************72
c
cc test02() tests the second polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
      implicit none

      integer, parameter :: degree = 3
      integer, parameter :: order = degree + 1

      logical fail
      double precision fzr(degree)
      double precision fzi(degree)
      double precision pr(order)
      double precision pi(order)
      double precision zr(degree)
      double precision zi(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02():'
      write ( *, '(a)' ) '  Zeros on imaginary axis.'
      write ( *, '(a,i2)' ) '  Degree ', degree

      pr(1) = 1.0
      pr(2) = 0.0
      pr(3) = -10001.0001D0
      pr(4) = 0.0
      pi(1) = 0.0
      pi(2) = -10001.0001d0
      pi(3) = 0.0
      pi(4) = 1.0

      call print_coefficients ( order, pr, pi )
      call cpoly ( pr, pi, degree, zr, zi, fail )

      if ( fail ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  cpoly() failed on this example.'
      else
        call print_zeros ( degree, zr, zi )
        call compute_fz ( degree, zr, zi, pr, pi, fzr, fzi )
        call print_fz ( degree, fzr, fzi )
      end if

      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
cc test03() tests the third polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
      implicit none

      integer, parameter :: degree = 10
      integer, parameter :: order = degree + 1

      logical fail
      double precision fzr(degree)
      double precision fzi(degree)
      double precision pr(order)
      double precision pi(order)
      double precision zr(degree)
      double precision zi(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test03():'
      write ( *, '(a)' ) '  zeros at 1+I,1/2*(1+I)....1/(2**-9)*(1+I)'
      write ( *, '(a,i2)' ) '  Degree ', degree

      pr(1) = 1.0
      pr(2) = -1.998046875
      pr(3) = 0.0
      pr(4) = 0.7567065954208374D0
      pr(5) = -0.2002119533717632D0
      pr(6) = 1.271507365163416D-2
      pr(7) = 0.0
      pr(8) = -1.154642632172909D-5
      pr(9) = 1.584803612786345D-7
      pr(10) = -4.652065399568528D-10
      pr(11) = 0.0
      pi(1) = 0.0
      pi(2) = pr(2)
      pi(3) = 2.658859252929688d0
      pi(4) = -7.567065954208374d-1
      pi(5) = 0.0
      pi(6) = pr(6)
      pi(7) = -7.820779428584501d-4
      pi(8) = -pr(8)
      pi(9) = 0.0
      pi(10) = pr(10)
      pi(11) = 9.094947017729282d-13

      call print_coefficients ( order, pr, pi )
      call cpoly ( pr, pi, degree, zr, zi, fail )

      if ( fail ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  cpoly() failed on this example.'
      else
        call print_zeros ( degree, zr, zi )
        call compute_fz ( degree, zr, zi, pr, pi, fzr, fzi )
        call print_fz ( degree, fzr, fzi )
      end if

      return
      end
      subroutine test04 ( )

c*********************************************************************72
c
cc test04() tests the fourth polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
      implicit none

      integer, parameter :: degree = 10
      integer, parameter :: order = degree + 1

      logical fail
      double precision fzr(degree)
      double precision fzi(degree)
      double precision pr(order)
      double precision pi(order)
      double precision zr(degree)
      double precision zi(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test04():'
      write ( *, '(a)' ) '  Multiple zeros.'
      write ( *, '(a,i2)' ) '  Degree ', degree

      pr(1) = 1.0
      pr(2) = -10.0
      pr(3) = 3.0
      pr(4) = 284.0
      pr(5) = -1293.0
      pr(6) = 2374.0
      pr(7) = -1587.0
      pr(8) = -920.0
      pr(9) = 2204.0
      pr(10) = -1344.0
      pr(11) = 288.0
      pi(1) = 0.0
      pi(2) = -10.0
      pi(3) = 100.0
      pi(4) = -334.0
      pi(5) = 200.0
      pi(6) = 1394.0
      pi(7)  = -3836.0
      pi(8) = 4334.0
      pi(9) = -2352.0
      pi(10) = 504.0
      pi(11) = 0.0

      call print_coefficients ( order, pr, pi )
      call cpoly ( pr, pi, degree, zr, zi, fail )

      if ( fail ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  cpoly() failed on this example.'
      else
        call print_zeros ( degree, zr, zi )
        call compute_fz ( degree, zr, zi, pr, pi, fzr, fzi )
        call print_fz ( degree, fzr, fzi )
      end if

      return
      end
      subroutine test05 ( )

c*********************************************************************72
c
cc test05() tests the fifth polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
      implicit none

      integer, parameter :: degree = 12
      integer, parameter :: order = degree + 1

      logical fail
      double precision fzr(degree)
      double precision fzi(degree)
      double precision pr(order)
      double precision pi(order)
      double precision zr(degree)
      double precision zi(degree)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test05():'
      write ( *, '(a)' ) '  Zeros on unit circle centered at 0+2i.'
      write ( *, '(a,i2)' ) '  Degree ', degree

      pr(1) = 1.0
      pr(2) = 0.0
      pr(3) = -264.0
      pr(4) = 0.0
      pr(5) = 7920.0
      pr(6) = 0.0
      pr(7) = -59136.0
      pr(8) = 0.0
      pr(9) = 126720.0
      pr(10) = 0.0
      pr(11) = -67584.0
      pr(12) = 0.0
      pr(13) = 4095.0
      pi(1) = 0.0
      pi(2) = -24.0
      pi(3) = 0.0
      pi(4) = 1760.0
      pi(5) = 0.0
      pi(6) = -25344.0
      pi(7) = 0.0
      pi(8) = 101376.0
      pi(9) = 0.0
      pi(10) = -112640.0
      pi(11) = 0.0
      pi(12) = 24576.0
      pi(13) = 0.0

      call print_coefficients ( order, pr, pi )
      call cpoly ( pr, pi, degree, zr, zi, fail )

      if ( fail ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  cpoly() failed on this example.'
      else
        call print_zeros ( degree, zr, zi )
        call compute_fz ( degree, zr, zi, pr, pi, fzr, fzi )
        call print_fz ( degree, fzr, fzi )
      end if

      return
      end
      subroutine print_coefficients ( order, pr, pi )

c*********************************************************************72
c
cc print_coefficients() prints the coefficients.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Input:
c
c    integer order: the order of the polynomial.
c
c    double precision pr(order), pi(order): the real and imaginary parts
c    of the coefficients.
c
      implicit none

      integer order

      integer i
      double precision pi(order)
      double precision pr(order)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Coefficients:'
      write ( *, '(50(2d26.16/))' ) ( pr(i), pi(i), i = 1, order )

      return
      end
      subroutine print_zeros ( degree, zr, zi )

c*********************************************************************72
c
cc print_zeros() prints the zeros.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Input:
c
c    integer degree: the degree of the polynomial.
c
c    double precision zr(degree), zi(degree): the real and imaginary parts
c    of the zeros.
c
      implicit none

      integer degree

      integer i
      double precision zi(degree)
      double precision zr(degree)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Zeros:'
      write ( *, '(50(2d26.16/))' ) ( zr(i), zi(i), i = 1, degree )

      return
      end
      subroutine compute_fz ( degree, zr, zi, pr, pi, fzr, fzi )

c*********************************************************************72
c
cc compute_fz() computes the function values at the zeros.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Input:
c
c    integer degree: the degree of the polynomial.
c
c    double precision pr(degree+1), pi(degree+1): the real and imaginary parts
c    of the polynomial coefficients.
c
c    double precision zr(degree), zi(degree): the real and imaginary parts
c    of the computed zeros.
c
c  Output:
c
c    double precision fzr(degree), fzi(degree): the real and imaginary parts
c    of the function values.
c
      implicit none

      integer degree

      integer i
      double precision fzi(degree)
      double precision fzr(degree)
      double precision pi(degree+1)
      double precision pr(degree+1)
      double precision qpi(degree+1)
      double precision qpr(degree+1)
      double precision zi(degree+1)
      double precision zr(degree+1)

      do i = 1, degree
        call polyev ( degree + 1, zr(i), zi(i), pr, pi, qpr, qpi, 
     &    fzr(i), fzi(i) )
      end do

      return
      end
      subroutine print_fz ( degree, fzr, fzi )

c*********************************************************************72
c
cc print_fz() prints the function values at the zeros.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    13 March 2024
c
c  Author:
c
c    Original Fortran77 version by Michael Jenkins, Joseph Traub.
c    This version by John Burkardt.
c
c  Input:
c
c    integer degree: the degree of the polynomial.
c
c    double precision fzr(degree), fzi(degree): the real and imaginary parts
c    of the function values.
c
      implicit none

      integer degree

      integer i
      double precision fzi(degree)
      double precision fzr(degree)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) '  Function values:'
      write ( *, '(50(2d26.16/))' ) ( fzr(i), fzi(i), i = 1, degree )

      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    12 January 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ',
     &  'May      ', 'June     ', 'July     ', 'August   ',
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *,
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' )
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end