# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>

# include "lagrange.h"

/******************************************************************************/

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

/******************************************************************************/
/*
  Purpose:

    lagrange_basis_antideriv() returns the antiderivative of a 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.

  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 );

  free ( pcof );

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

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

/******************************************************************************/
/*
  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:

    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.

    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 )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_coef(): Fatal error!\n" );
    printf ( "  0 <= IPOL < NPOL is required.\n" );
    exit ( 1 );
  }
/*
  Check that the abscissas are distinct.
*/
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_coef(): Fatal error!\n" );
    printf ( "  Two or more entries of XPOL are equal:\n" );
    exit ( 1 );
  }

  pcof = ( double * ) malloc ( npol * sizeof ( double ) );

  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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 )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_deriv(): Fatal error!\n" );
    printf ( "  0 <= IPOL < NPOL is required.\n" );
    exit ( 1 );
  }
/*
  Check that the abscissas are distinct.
*/
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_deriv(): Fatal error!\n" );
    printf ( "  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  Purpose:

    lagrange_basis_deriv2(): second 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 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 )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_deriv2(): Fatal error!\n" );
    printf ( "  0 <= IPOL < NPOL is required.\n" );
    exit ( 1 );
  }
/*
  Check that the abscissas are distinct.
*/
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_deriv2(): Fatal error!\n" );
    printf ( "  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    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.
    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 )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_value(): Fatal error!\n" );
    printf ( "  0 <= IPOL < NPOL is required.\n" );
    exit ( 1 );
  }
/*
  Check that the abscissas are distinct.
*/
  if ( ! r8vec_is_distinct ( npol, xpol ) )
  {
    printf ( "\n" );
    printf ( "lagrange_basis_value(): Fatal error!\n" );
    printf ( "  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;
}
/******************************************************************************/

double *r8poly_ant_coef ( int n, double *poly_cof )

/******************************************************************************/
/*
  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 = ( double * ) malloc ( ( n + 2 ) * sizeof ( double ) );
/*
  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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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.

    char *TITLE, a title.
*/
{
  int i;
  double mag;
  char plus_minus;

  if ( 0 < s_len_trim ( title ) )
  {
    printf ( "\n" );
    printf ( "%s\n", title );
  }
  printf ( "\n" );

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

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

  mag = fabs ( a[n] );

  if ( 2 <= n )
  {
    printf ( "  p(x) = %c%g * x^%d\n", plus_minus, mag, n );
  }
  else if ( n == 1 )
  {
    printf ( "  p(x) = %c%g * x\n", plus_minus, mag );
  }
  else if ( n == 0 )
  {
    printf ( "  p(x) = %c%g\n", plus_minus, mag );
  }

  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 )
      {
        printf ( "         %c%g * x^%d\n", plus_minus, mag, i );
      }
      else if ( i == 1 )
      {
        printf ( "         %c%g * x\n", plus_minus, mag );
      }
      else if ( i == 0 )
      {
        printf ( "         %c%g\n", plus_minus, mag );
      }
    }
  }
  return;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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[0:M], the coefficients of the polynomial.
    C[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;
}
/******************************************************************************/

bool r8vec_is_distinct ( int n, double x[] )

/******************************************************************************/
/*
  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:

    31 March 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 N elements of X
    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;
}
/******************************************************************************/

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

/******************************************************************************/
/*
  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:

    08 April 2009

  Author:

    John Burkardt

  Input:

    int N, the number of components of the vector.

    double A[N], the vector to be printed.

    char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );
  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %8d: %14g\n", i, a[i] );
  }

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

int s_len_trim ( char *s )

/******************************************************************************/
/*
  Purpose:

    s_len_trim() returns the length of a string to the last nonblank.

  Discussion:

    It turns out that I also want to ignore the '\n' character!

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    05 October 2014

  Author:

    John Burkardt

  Input:

    char *S, a pointer to 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;
  char *t;

  n = strlen ( s );
  t = s + strlen ( s ) - 1;

  while ( 0 < n )
  {
    if ( *t != ' ' && *t != '\n' )
    {
      return n;
    }
    t--;
    n--;
  }

  return n;
}