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

# include "barycentric_interp_1d.h"
# include "test_interp_1d.h"

int main ( );
void lagcheby1_interp_1d_test ( int prob, int nd );
void lagcheby2_interp_1d_test ( int prob, int nd );
void lageven_interp_1d_test ( int prob, int nd );
double r8_choose ( int n, int k );
double r8_mop ( int i );
int r8_nint ( double x );
double *r8vec_cheby1space_new ( int n, double a, double b );
double *r8vec_cheby2space_new ( int n, double a, double b );
double *r8vec_copy_new ( int n, double a1[] );
double *r8vec_midspace_new ( int n, double a, double b );
double r8vec_norm_affine ( int n, double v0[], double v1[] );
void r8vec2_print ( int n, double a1[], double a2[], char *title );
void timestamp ( );

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

int main ( )

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

    barycentric_interp_1d_test() tests barycentric_interp_1d().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 October 2012

  Author:

    John Burkardt
*/
{
  int i;
  int nd;
  int nd_test[6] = { 4, 8, 16, 32, 64, 1000 };
  int nd_test_num = 6;
  int prob;
  int prob_num;

  timestamp ( );
  printf ( "\n" );
  printf ( "barycentric_interp_1d_test():\n" );
  printf ( "  C version\n" );
  printf ( "  Test barycentric_interp_1d().\n" );
  printf ( "  These tests need the test_interp_1d() library.\n" );

  prob_num = p00_prob_num ( );

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    for ( i = 0; i < nd_test_num; i++ )
    {
      nd = nd_test[i];
      lagcheby1_interp_1d_test ( prob, nd );
    }
  }

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    for ( i = 0; i < nd_test_num; i++ )
    {
      nd = nd_test[i];
      lagcheby2_interp_1d_test ( prob, nd );
    }
  }

  for ( prob = 1; prob <= prob_num; prob++ )
  {
    for ( i = 0; i < nd_test_num; i++ )
    {
      nd = nd_test[i];
      lageven_interp_1d_test ( prob, nd );
    }
  }
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "barycentric_interp_1d_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void lagcheby1_interp_1d_test ( int prob, int nd )

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

    lagcheby1_interp_1d_test() tests lagcheby1_interp_1d().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    30 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int PROB, the problem index.

    Input, int ND, the number of data points to use.
