# include <cmath>
# include <cstdlib>
# include <cstring>
# include <iomanip>
# include <iostream>

using namespace std;

# include "lagrange.hpp"

//****************************************************************************80

double *lagrange_basis_antideriv ( int npol, int ipol, double *xpol )

//****************************************************************************80
//
//  Purpose:
//
//    lagrange_basis_antideriv() returns the antiderivative of a Lagrange basis polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 August 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int NPOL, the number of abscissas.
//    NPOL must be at least 1.
//
//    int 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 LAGRANGE_BASIS_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
//
{
  double *pcof;
  double *pcof_antideriv;
//
//  Get the standard polynomial coefficients of the ipol-th Lagrange
//  basis polynomial.
//
  pcof = lagrange_basis_coef ( npol, ipol, xpol );
//
//  Get the standard polynomial coefficients of the antiderivative of
//  the ipol-th Lagrange basis polynomial.
//
  pcof_antideriv = r8poly_ant_coef ( npol-1, pcof );

  delete [] pcof;

  return pcof_antideriv;
}
//****************************************************************************80

double *lagrange_basis_coef ( int npol, int ipol, double *xpol )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    11 August 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int NPOL, the number of abscissas.
//    NPOL must be at least 1.
//
//    int IPOL, the index of the polynomial to evaluate.
//    IPOL must be between 1 and NPOL.
//
//    double XPOL(NPOL), the abscissas of the
//    Lagrange polynomials.  The entries in XPOL must be distinct.
//
//  Output:
//
//    double 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
//
{
  int i;
  int indx;
  int j;
  double *pcof;
//
//   Make sure IPOL is legal.
//
  if ( ipol < 0 || npol <= ipol )
  {
    cout << "\n";
    cout << "lagrange_basis_coef(): Fatal error!\n";
    cout << "  0 <= IPOL < NPOL is required.\n";
    exit ( 1 );
  }
//
//  Check that the abscissas are distinct.
//
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    cout << "\n";
    cout << "lagrange_basis_coef(): Fatal error!\n";
    cout << "  Two or more entries of XPOL are equal:\n";
    exit ( 1 );
  }

  pcof = new double[npol];

  pcof[0] = 1.0;
  for ( j = 1; j < npol; j++ )
  {
    pcof[j] = 0.0;
  }

  indx = 0;

  for ( i = 0; i < npol; i++ )
  {
    if ( i != ipol )
    {
      indx = indx + 1;

      for ( j = indx; 0 <= j; j-- )
      {
        pcof[j] = - xpol[i] * pcof[j] / ( xpol[ipol] - xpol[i] );

        if ( 0 < j )
        {
          pcof[j] = pcof[j] + pcof[j-1] / ( xpol[ipol] - xpol[i] );
        }
      }
    }
  }

  return pcof;
}
//****************************************************************************80

double lagrange_basis_deriv ( int npol, int ipol, double *xpol, double xval )

//****************************************************************************80
//
//  Purpose:
//
//    lagrange_basis_deriv(): derivative of the IPOL-th Lagrange basis polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    06 August 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int NPOL, the number of abscissas.
//    NPOL must be at least 1.
//
//    int 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.
//
{
  double dpdx;
  int i;
  int j;
  double p;
//
//  Make sure IPOL is legal.
//
  if ( ipol < 0 || npol <= ipol )
  {
    cout << "\n";
    cout << "lagrange_basis_deriv(): Fatal error!\n";
    cout << "  0 <= IPOL < NPOL is required.\n";
    exit ( 1 );
  }
//
//  Check that the abscissas are distinct.
//
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    cout << "\n";
    cout << "lagrange_basis_deriv(): Fatal error!\n";
    cout << "  Two or more entries of XPOL are equal:\n";
    exit ( 1 );
  }
//
//  Evaluate the derivative, which can be found by summing up the result
//  of differentiating one factor at a time, successively.
//
  dpdx = 0.0;

  for ( i = 0; i < npol; i++ )
  {
    if ( i != ipol )
    {
      p = 1.0;
      for ( j = 0; j < npol; j++ )
      {
        if ( j == i )
        {
          p = p                      / ( xpol[ipol] - xpol[j] );
        }
        else if ( j != ipol )
        {
          p = p * ( xval - xpol[j] ) / ( xpol[ipol] - xpol[j] );
        }
      }
      dpdx = dpdx + p;
    }
  }

  return dpdx;
}
//****************************************************************************80

