program divdif_prb
!
!*******************************************************************************
!
!! DIVDIF_PRB calls a set of problems for DIVDIF.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DIVDIF_PRB'
  write ( *, '(a)' ) '  A set of tests for DIVDIF.'
 
  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09
  call test10
  call test11
  call test12
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DIVDIF_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests DATA_TO_DIF_DISPLAY;
!! TEST01 tests DIF_APPEND;
!! TEST01 tests DIF_ANTIDERIV;
!! TEST01 tests DIF_DERIV;
!! TEST01 tests DIF_SHIFT_ZERO.
!! TEST01 tests DIF_VAL;
!
  implicit none
!
  integer, parameter :: maxtab = 10
!
  real diftab1(maxtab)
  real diftab2(maxtab)
  real diftab3(maxtab)
  integer i
  integer ntab
  real xtab(maxtab)
  real xval
  real ytab(maxtab)
  real yval
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  For a divided difference polynomial:'
  write ( *, '(a)' ) '  DATA_TO_DIF_DISPLAY sets up a difference table'
  write ( *, '(a)' ) '    and displays intermediate calculations;'
  write ( *, '(a)' ) '  DIF_APPEND appends a new data point;'
  write ( *, '(a)' ) '  DIF_ANTIDERIV computes the antiderivative;'
  write ( *, '(a)' ) '  DIF_DERIV computes the derivative;'
  write ( *, '(a)' ) '  DIF_SHIFT_ZERO shifts all the abscissas to 0;'
  write ( *, '(a)' ) '  DIF_VAL evaluates at a point.'
  write ( *, '(a)' ) ' '
!
!  Set XTAB, YTAB to X, X**2.
!
  ntab = 4
  call rvec_identity ( ntab, xtab )
  ytab(1:ntab) = xtab(1:ntab)**2

  call data_to_dif_display ( diftab1, ntab, xtab, ytab )

  call dif_print ( diftab1, ntab, xtab, '  The divided difference polynomial:' )
!
!  Append (5,25) to the table.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Append the data (5,25) to the table.'
  write ( *, '(a)' ) ' '

  ntab = ntab + 1
  xval = 5.0E+00
  yval = 25.0E+00
 
  call dif_append ( diftab1, ntab, xtab, xval, ytab, yval )

  call dif_print ( diftab1, ntab, xtab, '  The divided difference polynomial:' )
!
!  Evaluate the polynomial at 2.5.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Evaluate the table at a point.'
  write ( *, '(a)' ) ' '

  xval = 2.5E+00
 
  call dif_val ( diftab1, ntab, xtab, xval, yval )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a,g14.6)' ) '  P( ', xval, ' ) = ', yval
!
!  Shift the base to zero.
!
  call dif_shift_zero ( diftab1, ntab, xtab )
 
  call dif_print ( diftab1, ntab, xtab, '  The table, rebased at 0:' )
!
!  Compute the derivative.
!
  call dif_deriv ( diftab1, diftab2, ntab, xtab )
 
  call dif_print ( diftab2, ntab-1, xtab, '  The derivative:' )
!
!  Compute the antiderivative.
!
  call dif_antideriv ( diftab1, diftab3, ntab, xtab )
 
  call dif_print ( diftab3, ntab+1, xtab, '  The antiderivative:' )
 
  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests DATA_TO_DIF;
!! TEST02 tests DIF_VAL.
!
  implicit none
!
  integer, parameter :: maxtab = 5
!
  real diftab(maxtab)
  real error
  integer j
  integer ntab
  real true
  real xtab(maxtab)
  real xval
  real ytab(maxtab)
  real yval
!
  xval = 2.5E+00
  true = exp ( xval )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Approximate Y = EXP(X) using orders 1 to 5'
  write ( *, '(a,g14.6)' ) '  Evaluate at X = ', xval
  write ( *, '(a,g14.6)' ) '  where EXP(X)=   ', true
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Order  Approximate Y     Error'
  write ( *, '(a)' ) ' '
 
  do ntab = 1, maxtab
 
    do j = 1, ntab
      xtab(j) = real ( j - 1 )
      ytab(j) = exp ( xtab(j) )
    end do
 
    call data_to_dif ( diftab, ntab, xtab, ytab )

    call dif_val ( diftab, ntab, xtab, xval, yval )
 
    error = yval - true
    write ( *, ' ( i6, 2g14.6 )' ) ntab, yval, error
 
  end do
 
  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests DIF_ROOT.
!
  implicit none
!
  integer, parameter :: maxtab = 6
!
  real abserr
  real, external :: cubic
  real diftab(maxtab)
  integer iprint
  integer maxstp
  real relerr
  real xroot
  real xtry1
  real xtry2
