subroutine data_to_dif ( diftab, ntab, xtab, ytab )
!
!*******************************************************************************
!
!! DATA_TO_DIF sets up a divided difference table from raw data.
!
!
!  Discussion:
!
!    Space can be saved by using a single array for both the DIFTAB and
!    YTAB dummy parameters.  In that case, the divided difference table will
!    overwrite the Y data without interfering with the computation.
!
!  Modified:
!
!    11 April 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real DIFTAB(NTAB), the divided difference coefficients
!    corresponding to the input (XTAB,YTAB).
!
!    Input, integer NTAB, the number of pairs of points
!    (XTAB(I),YTAB(I)) which are to be used as data.  The
!    number of entries to be used in DIFTAB, XTAB and YTAB.
!
!    Input, real XTAB(NTAB), the X values at which data was taken.
!    These values must be distinct.
!
!    Input, real YTAB(NTAB), the corresponding Y values.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  integer j
  logical rvec_distinct
  real xtab(ntab)
  real ytab(ntab)
!
  if ( .not. rvec_distinct ( ntab, xtab ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DATA_TO_DIF - Fatal error!'
    write ( *, '(a)' ) '  Two entries of XTAB are equal!'
    return
  end if
!
!  Copy the data values into DIFTAB.
!
  diftab(1:ntab) = ytab(1:ntab)
!
!  Compute the divided differences.
!
  do i = 2, ntab
    do j = ntab, i, -1

      diftab(j) = ( diftab(j) - diftab(j-1) ) / ( xtab(j) - xtab(j+1-i) )

    end do
  end do

  return
end
subroutine data_to_dif_display ( diftab, ntab, xtab, ytab )
!
!*******************************************************************************
!
!! DATA_TO_DIF_DISPLAY computes a divided difference table and shows how.
!
!
!  Modified:
!
!    02 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real DIFTAB(NTAB), the divided difference
!    coefficients corresponding to the input (XTAB,YTAB).
!
!    Input, integer NTAB, the number of pairs of points
!    (XTAB(I),YTAB(I)) which are to be used as data.  The
!    number of entries to be used in DIFTAB, XTAB and YTAB.
!
!    Input, real XTAB(NTAB), the X values at which data was taken.
!
!    Input, real YTAB(NTAB), the corresponding Y values.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  integer j
  logical rvec_distinct
  real xtab(ntab)
  real ytab(ntab)
!
  if ( .not. rvec_distinct ( ntab, xtab ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DATA_TO_DIF_DISPLAY - Fatal error!'
    write ( *, '(a)' ) '  Two entries of XTAB are equal!'
    return
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Divided difference table:'
  write ( *, '(a)' ) ' '
  write ( *, '(4x,5g14.6)' ) xtab(1:ntab)
  write ( *, '(a)' ) ' '
  write ( *, '(i3,1x,5g14.6)' ) 0, ytab(1:ntab)
!
!  Copy the data values into DIFTAB.
!
  diftab(1:ntab) = ytab(1:ntab)
!
!  Compute the divided differences.
!
  do i = 2, ntab
    do j = ntab, i, -1

      diftab(j) = ( diftab(j) - diftab(j-1) ) / ( xtab(j) - xtab(j+1-i) )

    end do

    write ( *, '(i3,1x,5g14.6)' ) i-1, diftab(i:ntab)

  end do

  return
end
subroutine data_to_rpoly ( c, ntab, xtab, ytab )
!
!*******************************************************************************
!
!! DATA_TO_RPOLY computes the coefficients of a polynomial interpolating data.
!
!
!  Discussion:
!
!    Space can be saved by using a single array for both the C and
!    YTAB parameters.  In that case, the coefficients will
!    overwrite the Y data without interfering with the computation.
!
!  Modified:
!
!    16 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real C(0:NTAB-1), the coefficients of the polynomial that passes
!    through the data (XTAB,YTAB).  C(0) is the constant term.
!
!    Input, integer NTAB, the number of data points.
!
!    Input, real XTAB(NTAB), YTAB(NTAB), the data values.
!
  implicit none
!
  integer ntab
!
  real c(0:ntab-1)
  integer i
  integer j
  logical rvec_distinct
  real xtab(ntab)
  real ytab(ntab)
!
  if ( .not. rvec_distinct ( ntab, xtab ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DATA_TO_POLY - Fatal error!'
    write ( *, '(a)' ) '  Two entries of XTAB are equal!'
    return
  end if

  call data_to_dif ( c, ntab, xtab, ytab )

  call dif_to_poly ( c, ntab, xtab )

  return
end
subroutine dif_antideriv ( diftab, anttab, ntab, xtab )
!
!*******************************************************************************
!
!! DIF_ANTIDERIV computes the antiderivative of a  divided difference polynomial.
!
!
!  Discussion:
!
!    The routine uses the divided difference representation (XTAB, DIFTAB)
!    of a polynomial to compute the divided difference representation
!    (XTAB, ANTTAB) of the antiderivative of the polynomial.
!
!  Definition:
!
!    The antiderivative of a polynomial P(X) is any polynomial Q(X)
!    with the property that d/dX Q(X) = P(X).
!
!    This routine chooses the antiderivative whose constant term is zero.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real DIFTAB(NTAB).
!
!    On input, (XTAB,DIFTAB) is a divided difference table, which may
!    be thought of as representing a polynomial through data points
!    (XTAB,YTAB).
!
!    On output, (XTAB,DIFTAB) represent the same polynomial, but the
!    numeric values in DIFTAB will have changed, because XTAB is changed.
!
!    Output, real ANTTAB(NTAB+1), the divided difference table computed for
!    the antiderivative of the polynomial represented by (XTAB,DIFTAB).
!    ANTTAB(1) is the arbitrary integration constant, and is set to zero,
!    but may be reset to any desired value.  The entries of ANTTAB have
!    been computed with X data based at zero.
!
!    Input, integer NTAB, the dimension of the array DIFTAB.
!
!    Input/output, real XTAB(NTAB+1).  (Note the extra entry!)
!
!    On input, entries 1 through NTAB of XTAB contain the X locations where
!    data was taken.  The NTAB+1 entry of XTAB may have any value.
!
!    On output, all the entries of XTAB have been set to zero.  This
!    corresponds to the changes in DIFTAB, and to the way ANTTAB was
!    computed.
!
  implicit none
!
  integer ntab
!
  real anttab(ntab+1)
  real diftab(ntab)
  real xtab(ntab+1)
!
!  Recompute XTAB and DIFTAB, so that XTAB(1:NTAB) is entirely zero.
!
  call dif_shift_zero ( diftab, ntab, xtab )
!
!  Append a final zero to XTAB.
!
  xtab(ntab+1) = 0.0E+00
!
!  Get the antiderivative of the standard form polynomial.
!
  call rpoly_ant_cof ( ntab-1, diftab, anttab )

  return
end
subroutine dif_append ( diftab, ntab, xtab, xval, ytab, yval )
!
!*******************************************************************************
!
!! DIF_APPEND adds a pair of data values to a divided difference table.
!
!
!  Modified:
!
!    27 March 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real DIFTAB(NTAB), the difference table.
!    On output, DIFTAB has been updated to include the new data.
!
!    Input, integer NTAB, the number of divided differences to be stored in
!    DIFTAB, XTAB, and YTAB.  This includes the NTAB-1 values stored
!    before input, and the extra one to be added by ADDDIF.  If NTAB
!    is 1, then no data is assumed to be prestored in the arrays on
!    input.
!
!    Input/output, real XTAB(NTAB).
!    On input, XTAB(1:NTAB-1) contain the X data values.
!    On output, XTAB(1) is set to XVAL, and XTAB(2:NTAB) contain the
!    original X data.
!
!    Input, real XVAL, the X data value to be inserted as XTAB(1).
!
!    Input/output, real YTAB(NTAB).
!    On input, YTAB(1:NTAB-1) are the Y values corresponding to XTAB(1:NTAB-1).
!    On output, YTAB(1:NTAB) are the Y values corresponding to XTAB(1:NTAB).
!
!    Input, real YVAL, the Y data value to be inserted as YTAB(1).
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  real xtab(ntab)
  real xval
  real ytab(ntab)
  real yval
!
!  Move the original data up one index.
!
  diftab(2:ntab) = diftab(1:ntab-1)
  xtab(2:ntab) = xtab(1:ntab-1)
  ytab(2:ntab) = ytab(1:ntab-1)
!
!  Insert the new data.
!
  diftab(1) = yval
  xtab(1) = xval
  ytab(1) = yval
!
!  Recompute the difference table.
!
  do i = 2, ntab
    diftab(i) = ( diftab(i) - diftab(i-1) ) / ( xtab(i) - xtab(1) )
  end do

  return
end
subroutine dif_basis ( diftab, ntab, xtab )
!
!*******************************************************************************
!
!! DIF_BASIS computes all Lagrange basis polynomials in divided difference form.
!
!
!  Definition:
!
!    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
!    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
!    XTAB(J) for J not equal to I, and 1 when J is equal to I.
!
!    The Lagrange basis polynomials have the property that the interpolating
!    polynomial through a set of NTAB data points (XTAB,YTAB) may be
!    represented as
!
!      P(X) = Sum ( I = 1 to N ) YTAB(I) * L(I,NTAB,XTAB)(X)
!
!    Higher order interpolation at selected points may be accomplished
!    using repeated X values, and scaled derivative values.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real DIFTAB(NTAB,NTAB), the set of divided difference tables.
!    Column I of DIFTAB contains the table for the I-th Lagrange basis
!    polynomial.
!
!    Input, integer NTAB, the number of X data points XTAB, and the number of
!    basis polynomials to compute.
!
!    Input, real XTAB(NTAB), the X values upon which the Lagrange basis
!    polynomials are to be based.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab,ntab)
  integer i
  integer j
  real xtab(ntab)
!
!  Initialize DIFTAB to the identity matrix.
!
  diftab(1:ntab,1:ntab) = 0.0E+00
  do i = 1, ntab
    diftab(i,i) = 1.0E+00
  end do
!
!  Compute each Lagrange basis polynomial.
!
  do i = 1, ntab
    call data_to_dif ( diftab(1,i), ntab, xtab, diftab(1,i) )
  end do

  return
end
subroutine dif_basis_i ( diftab, ival, ntab, xtab )
!
!*******************************************************************************
!
!! DIF_BASIS_I computes the I-th Lagrange basis polynomial in divided difference form.
!
!
!  Definition:
!
!    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
!    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
!    XTAB(J) for J not equal to I, and 1 when J is equal to I.
!
!    The Lagrange basis polynomials have the property that the interpolating
!    polynomial through a set of NTAB data points (XTAB,YTAB) may be
!    represented as
!
!      P(X) = Sum ( I = 1 to N ) YTAB(I) * L(I,NTAB,XTAB)(X)
!
!    Higher order interpolation at selected points may be accomplished
!    using repeated X values, and scaled derivative values.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real DIFTAB(NTAB), the divided difference table for the IVAL-th
!    Lagrange basis polynomial.
!
!    Input, integer IVAL, the index of the desired Lagrange basis polynomial.
!    IVAL should be between 1 and NTAB.
!
!    Input, integer NTAB, the number of data points XTAB.
!
!    Input, real XTAB(NTAB), the X values upon which the Lagrange basis
!    polynomial is to be based.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  integer ival
  real xtab(ntab)
!
!  Check IVAL.
!
  if ( ival < 1 .or. ival > ntab ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DIF_BASIS_I - Fatal error!'
    write ( *, '(a,i6)' ) '  IVAL must be between 1 and ', ntab
    write ( *, '(a,i6)' ) '  but your value is ', ival
    stop
  end if
!
!  Initialize DIFTAB to Delta(I,J).
!
  diftab(1:ntab) = 0.0E+00
  diftab(ival) = 1.0E+00
!
!  Compute the IVAL-th Lagrange basis polynomial.
!
  call data_to_dif ( diftab, ntab, xtab, diftab )

  return
end
subroutine dif_deriv ( diftab, dertab, ntab, xtab )
!
!*******************************************************************************
!
!! DIF_DERIV computes the derivative of a polynomial in divided difference form.
!
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real DIFTAB(NTAB), the divided difference
!    table for the polynomial whose derivative is desired.
!
!    On output, DIFTAB has been updated so that it is zero
!    based.  However, DIFTAB is still the divided difference
!    table for the same polynomial.
!
!    Output, real DERTAB(NTAB-1), contains the zero based divided
!    difference table for the polynomial which is the derivative of the
!    polynomial represented by DIFTAB.
!
!    Input, integer NTAB, the dimension of the arrays DIFTAB, DERTAB, and XTAB.
!
!    Input/output, real XTAB(NTAB).
!
!    On input, XTAB contains the abscissas for the divided difference table
!    stored in DIFTAB.
!
!    On output, XTAB has been set to all zeroes, and DIFTAB has been updated,
!    and DERTAB computed so that the original polynomial is still represented
!    by (XTAB(1:NTAB), DIFTAB(1:NTAB)), and its derivative by
!    (XTAB(1:NTAB-1), DERTAB(1:NTAB-1)).
!
  implicit none
!
  integer ntab
!
  real dertab(ntab-1)
  real diftab(ntab)
  integer i
  real xtab(ntab)
!
  call dif_shift_zero ( diftab, ntab, xtab )

  do i = 1, ntab-1
    dertab(i) = real ( i ) * diftab(i+1)
  end do

  return
end
subroutine dif_print ( diftab, ntab, xtab, title )
!
!*******************************************************************************
!
!! DIF_PRINT prints the polynomial represented by a divided difference table.
!
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real DIFTAB(NTAB), the divided difference table
!    for the polynomial.
!
!    Input, integer NTAB, the dimension of the arrays DIFTAB and XTAB.
!
!    Input, real XTAB(NTAB), the X values for the polynomial.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  character ( len = * ) title
  real xtab(ntab)
!
  if ( len_trim ( title ) > 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '
  write ( *, '( '' p(x) =                      '', g14.6 )' ) diftab(1)

  do i = 2, ntab
    write ( *, '( '' + ( x - '', g14.6, '') * ( '', g14.6 )' ) &
      xtab(i-1), diftab(i)
  end do

  write ( *, '(80a1)' ) ( ')', i = 1, ntab-1 )

  return
end
subroutine dif_root ( abserr, diftab, fxname, iprint, maxstp, maxtab, &
  relerr, xroot, xtry1, xtry2 )
!
!*******************************************************************************
!
!! DIF_ROOT seeks a zero of F(X) using divided difference techniques.
!
!
!  Discussion:
!
!    The method uses the idea of considering the related function
!
!      H(X) = 1 / F(X)
!
!    The iteration begins with two estimates for the root supplied by
!    the user.
!
!    From the most recent approximation to the root, X(K), the next
!    approximation X(K+1) is determined by:
!
!      X(K+1) = X(K) + H(X(K-R),...,X(K-1)) / H(X(K-R),...,X(K-1),X(K))
!
!    where K-R = 1 until the maximal order NTAB is reached.
!
!    Generally, the next iterate X(K+1) is the zero of a rational function
!    which passes through the previous data points.
!
!  Reference:
!
!    F M Larkin,
!    Root Finding by Divided Differences,
!    Numerische Mathematik,
!    Volume 37, pages 93-104, 1981.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ABSERR, a positive absolute error tolerance.  If an estimate
!    X for the root is found with ABS ( F(X) ) <= ABSERR, the iteration
!    is stopped.
!
!    Workspace, real DIFTAB(MAXTAB).
!
!    Input, external FXNAME, the name of the function routine which evaluates
!    F(X).  The form of FXNAME must be similar to the following function which
!    has F(X) = ( X - 1 ) * ( X + 1 ).
!
!    function parab ( x )
!
!      real parab
!      real x
!
!      parab = ( x - 1.0E+00 ) * ( x + 1.0E+00 )
!
!      return
!    end
!
!    Input, integer IPRINT, a switch controlling printed output:
!    0, only print error messages.
!    nonzero, also print a table of the iterative process.
!
!    Input, integer MAXSTP, the limit on how many iterations may be tried.
!
!    Input, integer MAXTAB, the limit on how high an order can be used in
!    the divided difference table.  MAXTAB must be at least 2, and probably
!    should not be too large.  Perhaps a value of 5 or 6 is reasonable,
!    20 is too large.
!
!    Input, real RELERR, a tolerance on the size of the change in the root
!    estimates.  If a step is taken, and the change in the root estimate is
!    less than RELERR, the iteration will stop.
!
!    Output, real XROOT, the point which the program has produced as an
!    approximate root.
!
!    Either ABS ( F(XROOT) ) <= ABSERR, or the maximum number of steps was
!    reached, or the current estimate of the root could not be significantly
!    improved.
!
!    Input, real XTRY1, XTRY2, two initial approximations to the root,
!    supplied by the user, which must be distinct.
!
  implicit none
!
  integer maxtab
!
  real abserr
  real diftab(maxtab)
  real froot
  real ftemp1
  real ftemp2
  real, external :: fxname
  integer iprint
  integer istep
  integer maxstp
  integer ntab
  real relerr
  real temp
  real xdelt
  real xold
  real xroot
  real xtab(maxtab)
  real xtry1
  real xtry2
  real yval
!
!  Make sure XTRY1 and XTRY2 are not equal.
!
  if ( xtry1 == xtry2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DIF_ROOT - Fatal error!'
    write ( *, '(a)' ) '  XTRY1 = XTRY2 on input.'
    stop
  end if
!
!  Make sure MAXTAB is at least 2.
!
  if ( maxtab < 2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DIF_ROOT - Fatal error!'
    write ( *, '(a)' ) '  MAXTAB < 2 on input!'
    stop
  end if

  xtab(1) = xtry1
  xtab(2) = xtry2
  ftemp1 = fxname ( xtry1 )
  ftemp2 = fxname ( xtry2 )

  if ( abs ( ftemp1 ) > abs ( ftemp2 ) ) then
    xtab(1) = xtry2
    xtab(2) = xtry1
    call r_swap ( ftemp1, ftemp2 )
  end if
!
!  Initialize the number of steps.
!
  istep = 0
!
!  Initialize the number of data points.
!
  ntab = 2

  if ( iprint > 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '   Step  NTAB    XROOT        F(XROOT)      XDELT'
    write ( *, '(a)' ) ' '
  end if
!
!  Initialize the divided difference table data.
!
  diftab(1) = 1.0E+00 / ftemp1
  diftab(2) = 1.0E+00 / ftemp2

  call data_to_dif ( diftab, ntab, xtab, diftab )
!
!  Initialize values used in the iteration.
!
  xroot = xtry1
  froot = ftemp1
  xdelt = xtry1 - xtry2
!
!  Does the starting data already satisfy the function norm
!  error tolerance ABSERR, or the interval norm error tolerance
!  RELERR?
!
  do

    if ( iprint > 0 ) then
      write ( *, '(i4, i2, 3g14.6)' ) istep, ntab, xroot, froot, xdelt
    end if

    if ( abs ( froot ) <= abserr ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIF_ROOT - Absolute convergence,'
      write ( *, '(a)' ) '  The function value meets the error tolerance.'
      exit
    end if

    if ( abs ( xdelt ) <= relerr ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIF_ROOT - Relative convergence.'
      write ( *, '(a)' ) '  The stepsize meets the error tolerance.'
      exit
    end if
!
!  Check the number of steps taken.
!
    if ( istep >= maxstp ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIF_ROOT - Nonconvergence!'
      write ( *, '(a)' ) '  The maximum number of steps was taken.'
      exit
    end if
!
!  Generate the next point, XVAL.
!
    xold = xroot
    istep = istep + 1
!
!  Algorithm breakdown: The divisor DIFTAB(NTAB) is zero.
!
    if ( diftab(ntab) == 0.0E+00 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DIF_ROOT - Fatal error!'
      write ( *, '(a,i6)' ) '  Algorithm using differences of order ', ntab
      write ( *, '(a)' ) '  A zero-divisor was computed.'
      write ( *, '(a)' ) '  The algorithm has broken down.'
      write ( *, '(a)' ) '  Examine the results.  They may be useful.'
      write ( *, '(a)' ) '  Perhaps a lower value of MAXTAB would help.'
      exit
    end if

    xroot = xtab(ntab) + ( diftab(ntab-1) / diftab(ntab) )
    xdelt = xroot - xold
    froot = fxname ( xroot )

    if ( abs ( froot ) <= abserr ) then
      cycle
    end if

    yval = 1.0E+00 / froot
!
!  If we are now using MAXTAB points, we have to remove an old
!  one before adding the new one.
!
    if ( ntab < maxtab ) then
      ntab = ntab + 1
    end if

    call dif_append ( diftab, ntab, xtab, xroot, diftab, yval )

  end do

  return
end
subroutine dif_shift_x ( diftab, ntab, xtab, xval )
!
!*******************************************************************************
!
!! DIF_SHIFT_X replaces one abscissa of a divided difference table with a new one.
!
!
!  Discussion:
!
!    The routine shifts the representation of a divided difference polynomial by
!    dropping the last X value in XTAB, and adding a new X value to the
!    beginning of the XTAB array, suitably modifying the coefficients stored
!    in DIFTAB.
!
!    The representation of the polynomial is changed, but the polynomial itself
!    is should be identical.
!
!  Modified:
!
!    27 March 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real DIFTAB(NTAB), the divided difference coefficients
!    corresponding to the XTAB array.  On output, this array has been
!    adjusted.
!
!    Input, integer NTAB, the number of divided difference coefficients, and
!    the number of entries in XTAB.
!
!    Input/output, real XTAB(NTAB), the X values used in the representation of
!    the divided difference polynomial.
!
!    After a call to this routine, the last entry of XTAB has been dropped,
!    the other entries have shifted up one index, and XVAL has been inserted
!    at the beginning of the array.
!
!    Input, real XVAL, a new X value which is to be used in the representation
!    of the polynomial.  On output, XTAB(1) equals XVAL and the representation
!    of the polynomial has been suitably changed.
!
!    Note that XVAL does not have to be distinct from any of the original XTAB
!    values.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  real xtab(ntab)
  real xval
!
!  Recompute the divided difference coefficients.
!
  do i = ntab-1, 1, -1
    diftab(i) = diftab(i) + ( xval - xtab(i) ) * diftab(i+1)
  end do
!
!  Shift the X values up one position.
!
  xtab(2:ntab) = xtab(1:ntab-1)

! do i = ntab, 2, -1
!   xtab(i) = xtab(i-1)
! end do

  xtab(1) = xval

  return
end
subroutine dif_shift_zero ( diftab, ntab, xtab )
!
!*******************************************************************************
!
!! DIF_SHIFT_ZERO shifts a divided difference table so that all abscissas are zero.
!
!
!  Discussion:
!
!    When the abscissas are changed, the coefficients naturally
!    must also be changed.
!
!    The resulting pair (XTAB, DIFTAB) still represents the
!    same polynomial, but the entries in DIFTAB are now the
!    standard polynomial coefficients.
!
!  Modified:
!
!    13 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real DIFTAB(NTAB), the divided difference table
!    for the polynomial.  On output, DIFTAB is also
!    the coefficient array for the standard representation
!    of the polynomial.
!
!    Input, integer NTAB, the length of the DIFTAB and XTAB arrays.
!
!    Input/output, real XTAB(NTAB), the X values that correspond to the
!    divided difference table.  On output, XTAB contains only zeroes.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  real xtab(ntab)
  real xval
!
  xval = 0.0E+00

  do i = 1, ntab
    call dif_shift_x ( diftab, ntab, xtab, xval )
  end do

  return
end
subroutine dif_to_rpoly ( diftab, ntab, xtab )
!
!*******************************************************************************
!
!! DIF_TO_RPOLY converts a divided difference polynomial to standard form.
!
!
!  Discussion:
!
!    The vector DIFTAB, containing the divided difference polynomial
!    coefficients is overwritten with the standard form polynomial
!    coefficients, but the abscissa vector XTAB is unchanged.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real DIFTAB(NTAB).
!
!    On input, DIFTAB contains the coefficients of the divided difference
!    polynomials, corresponding to the XTAB array.
!
!    On output, DIFTAB contains the standard form polyomial coefficients.
!    DIFTAB(1) is the constant term, and DIFTAB(NTAB) is the coefficient
!    of X**(NTAB-1).
!
!    Input, integer NTAB, the number of coefficients, and abscissas.
!
!    Input, real XTAB(NTAB), the X values used in the divided difference
!    representation of the polynomial.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  integer j
  real xtab(ntab)
!
!  Recompute the divided difference coefficients.
!
  do j = 1, ntab-1
    do i = 1, ntab-j
      diftab(ntab-i) = diftab(ntab-i) - xtab(ntab-i-j+1) * diftab(ntab-i+1)
    end do
  end do

  return
end
subroutine dif_val ( diftab, ntab, xtab, xval, yval )
!
!*******************************************************************************
!
!! DIF_VAL evaluates a divided difference polynomial at a point.
!
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real DIFTAB(NTAB), the divided difference polynomial coefficients.
!
!    Input, integer NTAB, the number of divided difference
!    coefficients in DIFTAB, and the number of points XTAB.
!
!    Input, real XTAB(NTAB), the X values upon which the
!    divided difference polynomial is based.
!
!    Input, real XVAL, the value where the polynomial is to be evaluated.
!
!    Output, real YVAL, the value of the polynomial at XVAL.
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  real xtab(ntab)
  real xval
  real yval
!
  yval = diftab(ntab)
  do i = 1, ntab-1
    yval = diftab(ntab-i) + ( xval - xtab(ntab-i) ) * yval
  end do

  return
end
subroutine nc_rule ( norder, a, b, xtab, weight )
!
!*******************************************************************************
!
!! NC_RULE computes the coefficients of a Newton-Cotes quadrature rule.
!
!
!  Definition:
!
!    For the interval [A,B], the Newton-Cotes quadrature rule estimates
!
!      Integral ( A <= X <= B ) F(X) dX
!
!    using NORDER equally spaced abscissas XTAB(I) and a weight vector
!    WEIGHT(I):
!
!      Sum ( I = 1 to N ) WEIGHT(I) * F ( XTAB(I) ).
!
!    For the CLOSED rule, the abscissas include the points A and B.
!    For the OPEN rule, the abscissas do not include A and B.
!
!  Modified:
!
!    25 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NORDER, the order of the rule.
!
!    Input, real A, B, the left and right endpoints of the interval
!    over which the quadrature rule is to be applied.
!
!    Input, real XTAB(NORDER), the abscissas of the rule.
!
!    Output, real WEIGHT(NORDER), the weights of the rule.
!
  implicit none
!
  integer norder
!
  real a
  real b
  real poly_cof(norder)
  integer i
  real weight(norder)
  real xtab(norder)
  real yvala
  real yvalb
!
  do i = 1, norder
!
!  Compute the Lagrange basis polynomial which is 1 at XTAB(I),
!  and zero at the other nodes.
!
    call rpoly_basis_1 ( poly_cof, i, norder, xtab )
!
!  Evaluate the antiderivative of the polynomial at the left and
!  right endpoints.
!
    call rpoly_ant_val ( norder-1, poly_cof, a, yvala )

    call rpoly_ant_val ( norder-1, poly_cof, b, yvalb )

    weight(i) = yvalb - yvala

  end do

  return
end
subroutine ncc_rule ( norder, xtab, weight )
!
!*******************************************************************************
!
!! NCC_RULE computes the coefficients of a Newton-Cotes closed quadrature rule.
!
!
!  Definition:
!
!    For the interval [-1,1], the Newton-Cotes quadrature rule estimates
!
!      Integral ( -1 <= X <= 1 ) F(X) dX
!
!    using NORDER equally spaced abscissas XTAB(I) and a weight vector
!    WEIGHT(I):
!
!      Sum ( I = 1 to N ) WEIGHT(I) * F ( XTAB(I) ).
!
!    For the CLOSED rule, the abscissas include A and B.
!
!  Modified:
!
!    25 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NORDER, the order of the rule.
!
!    Output, real XTAB(NORDER), the abscissas of the rule.
!
!    Output, real WEIGHT(NORDER), the weights of the rule.
!
  implicit none
!
  integer norder
!
  real a
  real b
  integer i
  real weight(norder)
  real xtab(norder)
!
!  Compute a closed quadrature rule.
!
  a = -1.0E+00
  b = 1.0E+00

  do i = 1, norder
    xtab(i) = ( real ( norder - i ) * a + real ( i - 1 ) * b ) &
      / real ( norder - 1 )
  end do

  call nc_rule ( norder, a, b, xtab, weight )

  return
end
subroutine nco_rule ( norder, xtab, weight )
!
!*******************************************************************************
!
!! NCO_RULE computes the coefficients of a Newton-Cotes open quadrature rule.
!
!
!  Definition:
!
!    For the interval [-1,1], the Newton-Cotes quadrature rule estimates
!
!      Integral ( -1 <= X <= 1 ) F(X) dX
!
!    using NORDER equally spaced abscissas XTAB(I) and a weight vector
!    WEIGHT(I):
!
!      Sum ( I = 1 to N ) WEIGHT(I) * F ( XTAB(I) ).
!
!    For the OPEN rule, the abscissas do not include A and B.
!
!  Modified:
!
!    25 September 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NORDER, the order of the rule.
!
!    Output, real XTAB(NORDER), the abscissas of the rule.
!
!    Output, real WEIGHT(NORDER), the weights of the  rule.
!
  implicit none
!
  integer norder
!
  real a
  real b
  integer i
  real weight(norder)
  real xtab(norder)
!
  a = -1.0E+00
  b = 1.0E+00

  do i = 1, norder
    xtab(i) = ( real ( norder + 1 - i ) * a + real ( i ) * b ) &
      / real ( norder + 1 )
  end do

  call nc_rule ( norder, a, b, xtab, weight )

  return
end
subroutine r_swap ( x, y )
!
!*******************************************************************************
!
!! R_SWAP swaps two real values.
!
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none
!
  real x
  real y
  real z
!
  z = x
  x = y
  y = z

  return
end
subroutine roots_to_dif ( diftab, ntab, xtab )
!
!*******************************************************************************
!
!! ROOTS_TO_DIF sets up a divided difference table for a polynomial from its roots.
!
!
!  Discussion:
!
!    This turns out to be a simple task, because of two facts:
!
!    * The divided difference polynomial of one smaller degree which
!      passes through the values ( ROOT(I), 0 ) is the zero polynomial,
!      and hence has a zero divided difference table.
!
!    * We want a polynomial of one degree higher, but we don't want it
!      to pass through an addditional point.  Instead, we specify that
!      the polynomial is MONIC.  This means that the divided difference
!      table is almost the same as for the zero polynomial, except that
!      there is one more pair of entries, an arbitrary X value, and
!      a Y value of 1.
!
!  Modified:
!
!    18 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real DIFTAB(NTAB), the divided difference coefficients
!    for the polynomial.
!
!    Input, integer NTAB, is the size of the divided difference table.
!    This must be one more than the number of roots.
!
!    Input/output, real XTAB(NTAB).
!    On input, XTAB contains in the first NTAB-1 entries the roots of
!    the polynomial.
!    On output, an extra value has been appended in XTAB(NTAB).
!
  implicit none
!
  integer ntab
!
  real diftab(ntab)
  integer i
  integer j
  logical rvec_distinct
  real xtab(ntab)
!
!  Build the appropriate difference table for the polynomial
!  through ( ROOTS(I), 0 ) of degree NTAB-2.
!
  diftab(1:ntab-1) = 0.0E+00
!
!  Append the extra data to make a monic polynomial of degree NTAB-1
!  which is zero at the NTAB-1 roots.
!
  xtab(ntab) = 1492.0E+00
  diftab(ntab) = 1.0E+00

  return
end
subroutine roots_to_rpoly ( n, x, c )
!
!*******************************************************************************
!
!! ROOTS_TO_RPOLY converts polynomial roots to polynomial coefficients.
!
!
!  Modified:
!
!    17 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of roots specified.
!
!    Input, real X(N), the roots.
!
!    Output, real C(0:N), the coefficients of the polynomial.
!
  implicit none
!
  integer n
!
  real c(0:n)
  integer i
  integer j
  real x(n)
!
!  Initialize C to (0, 0, ..., 0, 1).
!  Essentially, we are setting up a divided difference table.
!
  c(0:n-1) = 0.0E+00
  c(n) = 1.0E+00
!
!  Convert to standard polynomial form by shifting the abscissas
!  of the divided difference table to 0.
!
  do j = 1, n
    do i = 1, n+1-j
      c(n-i) = c(n-i) - x(n+1-i-j+1) * c(n-i+1)
    end do
  end do

  return
end
subroutine rpoly_ant_cof ( n, poly_cof, poly_cof2 )
!
!*******************************************************************************
!
!! RPOLY_ANT_COF integrates a polynomial in standard form.
!
!
!  Definition:
!
!    The antiderivative of a polynomial P(X) is any polynomial Q(X)
!    with the property that d/dX Q(X) = P(X).
!
!    This routine chooses the antiderivative whose constant term is zero.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the original polynomial.
!
!    Input, real POLY_COF(0:N), the polynomial coefficients.
!    POLY_COF(0) is the constant term, and POLY_COF(N) is the
!    coefficient of X**N.
!
!    Output, real POLY_COF2(0:N+1), the coefficients of the antiderivative
!    polynomial, in standard form.  The constant term is set to zero.
!
  implicit none
!
  integer n
!
  real poly_cof(0:n)
  real poly_cof2(0:n+1)
  integer i
!
!  Set the constant term.
!
  poly_cof2(0) = 0.0E+00
!
!  Integrate the polynomial.
!
  do i = 1, n+1
    poly_cof2(i) = poly_cof(i-1) / real ( i )
  end do

  return
end
subroutine rpoly_ant_val ( n, poly_cof, xval, yval )
!
!*******************************************************************************
!
!! RPOLY_ANT_VAL evaluates the antiderivative of a polynomial in standard form.
!
!
!  Discussion:
!
!    The constant term of the antiderivative is taken to be zero.
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the polynomial.
!
!    Input, real POLY_COF(0:N), the polynomial coefficients.  POLY_COF(0)
!    is the constant term, and POLY_COF(N) is the coefficient of
!    X**N.
!
!    Input, real XVAL, the point where the antiderivative is to be
!    evaluated.
!
!    Output, real YVAL, the value of the antiderivative of the polynomial
!    at XVAL.
!
  implicit none
!
  integer n
!
  integer i
  real poly_cof(0:n)
  real xval
  real yval
!
  yval = poly_cof(n) * xval / real ( n + 1 )

  do i = n-1, 0, -1
    yval = ( yval + poly_cof(i) / real ( i + 1 ) ) * xval
  end do

  return
end
subroutine rpoly_basis ( poly_cof, ntab, xtab )
!
!*******************************************************************************
!
!! RPOLY_BASIS computes all Lagrange basis polynomial in standard form.
!
!
!  Definition:
!
!    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
!    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
!    XTAB(J) for J not equal to I, and 1 when J is equal to I.
!
!    The Lagrange basis polynomials have the property that the interpolating
!    polynomial through a set of NTAB data points (XTAB,YTAB) may be
!    represented as
!
!      P(X) = Sum ( I = 1 to N ) YTAB(I) * L(I,NTAB,XTAB)(X)
!
!    Higher order interpolation at selected points may be accomplished
!    using repeated X values, and scaled derivative values.
!
!  Modified:
!
!    31 August 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real POLY_COF(NTAB,NTAB), the polynomial coefficients for the
!    I-th Lagrange basis polynomial are stored in column I.  POLY_COF(1,I)
!    is the constant term, and POLY_COF(1,NTAB) is the coefficient of
!    X**(NTAB-1).
!
!    Input, integer NTAB, the number of data points XTAB.
!
!    Input, real XTAB(NTAB), the X values upon which the Lagrange basis
!    polynomial is to be based.
!
  implicit none
!
  integer ntab
!
  integer i
  integer j
  real poly_cof(ntab,ntab)
  real xtab(ntab)
!
!  Initialize POLY_COF to the identity matrix.
!
  poly_cof(1:ntab,1:ntab) = 0.0E+00
  do i = 1, ntab
    poly_cof(i,i) = 1.0E+00
  end do
!
!  Compute the divided difference table for the IVAL-th Lagrange basis
!  polynomial.
!
  do i = 1, ntab
    call data_to_dif ( poly_cof(1,i), ntab, xtab, poly_cof(1,i) )
  end do
!
!  Convert the divided difference table coefficients to standard polynomial
!  coefficients.
!
  do i = 1, ntab
    call dif_to_rpoly ( poly_cof(1,i), ntab, xtab )
  end do

  return
end
subroutine rpoly_basis_1 ( poly_cof, ival, ntab, xtab )
!
!*******************************************************************************
!
!! RPOLY_BASIS_1 computes the I-th Lagrange basis polynomial in standard form.
!
!
!  Definition:
!
!    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
!    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
!    XTAB(J) for J not equal to I, and 1 when J is equal to I.
!
!    The Lagrange basis polynomials have the property that the interpolating
!    polynomial through a set of NTAB data points (XTAB,YTAB) may be
!    represented as
!
!      P(X) = Sum ( I = 1 to N ) YTAB(I) * L(I,NTAB,XTAB)(X)
!
!    Higher order interpolation at selected points may be accomplished
!    using repeated X values, and scaled derivative values.
!
!  Modified:
!
!    12 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real POLY_COF(0:NTAB-1), the polynomial coefficients for the
!    IVAL-th Lagrange basis polynomial.
!
!    Input, integer IVAL, the index of the desired Lagrange basis polynomial.
!    IVAL should be between 1 and NTAB.
!
!    Input, integer NTAB, the number of data points XTAB.
!
!    Input, real XTAB(NTAB), the X values upon which the Lagrange basis
!    polynomial is to be based.
!
  implicit none
!
  integer ntab
!
  integer i
  integer ival
  real poly_cof(0:ntab-1)
  real xtab(ntab)
!
!  Check IVAL.
!
  if ( ival < 1 .or. ival > ntab ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RPOLY_BASE_1 - Fatal error!'
    write ( *, '(a,i6)' ) '  IVAL must be between 1 and ', ntab
    write ( *, '(a,i6)' ) '  but your value is ', ival
    stop
  end if
!
!  Initialize POLY_COF to the IVAL-th column of the identity matrix.
!
  poly_cof(0:ntab-1) = 0.0E+00
  poly_cof(ival-1) = 1.0E+00
!
!  Compute the divided difference table for the IVAL-th Lagrange basis
!  polynomial.
!
  call data_to_dif ( poly_cof, ntab, xtab, poly_cof )
!
!  Convert the divided difference table coefficients to standard polynomial
!  coefficients.
!
  call dif_to_rpoly ( poly_cof, ntab, xtab )

  return
end
subroutine rpoly_degree ( na, a, degree )
!
!*******************************************************************************
!
!! RPOLY_DEGREE returns the degree of a polynomial.
!
!
!  Discussion:
!
!    The degree of a polynomial is the index of the highest power
!    of X with a nonzero coefficient.
!
!    The degree of a constant polynomial is 0.  The degree of the
!    zero polynomial is debatable, but this routine returns the
!    degree as 0.
!
!  Modified:
!
!    20 March 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NA, the dimension of A.
!
!    Input, real A(0:NA), the coefficients of the polynomials.
!
!    Output, integer DEGREE, the degree of A.
!
  implicit none
!
  integer na
!
  real a(0:na)
  integer degree
!
  degree = na

  do while ( degree > 0 )

    if ( a(degree) /= 0.0E+00 ) then
      return
    end if

    degree = degree - 1

  end do

  return
end
subroutine rpoly_der_cof ( n, poly_cof, poly_cof2 )
!
!*******************************************************************************
!
!! RPOLY_DER_COF computes the coefficients of the derivative of a polynomial.
!
!
!  Modified:
!
!    12 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the original polynomial.
!
!    Input, real POLY_COF(0:N), the coefficients of the polynomial to
!    be differentiated.  POLY_COF(0) is the constant term, and
!    POLY_COF(N) is the coefficient of X**N.
!
!    Output, real POLY_COF2(0:N-1), the coefficients of the derivative of
!    the polynomial.
!
  implicit none
!
  integer n
!
  real poly_cof(0:n)
  real poly_cof2(0:n-1)
  integer i
!
  do i = 0, n-1
    poly_cof2(i) = real ( i + 1 ) * poly_cof(i+1)
  end do

  return
end
subroutine rpoly_der_val ( n, poly_cof, xval, yval )
!
!*******************************************************************************
!
!! RPOLY_DER_VAL evaluates the derivative of a polynomial in standard form.
!
!
!  Definition:
!
!    A polynomial in standard form, with coefficients POLY_COF(*),
!    may be written:
!
!      P(X) = POLY_COF(0)
!           + POLY_COF(1) * X
!           ...
!           + POLY_COF(N) * X**N
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the polynomial.
!
!    Input, real POLY_COF(0:N), the polynomial coefficients.  POLY_COF(0)
!    is the constant term, and POLY_COF(N) is the coefficient of
!    X**N.
!
!    Input, real XVAL, a value where the derivative of the polynomial
!    is to be evaluated.
!
!    Output, real YVAL, the value of the derivative of the polynomial
!    at XVAL.
!
  implicit none
!
  integer n
!
  integer i
  real poly_cof(0:n)
  real xval
  real yval
!
  yval = real ( n ) * poly_cof(n)
  do i = n-1, 1, -1
    yval = yval * xval + real ( i ) * poly_cof(i)
  end do

  return
end
subroutine rpoly_print ( n, a, title )
!
!*******************************************************************************
!
!! RPOLY_PRINT prints out a polynomial.
!
!
!  Discussion:
!
!    The power sum form is:
!
!      p(x) = a(0) + a(1)*x + ... + a(n-1)*x**(n-1) + a(n)*x**(n)
!
!  Modified:
!
!    06 August 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the dimension of A.
!
!    Input, real A(0:N), the polynomial coefficients.
!    A(0) is the constant term and
!    A(N) is the coefficient of X**N.
!
!    Input, character ( len = * ) TITLE, an optional title.
!
  implicit none
!
  integer n
!
  real a(0:n)
  integer i
  real mag
  integer n2
  character plus_minus
  character ( len = * ) title
!
  if ( len_trim ( title ) > 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, * ) ' '

  call rpoly_degree ( n, a, n2 )

  if ( a(n2) < 0.0E+00 ) then
    plus_minus = '-'
  else
    plus_minus = ' '
  end if

  mag = abs ( a(n2) )

  if ( n2 >= 2 ) then
    write ( *, '( '' p(x) = '', a1, g14.6, '' * x ^ '', i3 )' ) &
      plus_minus, mag, n2
  else if ( n2 == 1 ) then
    write ( *, '( '' p(x) = '', a1, g14.6, '' * x'' )' ) plus_minus, mag
  else if ( n2 == 0 ) then
    write ( *, '( '' p(x) = '', a1, g14.6 )' ) plus_minus, mag
  end if

  do i = n2-1, 0, -1

    if ( a(i) < 0.0E+00 ) then
      plus_minus = '-'
    else
      plus_minus = '+'
    end if

    mag = abs ( a(i) )

    if ( mag /= 0.0E+00 ) then

      if ( i >= 2 ) then
        write ( *, ' ( ''        '', a1, g14.6, '' * x ^ '', i3 )' ) &
          plus_minus, mag, i
      else if ( i == 1 ) then
        write ( *, ' ( ''        '', a1, g14.6, '' * x'' )' ) plus_minus, mag
      else if ( i == 0 ) then
        write ( *, ' ( ''        '', a1, g14.6 )' ) plus_minus, mag
      end if
    end if

  end do

  return
end
subroutine rpoly_shift ( scale, shift, n, poly_cof )
!
!*******************************************************************************
!
!! RPOLY_SHIFT adjusts the coefficients of a polynomial for a new argument.
!
!
!  Discussion:
!
!
!    Assuming P(X) is a polynomial in the argument X, of the form:
!
!      P(X) =
!          C(N) * X**N
!        + ...
!        + C(1) * X
!        + C(0),
!
!    and that Z is related to X by the formula:
!
!      Z = SCALE * X + SHIFT
!
!    then this routine computes coefficients C for the polynomial Q(Z):
!
!      Q(Z) =
!          C(N) * Z**N
!        + ...
!        + C(1) * Z
!        + C(0)
!
!    so that:
!
!      Q(Z(X)) = P(X)
!
!  Example:
!
!    P(X) = 2 * X**2 - X + 6
!
!    Z = 2.0E+00 * X + 3.0E+00
!
!    Q(Z) = 0.5 *         Z**2 -  3.5 * Z + 12
!
!    Q(Z(X)) = 0.5 * ( 4.0E+00 * X**2 + 12.0E+00 * X +  9 )
!            - 3.5 * (               2.0E+00 * X +  3 )
!                                            + 12
!
!            = 2.0E+00         * X**2 -  1.0E+00 * X +  6
!
!            = P(X)
!
!  Reference:
!
!    Press, Flannery, Teukolsky, Vetterling,
!    Numerical Recipes: The Art of Scientific Computing,
!    Cambridge University Press.
!
!  Modified:
!
!    05 October 1999
!
!  Parameters:
!
!    Input, real SHIFT, SCALE, the shift and scale applied to X,
!    so that Z = SCALE * X + SHIFT.
!
!    Input, integer N, the number of coefficients.
!
!    Input/output, real POLY_COF(0:N).
!    On input, the coefficient array in terms of the X variable.
!    On output, the coefficient array in terms of the Z variable.
!
  implicit none
!
  integer n
!
  integer i
  integer j
  real poly_cof(0:n)
  real scale
  real shift
!
  do i = 1, n
    poly_cof(i:n) = poly_cof(i:n) / scale
  end do

  do i = 0, n-1
    do j = n-1, i, -1
      poly_cof(j) = poly_cof(j) - shift * poly_cof(j+1)
    end do
  end do

  return
end
subroutine rpoly_val ( n, poly_cof, xval, yval )
!
!*******************************************************************************
!
!! RPOLY_VAL evaluates a polynomial in standard form.
!
!
!  Definition:
!
!    A polynomial in standard form, with coefficients POLY_COF(*),
!    may be written:
!
!      P(X) = POLY_COF(0)
!           + POLY_COF(1) * X
!           ...
!           + POLY_COF(N) * X**N
!
!  Modified:
!
!    11 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the degree of the polynomial.
!
!    Input, real POLY_COF(0:N), the polynomial coefficients.  POLY_COF(0)
!    is the constant term, and POLY_COF(N) is the coefficient of
!    X**N.
!
!    Input, real XVAL, a value where the polynomial is to be evaluated.
!
!    Output, real YVAL, the value of the polynomial at XVAL.
!
  implicit none
!
  integer n
!
  integer i
  real poly_cof(0:n)
  real xval
  real yval
!
  yval = poly_cof(n)
  do i = n-1, 0, -1
    yval = yval * xval + poly_cof(i)
  end do

  return
end
function rvec_distinct ( n, x )
!
!*******************************************************************************
!
!! RVEC_DISTINCT is true if the entries in a real vector are distinct.
!
!
!  Modified:
!
!    16 September 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(N), the vector to be checked.
!
!    Output, logical RVEC_DISTINCT is .TRUE. if all N elements of X
!    are distinct.
!
  implicit none
!
  integer n
!
  integer i
  integer j
  logical rvec_distinct
  real x(n)
!
  rvec_distinct = .false.

  do i = 2, n
    do j = 1, i - 1
      if ( x(i) == x(j) ) then
        return
      end if
    end do
  end do

  rvec_distinct = .true.

  return
end
subroutine rvec_identity ( n, a )
!
!*******************************************************************************
!
!! RVEC_IDENTITY sets a real vector to the identity vector A(I)=I.
!
!
!  Modified:
!
!    09 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of elements of A.
!
!    Output, real A(N), the array to be initialized.
!
  implicit none
!
  integer n
!
  real a(n)
  integer i
!
  do i = 1, n
    a(i) = real ( i )
  end do

  return
end
subroutine rvec_print ( n, a, title )
!
!*******************************************************************************
!
!! RVEC_PRINT prints a real vector.
!
!
!  Modified:
!
!    16 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of components of the vector.
!
!    Input, real A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed first.
!    TITLE may be blank.
!
  implicit none
!
  integer n
!
  real a(n)
  integer i
  character ( len = * ) title
!
  if ( title /= ' ' ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(i6,g14.6)' ) i, a(i)
  end do

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

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

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

  return
end