double lagrange_basis_deriv2 ( int npol, int ipol, double *xpol, double xval )

//****************************************************************************80
//
//  Purpose:
//
//    lagrange_basis_deriv2(): second derivative of the IPOL-th Lagrange basis polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 August 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int NPOL, the number of abscissas.
//    NPOL must be at least 1.
//
//    int 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.
//
{
  double d2pdx2;
  int i;
  int j;
  int k;
  double p;
//
//   Make sure IPOL is legal.
//
  if ( ipol < 0 || npol <= ipol )
  {
    cout << "\n";
    cout << "lagrange_basis_deriv2(): Fatal error!\n";
    cout << "  0 <= IPOL < NPOL is required.\n";
    exit ( 1 );
  }
//
//  Check that the abscissas are distinct.
//
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    cout << "\n";
    cout << "lagrange_basis_deriv2(): Fatal error!\n";
    cout << "  Two or more entries of XPOL are equal:\n";
    exit ( 1 );
  }
//
//  Evaluate the second derivative, by summing up the result
//  of differentiating with respect to every pair of components,
//  omitting IPOL.
//
  d2pdx2 = 0.0;

  for ( i = 0; i < npol; i++ )
  {
    if ( i != ipol )
    {
      for ( j = 0; j < npol; j++ )
      {
        if ( j != ipol && j != i )
        {
          p = 1.0;
          for ( k = 0; k < npol; k++ )
          {
            if ( k != ipol )
            {
              if ( k == i || k == j )
              {
                p = p                      / ( xpol[ipol] - xpol[k] );
              }
              else
              {
                p = p * ( xval - xpol[k] ) / ( xpol[ipol] - xpol[k] );
              }
            }
          }
        }
      }
      d2pdx2 = d2pdx2 + p;
    }
  }

  return d2pdx2;
}
//****************************************************************************80

double lagrange_basis_value ( int npol, int ipol, double *xpol, double xval )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    11 August 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int NPOL, the number of abscissas.
//    NPOL must be at least 1.
//
//    int IPOL, the index of the polynomial to evaluate.
//    0 <= IPOL < NPOL.
//
//    double XPOL(NPOL), the abscissas of the Lagrange
//    polynomials.  The entries in XPOL must be distinct.
//
//    double XVAL, the point at which the IPOL-th
//    Lagrange polynomial is to be evaluated.
//
//  Output:
//
//    double lagrange_basis_val: the value of the IPOL-th Lagrange
//    polynomial at XVAL.
//
{
  int i;
  double pval;
//
//  Make sure IPOL is legal.
//
  if ( ipol < 0 || npol <= ipol )
  {
    cout << "\n";
    cout << "lagrange_basis_value(): Fatal error!\n";
    cout << "  0 <= IPOL < NPOL is required.\n";
    exit ( 1 );
  }
//
//  Check that the abscissas are distinct.
//
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    cout << "\n";
    cout << "lagrange_basis_value(): Fatal error!\n";
    cout << "  Two or more entries of XPOL are equal:\n";
    exit ( 1 );
  }
//
//  Evaluate the polynomial.
//
  pval = 1.0;

  for ( i = 0; i < npol; i++ )
  {
    if ( i != ipol )
    {
      pval = pval * ( xval - xpol[i] ) / ( xpol[ipol] - xpol[i] );
    }
  }

  return pval;
}
//****************************************************************************80

double *r8poly_ant_coef ( int n, double poly_cof[] )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    int N, the order of the polynomial.
//
//    double 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:
//
//    double POLY_COF2[0:N+1], the coefficients of
//    the antiderivative polynomial, in standard form.  The constant
//    term is set to zero.
//
{
  int i;
  double *poly_cof2;

  poly_cof2 = new double[n+2];
//
//  Set the constant term.
//
  poly_cof2[0] = 0.0;
//
//  Integrate the polynomial.
//
  for ( i = 1; i < n + 2; i++ )
  {
    poly_cof2[i] = poly_cof[i-1] / ( double ) ( i );
  }

  return poly_cof2;
}
//****************************************************************************80