!
!  Seek a root of the function F(X) = (X+3)*(X+1)*(X-1)
!  given starting estimates of 0.5 and 2.0.
!
  abserr = 0.00001E+00
  iprint = 1
  maxstp = 15
  relerr = 0.00001E+00
  xtry1 = 0.5E+00
  xtry2 = 2.0E+00
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  DIF_ROOT seeks a zero of F(x).'
  write ( *, '(a)' ) '  F(X) = (X+3)*(X+1)*(X-1)'
  write ( *, '(a)' ) ' '

  call dif_root ( abserr, diftab, cubic, iprint, maxstp, maxtab, &
    relerr, xroot, xtry1, xtry2 )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Estimated root = ', xroot
  write ( *, '(a,g14.6)' ) '  F(X) = ', cubic ( xroot )
 
  return
end
function cubic ( x )
!
!*******************************************************************************
!
!! CUBIC is the function whose root is sought.
!
  implicit none
!
  real cubic
  real x
!
  cubic = ( x + 3.0E+00 ) * ( x + 1.0E+00 ) * ( x - 1.0E+00 )
 
  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests DIF_BASIS.
!
  implicit none
!
  integer, parameter :: ntab = 5
!
  real diftab(ntab,ntab)
  integer i
  integer j
  integer, parameter :: nstep = 9
  real xhi
  real xlo
  real xtab(ntab)
  real xval
  real yval
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  DIF_BASIS computes Lagrange basis polynomials'
  write ( *, '(a)' ) '    in difference form.'
  write ( *, '(a)' ) ' '
!
!  Set the base points.
!
  call rvec_identity ( ntab, xtab )

  call rvec_print ( ntab, xtab, '  The base points:' )
!
!  Get the difference tables for the basis polynomials and print them.
!
  call dif_basis ( diftab, ntab, xtab )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'The table of difference vectors defining the basis'
  write ( *, '(a)' ) 'polynomials.  Each column represents a polynomial.'
  write ( *, '(a)' ) ' '
  do i = 1, ntab
    write ( *, '(5g14.6)' ) diftab(i,1:ntab)
  end do
!
!  Evaluate basis polynomial 3 at a set of points.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Evaluate basis polynomial #3 at a set of points.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      X        Y'
  write ( *, '(a)' ) ' '
  xhi = real ( ntab )
  xlo = 1.0E+00
 
  do i = 1, nstep
 
    xval = ( real ( nstep - i ) * xlo + real ( i - 1 ) * xhi ) &
      / real ( nstep - 1 )
 
    call dif_val ( diftab(1,3), ntab, xtab, xval, yval )
 
    write ( *, '(2g14.6)' ) xval, yval
 
  end do
 
  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests POLY_BASIS.
!
  implicit none
!
  integer, parameter :: ntab = 5
!
  real polcof(ntab,ntab)
  integer i
  integer j
  integer, parameter :: nstep = 9
  real xhi
  real xlo
  real xtab(ntab)
  real xval
  real yval
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  POLY_BASIS computes Lagrange basis polynomials'
  write ( *, '(a)' ) '  in standard form.'
  write ( *, '(a)' ) ' '
!
!  Set the base points.
!
  call rvec_identity ( ntab, xtab )
!
!  Get the difference tables for the basis polynomials and print them.
!
  call rpoly_basis ( polcof, ntab, xtab )

  do i = 1, ntab
    write ( *, '(5g14.6)' ) polcof(i,1:ntab)
  end do
!
!  Print basis polynomial 3 in polynomial form.
!
  call rpoly_print ( ntab-1, polcof(1,3), &
    '  Basis polynomial 3 in standard form:' )
!
!  Evaluate basis polynoimial 3 at a set of points.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Evaluate basis polynomial 3 at a set of points.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '      X        Y'
  write ( *, '(a)' ) ' '
  xhi = real ( ntab )
  xlo = 1.0E+00
 
  do i = 1, nstep
 
    xval = ( real ( nstep - i ) * xlo + real ( i - 1 ) * xhi ) &
      / real ( nstep - 1 )
 
    call rpoly_val ( ntab-1, polcof(1,3), xval, yval )
 
    write ( *, '(2g14.6)' ) xval, yval
 
  end do
 
  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests DIF_TO_RPOLY;
!! TEST06 tests DIF_SHIFT_ZERO.
!
  implicit none
!
  integer, parameter :: maxtab = 10
!
  real diftab1(maxtab)
  real diftab2(maxtab)
  integer i
  integer ntab
  real x
  real xtab1(maxtab)
  real xtab2(maxtab)
  real ytab1(maxtab)
  real ytab2(maxtab)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  DIF_TO_RPOLY converts a difference table to a'
  write ( *, '(a)' ) '    polynomial;'
  write ( *, '(a)' ) '  DIF_SHIFT_ZERO shifts a divided difference '
  write ( *, '(a)' ) '    table to all zero abscissas;'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  These are equivalent operations!'
  write ( *, '(a)' ) ' '
