subroutine lagrange_basis_antideriv ( npol, ipol, xpol, pcof_antideriv )

!*****************************************************************************80
!
!! lagrange_basis_antideriv() returns the antiderivative of a Lagrange basis polynomial.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    05 August 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer NPOL, the number of abscissas.
!    NPOL must be at least 1.
!
!    integer IPOL, the index of the polynomial to evaluate.
!    IPOL must be between 1 and NPOL.
!
!    real XPOL(NPOL), the abscissas of the
!    Lagrange polynomials.  The entries in XPOL must be distinct.
!
!  Output:
!
!    real PCOF_ANTIDERIV(1:NPOL+1), the standard polynomial
!    coefficients of the antiderivative of the IPOL-th Lagrange polynomial:
!      L(IPOL)(X) = SUM ( 0 <= I <= NPOL ) PCOF_ANTIDERIV(I+1) * X^I
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer npol

  integer ipol
  real ( kind = rk8 ) pcof(npol)
  real ( kind = rk8 ) pcof_antideriv(1:npol+1)
  real ( kind = rk8 ) xpol(npol)
!
!  Get the standard polynomial coefficients of the ipol-th Lagrange
!  basis polynomial.
!
  call lagrange_basis_coef ( npol, ipol, xpol, pcof )
!
!  Get the standard polynomial coefficients of the antiderivative of
!  the ipol-th Lagrange basis polynomial.
!
  call r8poly_ant_coef ( npol-1, pcof, pcof_antideriv )

  return
end
subroutine lagrange_basis_coef ( npol, ipol, xpol, pcof )

!*****************************************************************************80
!
!! lagrange_basis_coef() returns the coefficients of a Lagrange polynomial.
!
!  Discussion:
!
!    Given distinct abscissas XPOL(1:NPOL), the IPOL-th Lagrange
!    polynomial L(IPOL)(X) is defined as the polynomial of degree
!    NPOL - 1 which is 1 at XPOL(IPOL) and 0 at the NPOL - 1 other
!    abscissas.
!
!    A formal representation is:
!
!      L(IPOL)(X) = Product ( 1 <= I <= NPOL, I /= IPOL )
!       ( X - X(I) ) / ( X(IPOL) - X(I) )
!
!    However sometimes it is desirable to be able to write down
!    the standard polynomial coefficients of L(IPOL)(X).
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 1999
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer NPOL, the number of abscissas.
!    NPOL must be at least 1.
!
!    integer IPOL, the index of the polynomial to evaluate.
!    IPOL must be between 1 and NPOL.
!
!    real ( kind = rk8 ) XPOL(NPOL), the abscissas of the
!    Lagrange polynomials.  The entries in XPOL must be distinct.
!
!  Output:
!
!    real ( kind = rk8 ) PCOF(0:NPOL-1), the standard polynomial
!    coefficients of the IPOL-th Lagrange polynomial:
!      L(IPOL)(X) = SUM ( 0 <= I <= NPOL-1 ) PCOF(I) * X^I
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer npol

  integer i
  integer indx
  integer ipol
  integer j
  real ( kind = rk8 ) pcof(0:npol-1)
  logical r8vec_is_distinct
  real ( kind = rk8 ) xpol(npol)