void r8poly_print ( int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8poly_print() prints a polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 July 2015
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the dimension of A.
//
//    double A[N+1], the polynomial coefficients.
//    A(0) is the constant term and
//    A(N) is the coefficient of X^N.
//
//    string TITLE, a title.
//
{
  int i;
  double mag;
  char plus_minus;

  if ( 0 < title.length ( ) )
  {
    cout << "\n";
    cout << title << "\n";
  }
  cout << "\n";

  if ( n < 0 )
  {
    cout << "  p(x) = 0\n";
    return;
  }

  if ( a[n] < 0.0 )
  {
    plus_minus = '-';
  }
  else
  {
    plus_minus = ' ';
  }

  mag = fabs ( a[n] );

  if ( 2 <= n )
  {
    cout << "  p(x) = " << plus_minus
         << setw(14) << mag << " * x ^ " << n << "\n";
  }
  else if ( n == 1 )
  {
    cout << "  p(x) = " << plus_minus
         << setw(14) << mag << " * x\n";
  }
  else if ( n == 0 )
  {
    cout << "  p(x) = " << plus_minus
         << setw(14) << mag << "\n";
  }

  for ( i = n - 1; 0 <= i; i-- )
  {
    if ( a[i] < 0.0 )
    {
      plus_minus = '-';
    }
    else
    {
      plus_minus = '+';
    }

    mag = fabs ( a[i] );

    if ( mag != 0.0 )
    {
      if ( 2 <= i )
      {
        cout << "         " << plus_minus
             << setw(14) << mag << " * x ^ " << i << "\n";
      }
      else if ( i == 1 )
      {
        cout << "         " << plus_minus
             << setw(14) << mag << " * x\n";
      }
      else if ( i == 0 )
      {
        cout << "         " << plus_minus
             << setw(14) << mag << "\n";
      }
    }
  }

  return;
}
//****************************************************************************80

double r8poly_value ( int m, double c[], double x )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    int M, the degree of the polynomial.
//
//    double C[M+1], the coefficients of the polynomial.
//    A[0] is the constant term.
//
//    double X, the point at which the polynomial is to be evaluated.
//
//  Output:
//
//    double R8POLY_VALUE, the value of the polynomial at X.
//
{
  int i;
  double value;
  double xi;

  value = c[0];
  xi = 1.0;

  for ( i = 1; i <= m; i++ )
  {
    xi = xi * x;
    value = value + c[i] * xi;
  }

  return value;
}
//****************************************************************************80

bool r8vec_is_distinct ( int n, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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:
//
//    05 April 2018
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of entries in the vector.
//
//    double X[N], the vector to be checked.
//
//  Output:
//
//    bool R8VEC_IS_DISTINCT is true if all entries are distinct.
//
{
  int i;
  int j;
  bool value;

  value = true;

  for ( i = 1; i < n; i++ )
  {
    for ( j = 0; j < i; j++ )
    {
      if ( x[i] == x[j] )
      {
        value = false;
        break;
      }
    }
  }
  return value;
}
//****************************************************************************80

void r8vec_print ( int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    r8vec_print() prints an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of components of the vector.
//
//    double A[N], the vector to be printed.
//
//    string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8)  << i
         << ": " << setw(14) << a[i]  << "\n";
  }

  return;
}
//****************************************************************************80

int s_len_trim ( string s )

//****************************************************************************80
//
//  Purpose:
//
//    s_len_trim() returns the length of a string to the last nonblank.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 October 2014
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    string S, a string.
//
//  Output:
//
//    int S_LEN_TRIM, the length of the string to the last nonblank.
//    If S_LEN_TRIM is 0, then the string is entirely blank.
//
{
  int n;

  n = s.length ( );

  while ( 0 < n )
  {
    if ( s[n-1] != ' ' && s[n-1] != '\n' )
    {
      return n;
    }
    n = n - 1;
  }

  return n;
}