!
!  Set XTAB, YTAB to X, F(X)
!
  ntab = 4

  call rvec_identity ( ntab, xtab1 )
  ytab1(1:ntab) = xtab1(1:ntab)**3 - 2.0E+00 * xtab1(1:ntab)**2 &
    + 3.0E+00 * xtab1(1:ntab) - 4.0E+00

  call rvec_identity ( ntab, xtab2 )
  ytab2(1:ntab) = xtab2(1:ntab)**3 - 2.0E+00 * xtab2(1:ntab)**2 &
    + 3.0E+00 * xtab2(1:ntab) - 4.0E+00
!
!  Compute and display the finite difference table.
!
  call data_to_dif_display ( diftab1, ntab, xtab1, ytab1 )
 
  call data_to_dif_display ( diftab2, ntab, xtab2, ytab2 )
!
!  Examine corresponding polynomial.
!
  call dif_print ( diftab1, ntab, xtab1, '  The divided difference polynomial:' )
!
!  Shift to zero using ZERDIF.
!
  call dif_shift_zero ( diftab1, ntab, xtab1 )
 
  call rpoly_print ( ntab-1, diftab1, '  Using DIF_SHIFT_ZERO' )
!
!  Shift to zero using DIFPOL.
!
  call dif_to_rpoly ( diftab2, ntab, xtab2 )
 
  call rpoly_print ( ntab-1, diftab2, '  Using DIF_TO_RPOLY' )
 
  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests NCC_RULE;
!
  implicit none
!
  integer, parameter :: norder = 8
!
  integer i
  real weight(norder)
  real xtab(norder)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  NCC_RULE computes closed Newton Cotes formulas;'
  write ( *, '(a)') ' '

  call ncc_rule ( norder, xtab, weight )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Newton-Cotes Closed Quadrature Rule:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    Abscissa       Weight'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i3,2g14.6)' ) i, xtab(i), weight(i)
  end do

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests NCO_RULE.
!
  implicit none
!
  integer, parameter :: norder = 8
!
  integer i
  real weight(norder)
  real xtab(norder)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  NCO_RULE computes open Newton Cotes formulas.'
  write ( *, '(a)' ) ' '

  call nco_rule ( norder, xtab, weight )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Newton-Cotes Open Quadrature Rule:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    Abscissa       Weight'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i3,2g14.6)' ) i, xtab(i), weight(i)
  end do

  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 tests RPOLY_ANT_COF;
!! TEST09 tests RPOLY_ANT_VAL;
!! TEST09 tests RPOLY_DER_COF;
!! TEST09 tests RPOLY_DER_VAL;
!! TEST09 tests RPOLY_PRINT;
!! TEST09 tests RPOLY_VAL.
!
  implicit none
!
  integer, parameter :: n = 4
!
  integer i
  real poly_cof(0:n)
  real poly_cof2(0:n+1)
  real poly_cof3(0:n-1)
  real xval
  real yval
  real yval2
  real yval3
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  RPOLY_ANT_COF computes the coefficients of the'
  write ( *, '(a)' ) '    antiderivative of a polynomial;'
  write ( *, '(a)' ) '  RPOLY_ANT_VAL evaluates the antiderivative of'
  write ( *, '(a)' ) '    a polynomial;'
  write ( *, '(a)' ) '  RPOLY_DER_COF computes the coefficients of the'
  write ( *, '(a)' ) '    derivative of a polynomial;'
  write ( *, '(a)' ) '  RPOLY_DER_VAL evaluates the derivative of'
  write ( *, '(a)' ) '    a polynomial;'
  write ( *, '(a)' ) '  RPOLY_PRINT prints a polynomial;'
  write ( *, '(a)' ) '  RPOLY_VAL evaluates a polynomial.'

  do i = 0, n
    poly_cof(i) = real ( i + 1 )
  end do

  call rpoly_print ( n, poly_cof, '  Our initial polynomial:' )

  call rpoly_ant_cof ( n, poly_cof, poly_cof2 )

  call rpoly_print ( n+1, poly_cof2, '  The antiderivative polynomial:' )

  call rpoly_der_cof ( n, poly_cof, poly_cof3 )

  call rpoly_print ( n-1, poly_cof3, '  The derivative polynomial:' )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Evaluate the polynomial, antiderivative and'
  write ( *, '(a)' ) '  derivative, using only the original polynomial'
  write ( *, '(a)' ) '  coefficients:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  X   P(X)   Anti_P(X)     P''(X)'
  write ( *, '(a)' ) ' '

  do i = 0, 2

    xval = real ( i )

    call rpoly_val ( n, poly_cof, xval, yval )

    call rpoly_ant_val ( n, poly_cof, xval, yval2 )

    call rpoly_der_val ( n, poly_cof, xval, yval3 )

    write ( *, '(4g14.6)' ) xval, yval, yval2, yval3

  end do

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests ROOTS_TO_DIF;
!! TEST10 tests DIF_TO_RPOLY.
!
  implicit none
