program main

c*******************************************************************72
c
cc zero_laguerre_test() tests zero_laguerre().
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    26 March 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'zero_laguerre_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test zero_laguerre().'

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

      stop
      end
      subroutine test01 ( )

c*******************************************************************72
c
cc test01() runs the tests on a polynomial function.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    18 January 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision abserr
      integer degree
      double precision, external :: func01
      integer ierror
      integer k
      integer kmax
      double precision x
      double precision x0
c
c  Give a starting point.
c
      x0 = 1.0D+00
c
c  The polynomial degree.
c
      degree = 3
c
c  Set the error tolerance.
c
      abserr = 0.00001D+00
c
c  KMAX is the maximum number of iterations.
c
      kmax = 30

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test01():'
      write ( *, '(a)' ) '  (Polynomial function F(X))'
      write ( *, '(a)' ) '  Find a root of F(X)=(X+3)*(X+3)*(X-2)=0'

      call laguerre ( x0, degree, abserr, kmax, func01, x, ierror, k )

      write ( *, '(a)' ) ' '
      if ( ierror /= 0 ) then
        write ( *, '(a,i2)' ) '  Iteration failed, ierror = ', ierror
      else
        write ( *, '(a,i2)' ) '  Iteration steps taken: ', k
        write ( *, '(a,g14.6)' ) '  Estimated root X = ', x
        write ( *, '(a,g14.6)' ) '  F(X) = ', func01 ( x, 0 )
      end if

      return
      end
      function func01 ( x, ider )

c*******************************************************************72
c
cc func01() computes the function value for the first test.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    17 December 1998
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real X, the point at which the evaluation is to take place.
c
c    integer IDER, specifies what is to be evaluated:
c    0, evaluate the function.
c    1, evaluate the first derivative.
c    2, evaluate the second derivative.
c    3, evaluate the third derivative.
c
c  Output:
c
c    real FUNC01, the value of the function or derivative.
c
      implicit none

      double precision func01
      integer ider
      double precision x

      if ( ider == 0 ) then
        func01 = ( x + 3.0D+00 )**2 * ( x - 2.0E+00 )
      else if ( ider == 1 ) then
        func01 = ( x + 3.0D+00 ) * ( 3.0E+00 * x - 1.0E+00 )
      else if ( ider == 2 ) then
        func01 = 6.0D+00 * x + 8.0E+00
      else if ( ider == 3 ) then
        func01 = 6.0D+00
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'func01(): Fatal errorc'
        write ( *, '(a,i8)' ) '  Derivative of order IDER = ', ider
        write ( *, '(a)' ) '  was requested.'
        stop
      end if

      return
      end
      subroutine test02 ( )

c********************************************************************72
c
cc test02() runs the tests on the Newton polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    26 March 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision abserr
      integer degree
      double precision, external :: func02
      integer ierror
      integer k
      integer kmax
      double precision x
      double precision x0
c
c  Give a starting point.
c
      x0 = 1.0D+00
c
c  The polynomial degree.
c
      degree = 3
c
c  Set the error tolerance.
c
      abserr = 0.00001D+00
c
c  KMAX is the maximum number of iterations.
c
      kmax = 30

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02():'
      write ( *, '(a)' ) '  (Polynomial function F(X))'
      write ( *, '(a)' ) '  p(x) = x^3 - 2x - 5'

      call laguerre ( x0, degree, abserr, kmax, func02, x, ierror, k )

      write ( *, '(a)' ) ' '
      if ( ierror /= 0 ) then
        write ( *, '(a,i2)' ) '  Iteration failed, ierror = ', ierror
      else
        write ( *, '(a,i2)' ) '  Iteration steps taken: ', k
        write ( *, '(a,g14.6)' ) '  Estimated root X = ', x
        write ( *, '(a,g14.6)' ) '  F(X) = ', func02 ( x, 0 )
      end if

      return
      end
      function func02 ( x, ider )

c*********************************************************************72
c
cc func02() computes the function value for the Newton polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    26 March 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    real X, the point at which the evaluation is to take place.
c
c    integer IDER, specifies what is to be evaluated:
c    0, evaluate the function.
c    1, evaluate the first derivative.
c    2, evaluate the second derivative.
c    3, evaluate the third derivative.
c
c  Output:
c
c    real FUNC02, the value of the function or derivative.
c
      implicit none

      double precision func02
      integer ider
      double precision x

      if ( ider == 0 ) then
        func02 = x**3 - 2.0D+00 * x - 5.0D+00
      else if ( ider == 1 ) then
        func02 = 3.0D+00 * x**2 - 2.0D+00
      else if ( ider == 2 ) then
        func02 = 6.0D+00 * x 
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'func02(): Fatal errorc'
        write ( *, '(a,i8)' ) '  Derivative of order IDER = ', ider
        write ( *, '(a)' ) '  was requested.'
        stop
      end if

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
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, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end