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

# include "barycentric_interp_1d.h"
# include "r8lib.h"

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

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

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

    LAGCHEBY1_INTERP_1D evaluates the Lagrange Chebyshev 1 interpolant.

  Discussion:

    The weight vector WD computed below is only valid if the data points
    XD are, as expected, the Chebyshev Type 1 points for [-1,+1], or a linearly 
    mapped version for [A,B].  The XD values may be computed by:

      xd = r8vec_cheby1space ( nd, a, b )

    for instance.

    Thanks to John Ferrier for pointing out that DENOM and NUMER needed
    to be initialized, 16 September 2013.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 September 2013

  Author:

    John Burkardt

  Reference:

    Jean-Paul Berrut, Lloyd Trefethen,
    Barycentric Lagrange Interpolation,
    SIAM Review,
    Volume 46, Number 3, September 2004, pages 501-517.

  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 LAGCHEBY1_INTERP_1D[NI], the interpolated values.
*/
{
  double *denom;
  int *exact;
  int i;
  int j;
  double *numer;
  double pi = 3.141592653589793;
  double t;
  double theta;
  double wd;
  double *yi;

  denom = ( double * ) malloc ( ni * sizeof ( double ) );
  exact = ( int * ) malloc ( ni * sizeof ( int ) );
  numer = ( double * ) malloc ( ni * sizeof ( double ) );
  yi = ( double * ) malloc ( ni * sizeof ( double ) );

  for ( i = 0; i < ni; i++ )
  {
    exact[i] = -1;
    denom[i] = 0.0;
    numer[i] = 0.0;
  }

  for ( j = 0; j < nd; j++ )
  {
    theta = ( double ) ( 2 * j - 1 ) * pi / ( double ) ( 2 * nd );
    wd = r8_mop ( j + 1 ) * sin ( theta );

    for ( i = 0; i < ni; i++ )
    {
      if ( xi[i] == xd[j] )
      {
        exact[i] = j;
        numer[i] = yd[j];
        denom[i] = 1.0;
      }

      if ( exact[i] == -1 )
      {
        t = wd / ( xi[i] - xd[j] );
        numer[i] = numer[i] + t * yd[j];
        denom[i] = denom[i] + t;
      }
    }
  }

  for ( i = 0; i < ni; i++ )
  {
    yi[i] = numer[i] / denom[i];
  }

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

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

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

    LAGCHEBY2_INTERP_1D evaluates the Lagrange Chebyshev 2 interpolant.

  Discussion:

    The weight vector WD computed below is only valid if the data points
    XD are, as expected, the Chebyshev Type 2 points for [-1,+1], or a linearly 
    mapped version for [A,B].  The XD values may be computed by:

      xd = r8vec_cheby2space ( nd, a, b )

    for instance.

    Thanks to John Ferrier for pointing out that DENOM and NUMER needed
    to be initialized, 16 September 2013.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 September 2013

  Author:

    John Burkardt

  Reference:

    Jean-Paul Berrut, Lloyd Trefethen,
    Barycentric Lagrange Interpolation,
    SIAM Review,
    Volume 46, Number 3, September 2004, pages 501-517.

  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 LAGCHEBY2_INTERP_1D[NI], the interpolated values.
*/
{
  double *denom;
  int *exact;
  int i;
  int j;
  double *numer;
  double t;
  double wd;
  double *yi;

  denom = ( double * ) malloc ( ni * sizeof ( double ) );
  exact = ( int * ) malloc ( ni * sizeof ( int ) );
  numer = ( double * ) malloc ( ni * sizeof ( double ) );
  yi = ( double * ) malloc ( ni * sizeof ( double ) );

  for ( i = 0; i < ni; i++ )
  {
    exact[i] = -1;
    denom[i] = 0.0;
    numer[i] = 0.0;
  }

  for ( j = 0; j < nd; j++ )
  {
    wd = r8_mop ( j );
    if ( j == 0 || j == nd - 1 )
    {
      wd = 0.5 * wd;
    }

    for ( i = 0; i < ni; i++ )
    {
      if ( xi[i] == xd[j] )
      {
        exact[i] = j;
        numer[i] = yd[j];
        denom[i] = 1.0;
      }

      if ( exact[i] == -1 )
      {
        t = wd / ( xi[i] - xd[j] );
        numer[i] = numer[i] + t * yd[j];
        denom[i] = denom[i] + t;
      }
    }
  }

  for ( i = 0; i < ni; i++ )
  {
    yi[i] = numer[i] / denom[i];
  }

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

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

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

    LAGEVEN_VALUE_1D evaluates the Lagrange evenly-spaced interpolant.

  Discussion:

    The weight vector WD computed below is only valid if the data points
    XD are, as expected, evenly spaced in an interval [A,B] with
    spacing (B-A)/N.  The XD values might be computed by:

      xd[i] = ( ( 2 * nd - 2 * i + 1 ) * a 
              + (          2 * i - 1 ) * b ) 
              / ( 2 * nd             )

    for instance.

    Thanks to John Ferrier for pointing out that DENOM and NUMER needed
    to be initialized, 16 September 2013.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    16 September 2013

  Author:

    John Burkardt

  Reference:

    Jean-Paul Berrut, Lloyd Trefethen,
    Barycentric Lagrange Interpolation,
    SIAM Review,
    Volume 46, Number 3, September 2004, pages 501-517.

  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 LAGEVEN_INTERP_1D[NI], the interpolated values.
*/
{
  double *denom;
  int *exact;
  int i;
  int j;
  double *numer;
  double t;
  double wd;
  double *yi;

  denom = ( double * ) malloc ( ni * sizeof ( double ) );
  exact = ( int * ) malloc ( ni * sizeof ( int ) );
  numer = ( double * ) malloc ( ni * sizeof ( double ) );
  yi = ( double * ) malloc ( ni * sizeof ( double ) );
  
  for ( i = 0; i < ni; i++ )
  {
    exact[i] = -1;
    denom[i] = 0.0;
    numer[i] = 0.0;
  }

  for ( j = 0; j < nd; j++ )
  {
    wd = r8_mop ( j ) * r8_choose ( nd, j );

    for ( i = 0; i < ni; i++ )
    {
      if ( xi[i] == xd[j] )
      {
        exact[i] = j;
        numer[i] = yd[j];
        denom[i] = 1.0;
      }

      if ( exact[i] == -1 )
      {
        t = wd / ( xi[i] - xd[j] );
        numer[i] = numer[i] + t * yd[j];
        denom[i] = denom[i] + t;
      }
    }
  }

  for ( i = 0; i < ni; i++ )
  {
    yi[i] = numer[i] / denom[i];
  }

  return yi;
}