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

# include "lagrange_interp_1d.h"
# include "r8lib.h"

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

double *lagrange_basis_1d ( int nd, double xd[], int ni, double xi[] )

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

    LAGRANGE_BASIS_1D evaluates the Lagrange basis polynomials.

  Discussion:

    Given ND distinct abscissas, XD(1:ND),
    the I-th Lagrange basis polynomial LB(I)(T) is defined as the polynomial of
    degree ND - 1 which is 1 at XD(I) and 0 at the ND - 1
    other abscissas.

    A formal representation is:

      LB(I)(T) = Product ( 1 <= J <= ND, I /= J )
       ( T - T(J) ) / ( T(I) - T(J) )

    This routine accepts a set of NI values at which all the Lagrange
    basis polynomials should be evaluated.

    Given data values YD at each of the abscissas, the value of the
    Lagrange interpolating polynomial at each of the interpolation points
    is then simple to compute by matrix multiplication:

      YI(1:NI) = LB(1:NI,1:ND) * YD(1:ND)

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    13 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int ND, the number of data points.
    ND must be at least 1.

    Input, double XD[ND], the data points.

    Input, int NI, the number of interpolation points.

    Input, double XI[NI], the interpolation points.

    Output, double LAGRANGE_BASIS[NI*ND], the values
    of the Lagrange basis polynomials at the interpolation points.
*/
{
  int i;
  int j;
  int k;
  double *lb;
/*
  Evaluate the polynomial.
*/
  lb = ( double * ) malloc ( ni * nd * sizeof ( double ) );

  for ( j = 0; j < nd; j++ )
  {
    for ( i = 0; i < ni; i++ )
    {
      lb[i+j*ni] = 1.0;
    }
  }

  for ( i = 0; i < nd; i++ )
  {
    for ( j = 0; j < nd; j++ )
    {
      if ( j != i )
      {
        for ( k = 0; k < ni; k++ )
        {
          lb[k+i*ni] = lb[k+i*ni] * ( xi[k] - xd[j] ) / ( xd[i] - xd[j] );
        }
      }
    }
  }

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

double *lagrange_value_1d ( int nd, double xd[], double yd[], int ni, 
  double xi[] )

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

    LAGRANGE_VALUE_1D evaluates the Lagrange interpolant.

  Discussion:

    The Lagrange interpolant L(ND,XD,YD)(X) is the unique polynomial of
    degree ND-1 which interpolates the points (XD(I),YD(I)) for I = 1
    to ND.

    The Lagrange interpolant can be constructed from the Lagrange basis
    polynomials.  Given ND distinct abscissas, XD(1:ND), the I-th Lagrange 
    basis polynomial LB(ND,XD,I)(X) is defined as the polynomial of degree 
    ND - 1 which is 1 at  XD(I) and 0 at the ND - 1 other abscissas.

    Given data values YD at each of the abscissas, the value of the
    Lagrange interpolant may be written as

      L(ND,XD,YD)(X) = sum ( 1 <= I <= ND ) LB(ND,XD,I)(X) * YD(I)

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    13 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int ND, the number of data points.
    ND must be at least 1.

    Input, double XD[ND], the data points.

    Input, double YD[ND], the data values.

    Input, int NI, the number of interpolation points.

    Input, double XI[NI], the interpolation points.

    Output, double LAGRANGE_VALUE_1D[NI], the interpolated values.
*/
{
  double *lb;
  double *yi;

  lb = lagrange_basis_1d ( nd, xd, ni, xi );

  yi = r8mat_mv_new ( ni, nd, lb, yd );

  free ( lb );

  return yi;
}