*/
{
  double a;
  double b;
  double int_error;
  int ni;
  double *xd;
  double *xi;
  double *yd;
  double *yi;

  printf ( "\n" );
  printf ( "lagcheby1_interp_1d_test():\n" );
  printf ( "  lagcheby1_interp_1d() uses Chebyshev Type 1 spacing for data points.\n" );
  printf ( "  Interpolate data from test_interp_1d() problem #%d\n", prob );
  printf ( "  Number of data points = %d\n", nd );
/*
  Define the data.
*/
  a =  0.0;
  b = +1.0;
  xd = r8vec_cheby1space_new ( nd, a, b );
  yd = p00_f ( prob, nd, xd );

  if ( nd < 10 )
  {
    r8vec2_print ( nd, xd, yd, "  Data array:" );
  }
/*
  #1:  Does the interpolant match the function at the interpolation points?
*/
  ni = nd;
  xi = r8vec_copy_new ( ni, xd );
  yi = lagcheby1_interp_1d ( nd, xd, yd, ni, xi );

  int_error = r8vec_norm_affine ( ni, yi, yd ) / ( double ) ( ni );

  printf ( "\n" );
  printf ( "  L2 interpolation error averaged per interpolant node = %g\n", int_error );

  free ( xd );
  free ( xi );
  free ( yd );
  free ( yi );

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

void lagcheby2_interp_1d_test ( int prob, int nd )

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

    lagcheby2_interp_1d_test() tests lagcheby2_interp_1d().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    30 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int PROB, the problem index.

    Input, int ND, the number of data points to use.
*/
{
  double a;
  double b;
  double int_error;
  int ni;
  double *xd;
  double *xi;
  double *yd;
  double *yi;

  printf ( "\n" );
  printf ( "lagcheby2_interp_1d_test():\n" );
  printf ( "  lagcheby2_interp_1d() uses Chebyshev Type 2 spacing for data points.\n" );
  printf ( "  Interpolate data from test_interp_1d() problem #%d\n", prob );
  printf ( "  Number of data points = %d\n", nd );
/*
  Define the data.
*/
  a =  0.0;
  b = +1.0;
  xd = r8vec_cheby2space_new ( nd, a, b );
  yd = p00_f ( prob, nd, xd );

  if ( nd < 10 )
  {
    r8vec2_print ( nd, xd, yd, "  Data array:" );
  }
/*
  #1:  Does the interpolant match the function at the interpolation points?
*/
  ni = nd;
  xi = r8vec_copy_new ( ni, xd );
  yi = lagcheby2_interp_1d ( nd, xd, yd, ni, xi );

  int_error = r8vec_norm_affine ( ni, yi, yd ) / ( double ) ( ni );

  printf ( "\n" );
  printf ( "  L2 interpolation error averaged per interpolant node = %g\n", int_error );

  free ( xd );
  free ( xi );
  free ( yd );
  free ( yi );

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

void lageven_interp_1d_test ( int prob, int nd )

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

    lageven_interp_1d_test() tests lageven_interp_1d().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    30 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int PROB, the problem index.

    Input, int ND, the number of data points to use.
*/
{
  double a;
  double b;
  double int_error;
  int ni;
  double *xd;
  double *xi;
  double *yd;
  double *yi;

  printf ( "\n" );
  printf ( "lageven_interp_1d_test():\n" );
  printf ( "  lageven_interp_1d() uses even spacing for data points.\n" );
  printf ( "  Interpolate data from test_interp_1d() problem #%d\n", prob );
  printf ( "  Number of data points = %d\n", nd );
/*
  Define the data.
*/
  a =  0.0;
  b = +1.0;
  xd = r8vec_midspace_new ( nd, a, b );
  yd = p00_f ( prob, nd, xd );

  if ( nd < 10 )
  {
    r8vec2_print ( nd, xd, yd, "  Data array:" );
  }
/*
  #1:  Does the interpolant match the function at the interpolation points?
*/
  ni = nd;
  xi = r8vec_copy_new ( ni, xd );
  yi = lageven_interp_1d ( nd, xd, yd, ni, xi );

  int_error = r8vec_norm_affine ( ni, yi, yd ) / ( double ) ( ni );

  printf ( "\n" );
  printf ( "  L2 interpolation error averaged per interpolant node = %g\n", int_error );

  free ( xd );
  free ( xi );
  free ( yd );
  free ( yi );

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

double r8_choose ( int n, int k )

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

    r8_choose() computes the combinatorial coefficient C(N,K).

  Discussion:

    Real arithmetic is used, and C(N,K) is computed directly, via
    Gamma functions, rather than recursively.

    C(N,K) is the number of distinct combinations of K objects
    chosen from a set of N distinct objects.  A combination is
    like a set, in that order does not matter.

    C(N,K) = N! / ( (N-K)! * K! )

  Example:

    The number of combinations of 2 things chosen from 5 is 10.

    C(5,2) = ( 5 * 4 * 3 * 2 * 1 ) / ( ( 3 * 2 * 1 ) * ( 2 * 1 ) ) = 10.

    The actual combinations may be represented as:

      (1,2), (1,3), (1,4), (1,5), (2,3),
      (2,4), (2,5), (3,4), (3,5), (4,5).

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    25 July 2011

  Author:

    John Burkardt

  Parameters:

    Input, int N, the value of N.

    Input, int K, the value of K.

    Output, double R8_CHOOSE, the value of C(N,K)
*/
{
  double arg;
  double fack;
  double facn;
  double facnmk;
  double value;

  if ( n < 0 )
  {
    value = 0.0;
  }
  else if ( k == 0 )
  {
    value = 1.0;
  }
  else if ( k == 1 )
  {
    value = ( double ) ( n );
  }
  else if ( 1 < k && k < n - 1 )
  {
    arg = ( double ) ( n + 1 );
    facn = lgamma ( arg );

    arg = ( double ) ( k + 1 );
    fack = lgamma ( arg );

    arg = ( double ) ( n - k + 1 );
    facnmk = lgamma ( arg );

    value = r8_nint ( exp ( facn - fack - facnmk ) );
  }
  else if ( k == n - 1 )
  {
    value = ( double ) ( n );
  }
  else if ( k == n )
  {
    value = 1.0;
  }
  else
  {
    value = 0.0;
  }

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

double r8_mop ( int i )

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

    r8_mop() returns the I-th power of -1 as an R8 value.

  Discussion:

    An R8 is an double value.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 July 2008

  Author:

    John Burkardt

  Parameters:

    Input, int I, the power of -1.

    Output, double R8_MOP, the I-th power of -1.
*/
{
  double value;

  if ( ( i % 2 ) == 0 )
  {
    value = + 1.0;
  }
  else
  {
    value = - 1.0;
  }

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

int r8_nint ( double x )

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

    r8_nint() returns the nearest integer to an R8.

  Example:

        X         R8_NINT

      1.3         1
      1.4         1
      1.5         1 or 2
      1.6         2
      0.0         0
     -0.7        -1
     -1.1        -1
     -1.6        -2

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    05 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, the value.

    Output, int R8_NINT, the nearest integer to X.
*/
{
  int s;
  int value;

  if ( x < 0.0 )
  {
    s = - 1;
  }
  else
  {
    s = + 1;
  }
  value = s * ( int ) ( fabs ( x ) + 0.5 );

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

double *r8vec_cheby1space_new ( int n, double a, double b )

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

    r8vec_cheby1space_new() creates Chebyshev Type 1 values in [A,B].

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    17 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double A, B, the interval.

    Output, double R8VEC_CHEBY1SPACE_NEW[N], a vector of Type 2
    Chebyshev spaced data.
*/
{
  double c;
  int i;
  const double r8_pi = 3.141592653589793;
  double theta;
  double *x;

  x = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    theta = ( double ) ( n - i ) * r8_pi / ( double ) ( n + 1 );

    c = cos ( theta );

    x[i] = ( ( 1.0 - c ) * a  
           + ( 1.0 + c ) * b ) 
           /   2.0;
  }

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

double *r8vec_cheby2space_new ( int n, double a, double b )

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

    r8vec_cheby2space_new() creates Chebyshev Type 2 values in [A,B].

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    05 July 2017

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double A, B, the interval.

    Output, double R8VEC_CHEBY1SPACE_NEW[N], a vector of Type 1
    Chebyshev spaced data.
*/
{
  double c;
  int i;
  const double r8_pi = 3.141592653589793;
  double theta;
  double *x;

  x = ( double * ) malloc ( n * sizeof ( double ) );

  if ( n == 1 )
  {
    x[0] = ( a + b ) / 2.0;
  }
  else
  {
    for ( i = 0; i < n; i++ )
    {
      theta = ( double ) ( n - i - 1 ) * r8_pi / ( double ) ( n - 1 );

      c = cos ( theta );

      if ( ( n % 2 ) == 1 )
      {
        if ( 2 * i + 1 == n )
        {
          c = 0.0;
        }
      }

      x[i] = ( ( 1.0 - c ) * a  
             + ( 1.0 + c ) * b ) 
             /   2.0;
    }
  }

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

double *r8vec_copy_new ( int n, double a1[] )

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

    r8vec_copy_new() copies an R8VEC.

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 August 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vectors.

    Input, double A1[N], the vector to be copied.

    Output, double R8VEC_COPY_NEW[N], the copy of A1.
*/
{
  double *a2;
  int i;

  a2 = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a1[i];
  }
  return a2;
}
/******************************************************************************/

double *r8vec_midspace_new ( int n, double a, double b )

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

    r8vec_midspace_new() creates a vector of linearly spaced values.

  Discussion:

    An R8VEC is a vector of R8's.

    This function divides the interval [a,b] into n subintervals, and then
    returns the midpoints of those subintervals.

  Example:

    N = 5, A = 10, B = 20
    X = [ 11, 13, 15, 17, 19 ]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    03 June 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double A, B, the endpoints of the interval.

    Output, double R8VEC_MIDSPACE_NEW[N], a vector of linearly spaced data.
*/
{
  double *x;
  int i;

  x = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    x[i] = ( ( double ) ( 2 * n - 2 * i - 1 ) * a 
           + ( double ) (         2 * i + 1 ) * b ) 
           / ( double ) ( 2 * n );
  }

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

double r8vec_norm_affine ( int n, double v0[], double v1[] )

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

    r8vec_norm_affine() returns the affine L2 norm of an R8VEC.

  Discussion:

    The affine vector L2 norm is defined as:

      R8VEC_NORM_AFFINE(V0,V1)
        = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2 )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 October 2010

  Author:

    John Burkardt

  Parameters:

    Input, int N, the dimension of the vectors.

    Input, double V0[N], the base vector.

    Input, double V1[N], the vector whose affine L2 norm is desired.

    Output, double R8VEC_NORM_AFFINE, the affine L2 norm of V1.
*/
{
  int i;
  double value;

  value = 0.0;

  for ( i = 0; i < n; i++ )
  {
    value = value + ( v1[i] - v0[i] ) * ( v1[i] - v0[i] );
  }
  value = sqrt ( value );

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

void r8vec2_print ( int n, double a1[], double a2[], char *title )

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

    r8vec2_print() prints an R8VEC2.

  Discussion:

    An R8VEC2 is a dataset consisting of N pairs of real values, stored
    as two separate vectors A1 and A2.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 March 2009

  Author:

    John Burkardt

  Input:

    int N, the number of components of the vector.

    double A1[N], double A2[N], the vectors 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, "  %4d: %14g  %14g\n", i, a1[i], a2[i] );
  }

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

void timestamp ( )

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

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    24 September 2003

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}