!
  integer, parameter :: maxtab = 5
!
  real diftab(maxtab)
  integer ntab
  real xtab(maxtab)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  ROOTS_TO_DIF computes the divided difference'
  write ( *, '(a)' ) '    polynomial with given roots;'
  write ( *, '(a)' ) '  DIF_TO_RPOLY converts it to a standard form'
  write ( *, '(a)' ) '    polynomial.'
  write ( *, '(a)' ) ' '
 
  ntab = 2
  xtab(1) = 3.0E+00

  call roots_to_dif ( diftab, ntab, xtab )

  call dif_to_rpoly ( diftab, ntab, xtab )

  call rpoly_print ( ntab, diftab, '  The polynomial:' )
 
  ntab = 3
  xtab(1) = 3.0E+00
  xtab(2) = 1.0E+00

  call roots_to_dif ( diftab, ntab, xtab )

  call dif_to_rpoly ( diftab, ntab, xtab )

  call rpoly_print ( ntab, diftab, '  The polynomial:' )

  ntab = 4
  xtab(1) = 3.0E+00
  xtab(2) = 1.0E+00
  xtab(3) = 2.0E+00

  call roots_to_dif ( diftab, ntab, xtab )

  call dif_to_rpoly ( diftab, ntab, xtab )

  call rpoly_print ( ntab, diftab, '  The polynomial:' )

  ntab = 5
  xtab(1) = 3.0E+00
  xtab(2) = 1.0E+00
  xtab(3) = 2.0E+00
  xtab(4) = 4.0E+00

  call roots_to_dif ( diftab, ntab, xtab )

  call dif_to_rpoly ( diftab, ntab, xtab )

  call rpoly_print ( ntab, diftab, '  The polynomial:' )

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests ROOTS_TO_RPOLY.
!
  implicit none
!
  integer, parameter :: maxroot = 5
!
  real c(0:maxroot)
  integer nroot
  real roots(maxroot)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  ROOTS_TO_RPOLY computes polynomial coefficients'
  write ( *, '(a)' ) '    from roots.'
  write ( *, '(a)' ) ' '
 
  nroot = 1
  roots(1) = 3.0E+00

  call roots_to_rpoly ( nroot, roots, c )

  call rpoly_print ( nroot, c, '  The polynomial:' )
 
  nroot = 2
  roots(1) = 3.0E+00
  roots(2) = 1.0E+00

  call roots_to_rpoly ( nroot, roots, c )

  call rpoly_print ( nroot, c, '  The polynomial:' )

  nroot = 3
  roots(1) = 3.0E+00
  roots(2) = 1.0E+00
  roots(3) = 2.0E+00

  call roots_to_rpoly ( nroot, roots, c )

  call rpoly_print ( nroot, c, '  The polynomial:' )

  nroot = 4
  roots(1) = 3.0E+00
  roots(2) = 1.0E+00
  roots(3) = 2.0E+00
  roots(4) = 4.0E+00

  call roots_to_rpoly ( nroot, roots, c )

  call rpoly_print ( nroot, c, '  The polynomial:' )

  return
end
subroutine test12
!
!*******************************************************************************
!
!! TEST12 tests POLY_SHIFT.
!
  implicit none
!
  integer, parameter :: n = 2
!
  integer i
  real poly_cof(0:n)
  real scale
  real shift
!
  scale = 2.0E+00
  shift = +3.0E+00
  poly_cof(0) = +6.0E+00
  poly_cof(1) = -1.0E+00
  poly_cof(2) =  2.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  POLY_SHIFT shifts polynomial coefficients.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Polynomial coefficients for argument X'
  write ( *, '(a)' ) ' '

  do i = 0, n
    write ( *, '(i5,g14.6)' ) i, poly_cof(i)
  end do

  call rpoly_shift ( scale, shift, n, poly_cof )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) 'SCALE = ', scale
  write ( *, '(a,g14.6)' ) 'SHIFT = ', shift
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Polynomial coefficients for argument '
  write ( *, '(a)' ) '  Z = SCALE * X + SHIFT'
  write ( *, '(a)' ) ' '

  do i = 0, n
    write ( *, '(i5,g14.6)' ) i, poly_cof(i)
  end do

  return
end