!
!  Make sure IPOL is legal.
!
  if ( ipol < 1 .or. npol < ipol ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'lagrange_basis_coef(): Fatal error!'
    write ( *, '(a)' ) '  1 <= IPOL <= NPOL is required.'
    write ( *, '(a,i8)' ) '  IPOL = ', ipol
    write ( *, '(a,i8)' ) '  NPOL = ', npol
    stop 1
  end if
!
!  Check that the abscissas are distinct.
!
  if ( .not. r8vec_is_distinct ( npol, xpol ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'lagrange_basis_coef(): Fatal error!'
    write ( *, '(a)' ) '  Two or more entries of XPOL are equal:'
    stop 1
  end if

  pcof(0) = 1.0D+00
  pcof(1:npol-1) = 0.0D+00

  indx = 0

  do i = 1, npol

    if ( i /= ipol ) then

      indx = indx + 1

      do j = indx, 0, -1

        pcof(j) = -xpol(i) * pcof(j) / ( xpol(ipol) - xpol(i) )

        if ( 0 < j ) then
          pcof(j) = pcof(j) + pcof(j-1) / ( xpol(ipol) - xpol(i) )
        end if

      end do

    end if

  end do

  return
end
subroutine lagrange_basis_deriv ( npol, ipol, xpol, xval, dpdx )

!*****************************************************************************80
!
!! lagrange_basis_deriv(): derivative of the IPOL-th Lagrange basis polynomial.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    05 August 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer NPOL, the number of abscissas.
!    NPOL must be at least 1.
!
!    integer IPOL, the index of the polynomial to evaluate.
!    IPOL must be between 1 and NPOL.
!
!    real XPOL(NPOL), the abscissas of the Lagrange
!    polynomials.  The entries in XPOL must be distinct.
!
!    real XVAL, the evaluation point.
!
!  Output:
!
!    real DPDX, the derivative of the IPOL-th 
!    Lagrange polynomial at XVAL.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00)

  integer npol

  real ( kind = rk8 ) dpdx
  integer i
  integer ipol
  integer j
  real ( kind = rk8 ) p
  logical r8vec_is_distinct
  real ( kind = rk8 ) xpol(npol)
  real ( kind = rk8 ) xval
!
!  Make sure IPOL is legal.
!
  if ( ipol < 1 .or. npol < ipol ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'lagrange_basis_deriv(): Fatal error!'
    write ( *, '(a)' ) '  1 <= IPOL <= NPOL is required.'
    stop ( 1 )
  end if
!
!  Check that the abscissas are distinct.
!
  if ( .not. r8vec_is_distinct ( npol, xpol ) ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'lagrange_basis_deriv(): Fatal error!'
    write ( *, '(a)' ) '  Two or more entries of XPOL are equal:'
    stop ( 1 )
  end if
!
!  Evaluate the derivative, which can be found by summing up the result
!  of differentiating one factor at a time, successively.
!
  dpdx = 0.0D+00

  do i = 1, npol

    if ( i /= ipol ) then

      p = 1.0D+00
      do j = 1, npol

        if ( j == i ) then
          p = p                      / ( xpol(ipol) - xpol(j) )
        else if ( j /= ipol ) then
          p = p * ( xval - xpol(j) ) / ( xpol(ipol) - xpol(j) )
        end if

      end do

      dpdx = dpdx + p

    end if

  end do

  return
end
subroutine lagrange_basis_deriv2 ( npol, ipol, xpol, xval, d2pdx2 )

!*****************************************************************************80
!
!! lagrange_basis_deriv2(): second derivative of the IPOL-th Lagrange basis polynomial.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    05 August 2025
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer NPOL, the number of abscissas.
!    NPOL must be at least 1.
!
!    integer IPOL, the index of the polynomial to evaluate.
!    IPOL must be between 1 and NPOL.
!
!    real XPOL(NPOL), the abscissas of the Lagrange
!    polynomials.  The entries in XPOL must be distinct.
!
!    real XVAL, the evaluation point.
!
!  Output:
!
!    real D2PDX2, the second derivative of the IPOL-th 
!    Lagrange polynomial at XVAL.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00)

  integer npol

  real ( kind = rk8 ) d2pdx2
  integer i
  integer ipol
  integer j
  integer k
  real ( kind = rk8 ) p
  logical r8vec_is_distinct
  real ( kind = rk8 ) xpol(npol)
  real ( kind = rk8 ) xval
!
!  Make sure IPOL is legal.
!
  if ( ipol < 1 .or. npol < ipol ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'lagrange_basis_deriv2(): Fatal error!'
    write ( *, '(a)' ) '  1 <= IPOL <= NPOL is required.'
    stop ( 1 )
  end if
!
!  Check that the abscissas are distinct.
!
  if ( .not. r8vec_is_distinct ( npol, xpol ) ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'lagrange_basis_deriv2(): Fatal error!'
    write ( *, '(a)' ) '  Two or more entries of XPOL are equal:'
    stop ( 1 )
  end if
!
!  Evaluate the second derivative, by summing up the result
!  of differentiating with respect to every pair of components,
!  omitting IPOL.
!
  d2pdx2 = 0.0D+00

  do i = 1, npol

    if ( i /= ipol ) then

      do j = 1, npol

        if ( j /= ipol .and. j /= i ) then

          p = 1.0D+00
          do k = 1, npol

            if ( k /= ipol ) then

              if ( k == i .or. k == j ) then
                p = p                      / ( xpol(ipol) - xpol(k) )
              else
                p = p * ( xval - xpol(k) ) / ( xpol(ipol) - xpol(k) )
              end if

            end if

          end do

        end if

      end do

      d2pdx2 = d2pdx2 + p

    end if

  end do

  return
end
subroutine lagrange_basis_value ( npol, ipol, xpol, xval, pval, dpdx )

!*****************************************************************************80
!
!! lagrange_basis_value() evaluates the IPOL-th Lagrange polynomial.
!
!  Discussion:
!
!    Given NPOL distinct abscissas, XPOL(1:NPOL), the IPOL-th Lagrange
!    polynomial L(IPOL)(X) is defined as the polynomial of degree
!    NPOL - 1 which is 1 at XPOL(IPOL) and 0 at the NPOL - 1 other
!    abscissas.
!
!    A formal representation is:
!
!      L(IPOL)(X) = Product ( 1 <= I <= NPOL, I /= IPOL )
!       ( X - X(I) ) / ( X(IPOL) - X(I) )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 1999
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer NPOL, the number of abscissas.
!    NPOL must be at least 1.
!
!    integer IPOL, the index of the polynomial to evaluate.
!    IPOL must be between 1 and NPOL.
!
!    real ( kind = rk8 ) XPOL(NPOL), the abscissas of the Lagrange
!    polynomials.  The entries in XPOL must be distinct.
!
!    real ( kind = rk8 ) XVAL, the point at which the IPOL-th
!    Lagrange polynomial is to be evaluated.
!
!  Output:
!
!    real ( kind = rk8 ) PVAL, the value of the IPOL-th Lagrange
!    polynomial at XVAL.
!
!    real ( kind = rk8 ) DPDX, the derivative of the IPOL-th
!    Lagrange polynomial at XVAL.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer npol

  real ( kind = rk8 ) dpdx
  integer i
  integer ipol
  integer j
  real ( kind = rk8 ) p2
  real ( kind = rk8 ) pval
  logical r8vec_is_distinct
  real ( kind = rk8 ) xpol(npol)
  real ( kind = rk8 ) xval
!
!  Make sure IPOL is legal.
!
  if ( ipol < 1 .or. npol < ipol ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'lagrange_basis_value(): Fatal error!'
    write ( *, '(a)' ) '  1 <= IPOL <= NPOL is required.'
    write ( *, '(a,i8)' ) '  IPOL = ', ipol
    stop 1
  end if
!
!  Check that the abscissas are distinct.
!
  if ( .not. r8vec_is_distinct ( npol, xpol ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'lagrange_basis_value(): Fatal error!'
    write ( *, '(a)' ) '  Two or more entries of XPOL are equal:'
    stop 1
  end if
!
!  Evaluate the polynomial.
!
  pval = 1.0D+00

  do i = 1, npol

    if ( i /= ipol ) then

      pval = pval * ( xval - xpol(i) ) / ( xpol(ipol) - xpol(i) )

    end if

  end do
!
!  Evaluate the derivative, which can be found by summing up the result
!  of differentiating one factor at a time, successively.
!
  dpdx = 0.0D+00

  do i = 1, npol

    if ( i /= ipol ) then

      p2 = 1.0D+00
      do j = 1, npol

        if ( j == i ) then
          p2 = p2                      / ( xpol(ipol) - xpol(j) )
        else if ( j /= ipol ) then
          p2 = p2 * ( xval - xpol(j) ) / ( xpol(ipol) - xpol(j) )
        end if

      end do

      dpdx = dpdx + p2

    end if

  end do

  return
end
subroutine r8poly_ant_coef ( n, poly_cof, poly_cof2 )

!*****************************************************************************80
!
!! r8poly_ant_coef() integrates a polynomial in standard form.
!
!  Discussion:
!
!    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.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    15 October 2019
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the order of the polynomial.
!
!    real ( kind = rk8 ) POLY_COF(1:N+1), the polynomial coefficients.
!    POLY_COF(1) is the constant term, and POLY_COF(N+1) is the
!    coefficient of X^(N).
!
!  Output:
!
!    real ( kind = rk8 ) POLY_COF2(1:N+2), the coefficients of
!    the antiderivative polynomial, in standard form.  The constant
!    term is set to zero.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  integer i
  real ( kind = rk8 ) poly_cof(n+1)
  real ( kind = rk8 ) poly_cof2(n+2)
!
!  Set the constant term.
!
  poly_cof2(1) = 0.0D+00
!
!  Integrate the polynomial.
!
  do i = 2, n + 2
    poly_cof2(i) = poly_cof(i-1) / real ( i - 1, kind = rk8 )
  end do

  return
end
subroutine r8poly_print ( n, a, title )

!*****************************************************************************80
!
!! r8poly_print() prints 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)
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    15 July 2015
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the dimension of A.
!
!    real ( kind = rk8 ) A(0:N), the polynomial coefficients.
!    A(0) is the constant term and
!    A(N) is the coefficient of X^N.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(0:n)
  integer i
  real ( kind = rk8 ) mag
  character plus_minus
  character ( len = * ) title

  if ( 0 < len_trim ( title ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
  end if

  write ( *, '(a)' ) ' '

  if ( n < 0 ) then
    write ( *, '( ''  p(x) = 0'' )' )
    return
  end if

  if ( a(n) < 0.0D+00 ) then
    plus_minus = '-'
  else
    plus_minus = ' '
  end if

  mag = abs ( a(n) )

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

  do i = n - 1, 0, -1

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

    mag = abs ( a(i) )

    if ( mag /= 0.0D+00 ) then

      if ( 2 <= i ) 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
function r8poly_value ( m, c, x )

!*****************************************************************************80
!
!! r8poly_value() evaluates a polynomial using a naive method.
!
!  Discussion:
!
!    The polynomial
!
!      p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m
!
!    is to be evaluated at the value X.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    22 August 2017
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer M, the degree.
!
!    real ( kind = rk8 ) C(0:M), the polynomial coefficients.
!    C(I) is the coefficient of X^I.
!
!    real ( kind = rk8 ) X, the evaluation point.
!
!  Output:
!
!    real ( kind = rk8 ) R8POLY_VALUE, the polynomial value.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer m

  real ( kind = rk8 ) c(0:m)
  integer i
  real ( kind = rk8 ) r8poly_value
  real ( kind = rk8 ) value
  real ( kind = rk8 ) x
  real ( kind = rk8 ) xi

  value = c(0)
  xi = 1.0D+00
  do i = 1, m
    xi = xi * x
    value = value + c(i) * xi
  end do

  r8poly_value = value

  return
end
function r8vec_is_distinct ( n, a )

!*****************************************************************************80
!
!! r8vec_is_distinct() is true if the entries in an R8VEC are distinct.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    31 March 2018
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of entries in the vector.
!
!    real ( kind = rk8 ) A(N), the vector to be checked.
!
!  Output:
!
!    logical R8VEC_IS_DISTINCT is TRUE if the entries 
!    are distinct.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n)
  integer i
  integer j
  logical r8vec_is_distinct
  logical value

  value = .true.

  do i = 2, n
    do j = 1, i - 1
      if ( a(i) == a(j) ) then
        value = .false.
        exit
      end if
    end do
  end do

  r8vec_is_distinct = value

  return
end
subroutine r8vec_print ( n, a, title )

!*****************************************************************************80
!
!! r8vec_print() prints an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 September 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of components of the vector.
!
!    real ( kind = rk8 ) A(N), the vector to be printed.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a(n)
  integer i
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )
  write ( *, '(a)' ) ' '

  do i = 1, n
    write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i)
  end do

  return
end