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

using namespace std;

# include "divdif.H"

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

void data_to_dif ( int ntab, double xtab[], double ytab[], double diftab[] )

//****************************************************************************80
//
//  Purpose:
//
//    DATA_TO_DIF sets up a divided difference table from raw data.
//
//  Discussion:
//
//    Space can be saved by using a single array for both the DIFTAB and
//    YTAB dummy parameters.  In that case, the difference table will
//    overwrite the Y data without interfering with the computation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of pairs of points
//    (XTAB[I],YTAB[I]) which are to be used as data.  
//
//    Input, double XTAB[NTAB], the X values at which data was taken.
//    These values must be distinct.
//
//    Input, double YTAB[NTAB], the corresponding Y values.
//
//    Output, double DIFTAB[NTAB], the divided difference coefficients 
//    corresponding to the input (XTAB,YTAB).
//
{
  int i;
  int j;
//
//  Copy the data values into DIFTAB.
//
  for ( i = 0; i < ntab; i++ )
  {
    diftab[i] = ytab[i];
  }
//
//  Make sure the abscissas are distinct.
//
  for ( i = 0; i < ntab; i++ )
  {
    for ( j = i+1; j < ntab; j++ )
    {
      if ( xtab[i] - xtab[j] == 0.0 )
      {
        cout << "\n";
        cout << "DATA_TO_DIF - Fatal error!\n";
        cout << "  Two entries of XTAB are equal!\n";
        cout << "  XTAB[%d] = " << xtab[i] << "\n";
        cout << "  XTAB[%d] = " << xtab[j] << "\n";
        exit ( 1 );
      }
    }
  }
//
//  Compute the divided differences.
//
  for ( i = 1; i <= ntab-1; i++ )
  {
    for ( j = ntab-1; i <= j; j-- )
    {
      diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] );
    }
  }
 
  return;
}
//****************************************************************************80

void data_to_dif_display ( int ntab, double xtab[], double ytab[], 
  double diftab[] )

//****************************************************************************80
//
//  Purpose:
//
//    DATA_TO_DIF_DISPLAY sets up a divided difference table and prints out intermediate data.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of pairs of points
//    (XTAB[I],YTAB[I]) which are to be used as data.  
//
//    Input, double XTAB[NTAB], the X values at which data was taken.
//    These values must be distinct.
//
//    Input, double YTAB[NTAB], the corresponding Y values.
//
//    Output, double DIFTAB[NTAB], the divided difference coefficients 
//    corresponding to the input (XTAB,YTAB).
//
{
  int i;
  int j;

  if ( !r8vec_distinct ( ntab, xtab ) )
  {
    cout << "\n";
    cout << "DATA_TO_DIF_DISPLAY - Fatal error!\n";
    cout << "  Two entries of XTAB are equal!\n";
    return;
  }

  cout << "\n";
  cout << "  The divided difference table:\n";
  cout << "\n";
  cout << "        ";
  for ( i = 0; i < ntab; i++ )
  {
    cout << setw(10) << xtab[i] << "  ";
  }
  cout << "\n";
  cout << "\n";
  cout << setw(6) << 0 << "  ";
  for ( i = 0; i < ntab; i++ )
  {
    cout << setw(10) << ytab[i] << "  ";
  }
  cout << "\n";
//
//  Copy the data values into DIFTAB.
//
  for ( i = 0; i < ntab; i++ )
  {
    diftab[i] = ytab[i];
  }
//
//  Compute the divided differences.
//
  for ( i = 1; i <= ntab-1; i++ )
  {
    cout << setw(6) << i << "  ";
    for ( j = ntab-1; i <= j; j-- )
    {
      diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] );
    }
    for ( j = i; j < ntab; j++ )
    {
      cout << setw(10) << diftab[j] << "  ";
    }
    cout << "\n";
  }
 
  return;
}
//****************************************************************************80

void data_to_r8poly ( int ntab, double xtab[], double ytab[], double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    DATA_TO_R8POLY computes the coefficients of a polynomial interpolating data.
//
//  Discussion:
//
//    Space can be saved by using a single array for both the C and
//    YTAB parameters.  In that case, the coefficients will
//    overwrite the Y data without interfering with the computation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of data points.
//
//    Input, double XTAB[NTAB], YTAB[NTAB], the data values.
//
//    Output, double C[NTAB], the coefficients of the polynomial that passes
//    through the data (XTAB,YTAB).  C(0) is the constant term.
//
{
  if ( !r8vec_distinct ( ntab, xtab ) )
  {
    cout << "\n";
    cout << "DATA_TO_R8POLY - Fatal error!\n";
    cout << "  Two entries of XTAB are equal.\n";
    return;
  }

  data_to_dif ( ntab, xtab, ytab, c );

  dif_to_r8poly ( ntab, xtab, c, c );

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

void dif_antideriv ( int ntab, double xtab[], double diftab[], int *ntab2,
  double xtab2[], double diftab2[] )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_ANTIDERIV integrates a polynomial in divided difference form.
//
//  Discussion:
//
//    This routine uses the divided difference representation (XTAB, DIFTAB)
//    of a polynomial to compute the divided difference representation
//    (XTAB, ANTTAB) of the antiderivative of the polynomial.
//
//    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 GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the size of the input table.
//
//    Input, double XTAB[NTAB], the abscissas for the divided
//    difference table.
//
//    Input, double DIFTAB[NTAB], the divided difference table.
//
//    Output, int *NTAB2, the size of the output table, which is NTAB+1.
//
//    Input, double XTAB2[NTAB2], the abscissas for the divided
//    difference table for the antiderivative.
//
//    Output, double DIFTAB2[NTAB2], the divided difference 
//    table for the antiderivative.
//
{
  int i;
  double *xtab1;
  double *diftab1;
//
//  Using a temporary copy of the difference table, shift the
//  abscissas to zero.
//
  xtab1 = new double[ntab];
  diftab1 = new double[ntab];

  for ( i = 0; i < ntab; i++ )
  {
    xtab1[i] = xtab[i];
  }
  for ( i = 0; i < ntab; i++ )
  {
    diftab1[i] = diftab[i];
  }

  dif_shift_zero ( ntab, xtab1, diftab1 );

  cout << flush;
//
//  Append a final zero to XTAB.
//
  *ntab2 = ntab + 1;
  for ( i = 0; i < *ntab2; i++ )
  {
    xtab2[i] = 0.0;
  }
//
//  Get the antiderivative of the standard form polynomial.
//
  r8poly_ant_cof ( ntab, diftab1, diftab2 );

  delete [] xtab1;
  delete [] diftab1;

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

void dif_append ( int ntab, double xtab[], double diftab[], double xval, 
  double yval, int *ntab2, double xtab2[], double diftab2[] )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_APPEND adds a pair of data values to a divided difference table.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the size of the difference table.
//
//    Input, double XTAB[NTAB], the abscissas of the table.
//
//    Input, double DIFTAB[NTAB], the difference table. 
//
//    Input, double XVAL, the data abscissa to be added to the table.
//
//    Input, double YVAL, the data value to be added to the table.
//
//    Output, int *NTAB2, the updated size of the difference table.
//
//    Output, double XTAB2[*NTAB2], the updated abscissas.
//
//    Output, double DIFTAB2[*NTAB2], the updated difference table.
//
{
  int i;

  *ntab2 = ntab + 1;
//
//  Move the original data up one index.
//
  for ( i = *ntab2-1; 1 <= i; i-- )
  {
    diftab2[i] = diftab[i-1];
    xtab2[i] = xtab[i-1];
  }
//
//  Recompute the data.
//
  xtab2[0] = xval;
  diftab2[0] = yval;

  for ( i = 1; i < *ntab2; i++ )
  {
    diftab2[i] = ( diftab2[i] - diftab2[i-1] ) / ( xtab2[i] - xtab2[0] );
  }
 
  return;
}
//****************************************************************************80********

void dif_basis ( int ntab, double xtab[], double *diftab )

//****************************************************************************80********
//
//  Purpose:
//
//    DIF_BASIS computes all Lagrange basis polynomials in divided difference form.
//
//  Discussion:
//
//    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
//    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
//    XTAB(J) for J not equal to I, and 1 when J is equal to I.
//
//    The Lagrange basis polynomials have the property that the interpolating
//    polynomial through a set of NTAB data points (XTAB,YTAB) may be
//    represented as
//
//      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
//
//    Higher order interpolation at selected points may be accomplished
//    using repeated X values, and scaled derivative values.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of X data points XTAB, and the number of
//    basis polynomials to compute.
//
//    Input, double XTAB[NTAB], the X values upon which the Lagrange basis
//    polynomials are to be based.
//
//    Output, double DIFTAB[NTAB*NTAB], points to a list of NTAB * NTAB values, 
//    the set of divided difference tables, stored as consecutive rows.
//    Logical row I of DIFTAB contains the table for the I-th Lagrange basis
//    polynomial.
//
{
  int i;
  int j;
  double *pointer1;
  double *pointer2;

  pointer1 = diftab;

  for ( i = 0; i < ntab; i++ )
  {

    pointer2 = pointer1;

    for ( j = 0; j < ntab; j++ )
    {
      if ( j == i )
      {
        *pointer1 = 1.0;
      } 
      else
      {
        *pointer1 = 0.0;
      } 
      pointer1++;
    }

    data_to_dif ( ntab, xtab, pointer2, pointer2 );
  }

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

void dif_basis_i ( int ival, int ntab, double xtab[], double diftab[] )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_BASIS_I computes the I-th Lagrange basis polynomial in divided difference form.
//
//  Discussion:
//
//    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
//    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
//    XTAB(J) for J not equal to I, and 1 when J is equal to I.
//
//    The Lagrange basis polynomials have the property that the interpolating
//    polynomial through a set of NTAB data points (XTAB,YTAB) may be
//    represented as
//
//      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
//
//    Higher order interpolation at selected points may be accomplished
//    using repeated X values, and scaled derivative values.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int IVAL, the index of the desired Lagrange basis polynomial.
//    IVAL should be between 1 and NTAB.
//
//    Input, int NTAB, the number of data points XTAB.
//
//    Input, double XTAB[NTAB], the X values upon which the Lagrange basis
//    polynomial is to be based.
//
//    Output, double DIFTAB[NTAB], the divided difference table for the IVAL-th
//    Lagrange basis polynomial.
//
{
  int i;
//
//  Check IVAL.
//
  if ( ival < 1 || ntab < ival )
  {
    cout << "\n";
    cout << "DIF_BASIS_I - Fatal error!\n";
    cout << "  IVAL must be between 1 and " << ntab << ".\n";
    cout << "  but your value is " << ival << "\n";
    return;
  }
//
//  Initialize DIFTAB to Delta(I,J).
//
  for ( i = 0; i <= ntab-1; i++ )
  {
    diftab[i] = 0.0;
  }
  diftab[ival] = 1.0;
//
//  Compute the IVAL-th Lagrange basis polynomial.
//
  data_to_dif ( ntab, xtab, diftab, diftab );

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

void dif_deriv ( int ntab, double xtab[], double diftab[], int *ntab2,
  double xtab2[], double diftab2[] )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_DERIV computes the derivative of a polynomial in divided difference form.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the size of the input table.
//
//    Input, double XTAB[NTAB], the abscissas for the divided
//    difference table.
//
//    Input, double DIFTAB[NTAB], the divided difference table.
//
//    Output, int *NTAB2, the size of the output table, which is NTAB-1.
//
//    Input, double XTAB2[NTAB2], the abscissas for the divided
//    difference table for the derivative.
//
//    Output, double DIFTAB2[NTAB2], the divided difference 
//    table for the derivative.
//
{
  int i;
  double *xtab1;
  double *diftab1;
//
//  Using a temporary copy of the difference table, shift the
//  abscissas to zero.
//
  xtab1 = new double[ntab];
  diftab1 = new double[ntab];

  for ( i = 0; i < ntab; i++ )
  {
    xtab1[i] = xtab[i];
  }
  for ( i = 0; i < ntab; i++ )
  {
    diftab1[i] = diftab[i];
  }

  dif_shift_zero ( ntab, xtab1, diftab1 );
//
//  Construct the derivative.
//
  *ntab2 = ntab - 1;

  for ( i = 0; i < *ntab2; i++ )
  {
    xtab2[i] = 0.0;
  }

  for ( i = 0; i < *ntab2; i++ )
  {
    diftab2[i] = ( double ) ( i + 1 ) * diftab1[i+1];
  }

  delete [] xtab1;
  delete [] diftab1;

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

void dif_print ( int ntab, double xtab[], double diftab[], char *title  )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_PRINT prints the polynomial represented by a divided difference table.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the dimension of the arrays DIFTAB and XTAB.
//
//    Input, double XTAB[NTAB], the X values for the polynomial.
//
//    Input, double DIFTAB[NTAB], the divided difference table
//    for the polynomial.
//
//    Input, char *TITLE, a title to be printed first.
//    TITLE may be blank.
//
{
  int i;

  if ( s_len_trim ( title ) != 0 )
  {
    cout << "\n";
    cout << title << "\n";
  }
  cout << "\n";
  cout << "  p(x) =                       " 
       << setw(14) << diftab[0] << "\n";

  for ( i = 1; i < ntab; i++ )
  {
    cout << "       + ( x - " 
         << setw(10) << xtab[i-1] << ") * ( " 
         << setw(14) << diftab[i] << "\n";
  }

  cout << "        ";
  for ( i = 1; i < ntab; i++ )
  {
    cout << ")";
  }
  cout << "\n";

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

void dif_shift_zero ( int ntab, double xtab[], double diftab[] )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_SHIFT_ZERO shifts a divided difference table so that all abscissas are zero.
//
//  Discussion:
//
//    When the abscissas are changed, the coefficients naturally
//    must also be changed.
//
//    The resulting pair (XTAB, DIFTAB) still represents the
//    same polynomial, but the entries in DIFTAB are now the
//    standard polynomial coefficients.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the length of the DIFTAB and XTAB arrays.
//
//    Input/output, double XTAB[NTAB], the X values that correspond to the
//    divided difference table.  On output, XTAB contains only zeroes.
//
//    Input/output, double DIFTAB[NTAB], the divided difference table
//    for the polynomial.  On output, DIFTAB is also
//    the coefficient array for the standard representation
//    of the polynomial.
//
{
  int i;
  double xval;

  xval = 0.0;

  for ( i = 1; i <= ntab; i++ )
  {
    dif_shift_x ( ntab, xtab, diftab, xval );
  }

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

void dif_shift_x ( int ntab, double xtab[], double diftab[], double xval )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_SHIFT_X replaces one abscissa of a divided difference table with a new one.
//
//  Discussion:
//
//    This routine shifts the representation of a divided difference polynomial by
//    dropping the last X value in XTAB, and adding a new X value to the
//    beginning of the XTAB array, suitably modifying the coefficients stored
//    in DIFTAB.
//
//    The representation of the polynomial is changed, but the polynomial itself
//    should be identical.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of divided difference coefficients, and
//    the number of entries in XTAB.
//
//    Input/output, double XTAB[NTAB], the X values used in the representation of
//    the divided difference polynomial.
//    After a call to this routine, the last entry of XTAB has been dropped, the other
//    entries have shifted up one index, and XVAL has been inserted at the
//    beginning of the array.
//
//    Input/output, double DIFTAB[NTAB], the divided difference coefficients
//    corresponding to the XTAB array.  On output, this array has been
//    adjusted.
//
//    Input, double XVAL, a new X value which is to be used in the representation
//    of the polynomial.  On output, XTAB[0] equals XVAL and the representation
//    of the polynomial has been suitably changed.
//    Note that XVAL does not have to be distinct from any of the original XTAB
//    values.
//
{
  int i;
//
//  Recompute the divided difference coefficients.
//
  for ( i = ntab-2; 0 <= i; i-- )
  {
    diftab[i] = diftab[i] + ( xval - xtab[i] ) * diftab[i+1];
  }
//
//  Shift the X values up one position.
//
  for ( i = ntab-1; 0 < i; i-- )
  {
    xtab[i] = xtab[i-1];
  }

  xtab[0] = xval;

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

void dif_to_r8poly ( int ntab, double xtab[], double diftab[], double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_TO_R8POLY converts a divided difference table to a standard polynomial.
//
//  Discussion:
//
//    The vector DIFTAB, containing the divided difference polynomial
//    coefficients is overwritten with the standard form polynomial
//    coefficients, but the abscissa vector XTAB is unchanged.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of coefficients, and abscissas.
//
//    Input, double XTAB[NTAB], the X values used in the divided difference
//    representation of the polynomial.
//
//    Input, double DIFTAB[NTAB], the divided difference table.
//
//    Output, double C[NTAB], the standard form polyomial coefficients.
//    C[0] is the constant term, and C[NTAB-1] is the coefficient
//    of X**(NTAB-1).
//
{
  int i;
  int j;

  for ( i = 0; i < ntab; i++ )
  {
    c[i] = diftab[i];
  }
//
//  Recompute the divided difference coefficients.
//
  for ( j = 1; j <= ntab-1; j++ )
  {
    for ( i = 1; i <= ntab-j; i++ )
    {
      c[ntab-i-1] = c[ntab-i-1] - xtab[ntab-i-j] * c[ntab-i];
    }
  }

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

double dif_val ( int ntab, double xtab[], double diftab[], double xval )

//****************************************************************************80
//
//  Purpose:
//
//    DIF_VAL evaluates a divided difference polynomial at a point.
//
//  Discussion:
//
//    DATA_TO_DIF must be called first to set up the divided difference table.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, integer NTAB, the number of divided difference
//    coefficients in DIFTAB, and the number of points XTAB.
//
//    Input, double XTAB[NTAB], the X values upon which the
//    divided difference polynomial is based.
//
//    Input, double DIFTAB[NTAB], the divided difference table.
//
//    Input, double XVAL, a value of X at which the polynomial
//    is to be evaluated.
//
//    Output, double DIF_VAL, the value of the polynomial at XVAL.
//
{
  int i;
  double value;

  value = diftab[ntab-1];
  for ( i = 2; i <= ntab; i++ )
  {
    value = diftab[ntab-i] + ( xval - xtab[ntab-i] ) * value;
  }
 
  return value;
}
//****************************************************************************80

void nc_rule ( int norder, double a, double b, double xtab[], double weight[] )

//****************************************************************************80
//
//  Purpose:
//
//    NC_RULE computes the weights of a Newton-Cotes quadrature rule.
//
//  Discussion:
//
//    For the interval [A,B], the Newton-Cotes quadrature rule estimates
//
//      Integral ( A <= X <= B ) F(X) dX
//
//    using NORDER equally spaced abscissas XTAB(I) and a weight vector
//    WEIGHT(I):
//
//      Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ).
//
//    For the CLOSED rule, the abscissas include the points A and B.
//    For the OPEN rule, the abscissas do not include A and B.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NORDER, the order of the rule.
//
//    Input, double A, B, the left and right endpoints of the interval
//    over which the quadrature rule is to be applied.
//
//    Input, double XTAB[NORDER], the abscissas of the rule.
//
//    Output, double WEIGHT[NORDER], the weights of the rule.
//
{
  int i;
  double *poly_cof;
  double yvala;
  double yvalb;
//
//  Allocate temporary space for POLY_COF.
//
  poly_cof = new double[norder];

  for ( i = 1; i <= norder; i++ )
  {
//
//  Compute the Lagrange basis polynomial which is 1 at XTAB(I),
//  and zero at the other nodes.
//
    r8poly_basis_1 ( i, norder, xtab, poly_cof );
//
//  Evaluate the antiderivative of the polynomial at the left and
//  right endpoints.
//
    yvala = r8poly_ant_val ( norder-1, poly_cof, a );

    yvalb = r8poly_ant_val ( norder-1, poly_cof, b );

    weight[i-1] = yvalb - yvala;
  }
//
//  Free up POLY_COF space.
//
  delete [] poly_cof;

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

void ncc_rule ( int norder, double xtab[], double weight[] )

//****************************************************************************80
//
//  Purpose:
//
//    NCC_RULE computes the coefficients of a Newton-Cotes closed quadrature rule.
//
//  Discussion:
//
//    For the interval [-1,1], the Newton-Cotes quadrature rule estimates
//
//      Integral ( -1 <= X <= 1 ) F(X) dX
//
//    using NORDER equally spaced abscissas XTAB(I) and a weight vector
//    WEIGHT(I):
//
//      Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ).
//
//    For the CLOSED rule, the abscissas include A and B.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 September 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NORDER, the order of the rule.
//
//    Output, double XTAB[NORDER], the abscissas of the rule.
//
//    Output, double WEIGHT[NORDER], the weights of the rule.
//
{
  double a;
  double b;
  int i;
//
//  Compute a closed quadrature rule.
//
  a = -1.0;
  b = 1.0;

  for ( i = 1; i <= norder; i++ )
  {
    xtab[i-1] = ( ( double ) ( norder - i ) * a + ( double ) ( i - 1 ) * b )
      / ( double ) ( norder - 1 );
  }

  nc_rule ( norder, a, b, xtab, weight );

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

void nco_rule ( int norder, double xtab[], double weight[] )

//****************************************************************************80
//
//  Purpose:
//
//    NCO_RULE computes the coefficients of a Newton-Cotes open quadrature rule.
//
//  Discussion:
//
//    For the interval [-1,1], the Newton-Cotes quadrature rule estimates
//
//      Integral ( -1 <= X <= 1 ) F(X) dX
//
//    using NORDER equally spaced abscissas XTAB(I) and a weight vector
//    WEIGHT(I):
//
//      Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) ).
//
//    For the OPEN rule, the abscissas do not include A and B.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NORDER, the order of the rule.
//
//    Output, double XTAB[NORDER], the abscissas of the rule.
//
//    Output, double WEIGHT[NORDER], the weights of the  rule.
//
{
  double a;
  double b;
  int i;

  a = -1.0;
  b = 1.0;

  for ( i = 1; i <= norder; i++ )
  {
    xtab[i-1] = ( ( double ) ( norder + 1 - i ) * a + ( double ) ( i ) * b ) 
      / ( double ) ( norder + 1 );
  }

  nc_rule ( norder, a, b, xtab, weight );

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

void r8_swap ( double *x, double *y )

//****************************************************************************80
//
//  Purpose:
//
//    R8_SWAP swaps two R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, double *X, *Y.  On output, the values of X and
//    Y have been interchanged.
//
{
  double z;

  z = *x;
  *x = *y;
  *y = z;

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

void r8poly_ant_cof ( int n, double poly_cof[], double poly_cof2[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_ANT_COF integrates an R8POLY 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 GNU LGPL license. 
//
//  Modified:
//
//    13 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.
//
//    Input, double POLY_COF[N], the polynomial coefficients.
//    POLY_COF[0] is the constant term, and POLY_COF[N-1] is the
//    coefficient of X**(N-1).
//
//    Output, double POLY_COF2[N+1], the coefficients of the antiderivative
//    polynomial, in standard form.  The constant term is set to zero.
//
{
  int i;
//
//  Set the constant term.
//
  poly_cof2[0] = 0.0;
//
//  Integrate the polynomial.
//
  for ( i = 1; i <= n; i++ )
  {
    poly_cof2[i] = poly_cof[i-1] / ( double ) i;
  }

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

double r8poly_ant_val ( int n, double poly_cof[], double xval )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_ANT_VAL evaluates the antiderivative of an R8POLY in standard form.
//
//  Discussion:
//
//    The constant term of the antiderivative is taken to be zero.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.
//
//    Input, double POLY_COF[N], the polynomial coefficients.  POLY_COF[0]
//    is the constant term, and POLY_COF[N-1] is the coefficient of X**(N-1).
//
//    Input, double XVAL, the point where the antiderivative is to be
//    evaluated.
//
//    Output, double R8POLY_ANT_VAL, the value of the antiderivative of the polynomial
//    at XVAL.
//
{
  int i;
  double value;

  value = 0.0;

  for ( i = n-1; 0 <= i; i-- )
  {
    value = ( value + poly_cof[i] / ( double ) ( i + 1 ) ) * xval;
  }

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

void r8poly_basis ( int ntab, double xtab[], double *poly_cof )

//****************************************************************************80********
//
//  Purpose:
//
//    R8POLY_BASIS computes all Lagrange basis polynomials in standard form.
//
//  Discussion:
//
//    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
//    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
//    XTAB(J) for J not equal to I, and 1 when J is equal to I.
//
//    The Lagrange basis polynomials have the property that the interpolating
//    polynomial through a set of NTAB data points (XTAB,YTAB) may be
//    represented as
//
//      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
//
//    Higher order interpolation at selected points may be accomplished
//    using repeated X values, and scaled derivative values.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NTAB, the number of data points XTAB.
//
//    Input, double XTAB[NTAB], the X values upon which the Lagrange basis
//    polynomial is to be based.
//
//    Output, double *POLY_COF, points to NTAB * NTAB values.  The polynomial 
//    coefficients for the I-th Lagrange basis polynomial are stored in 
//    (logical) row I.  POLY_COF[0,*] is the constant term, and POLY_COF[NTAB-1,*] is 
//    the coefficient of X**(NTAB-1).
//
{
  int i;
  int j;
  double *pointer1;
  double *pointer2;

  pointer1 = poly_cof;

  for ( i = 0; i < ntab; i++ )
  {
    pointer2 = pointer1;

    for ( j = 0; j < ntab; j++ )
    {
      if ( j == i )
      {
        *pointer1 = 1.0;
      }
      else
      {
        *pointer1 = 0.0;
      } 
      pointer1++;
    }
//
//  Compute the divided difference table for the IVAL-th Lagrange basis
//  polynomial.
//
    data_to_dif ( ntab, xtab, pointer2, pointer2 );
//
//  Convert the divided difference table coefficients to standard polynomial
//  coefficients.
//
    dif_to_r8poly ( ntab, xtab, pointer2, pointer2 );

  }

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

void r8poly_basis_1 ( int ival, int ntab, double xtab[], double poly_cof[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_BASIS_1 computes the I-th Lagrange basis polynomial in standard form.
//
//  Discussion:
//
//    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
//    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
//    XTAB(J) for J not equal to I, and 1 when J is equal to I.
//
//    The Lagrange basis polynomials have the property that the interpolating
//    polynomial through a set of NTAB data points (XTAB,YTAB) may be
//    represented as
//
//      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
//
//    Higher order interpolation at selected points may be accomplished
//    using repeated X values, and scaled derivative values.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int IVAL, the index of the desired Lagrange basis polynomial.
//    IVAL should be between 1 and NTAB.
//
//    Input, int NTAB, the number of data points XTAB.
//
//    Input, double XTAB[NTAB], the X values upon which the Lagrange basis
//    polynomial is to be based.
//
//    Output, double POLY_COF[NTAB], the polynomial coefficients for the
//    IVAL-th Lagrange basis polynomial.
//
{
  int i;
//
//  Check IVAL.
//
  if ( ival < 1 || ntab < ival )
  {
    cout << "\n";
    cout << "R8POLY_BASIS_1 - Fatal error!\n";
    cout << "  IVAL must be between 1 and " << ntab << ".\n";
    cout << "  but your value is " << ival << ".\n";
    return;
  }
//
//  Initialize POLY_COF to the IVAL-th column of the identity matrix.
//
  for ( i = 0; i <= ntab-1; i++ )
  {
    poly_cof[i] = 0.0;
  }
  poly_cof[ival-1] = 1.0;
//
//  Compute the divided difference table for the IVAL-th Lagrange basis
//  polynomial.
//
  data_to_dif ( ntab, xtab, poly_cof, poly_cof );
//
//  Convert the divided difference table coefficients to standard polynomial
//  coefficients.
//
  dif_to_r8poly ( ntab, xtab, poly_cof, poly_cof );

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

void r8poly_der_cof ( int n, double poly_cof[], double poly_cof2[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_DER_COF computes the coefficients of the derivative of a real polynomial.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.
//
//    Input, double POLY_COF[N], the coefficients of the polynomial to
//    be differentiated.  POLY_COF[0] is the constant term, and
//    POLY_COF[N-1] is the coefficient of X**(N-1).
//
//    Output, double POLY_COF2[N-1], the coefficients of the derivative of
//    the polynomial.
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    poly_cof2[i] = ( double ) ( i + 1 ) * poly_cof[i+1];
  }

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

double r8poly_der_val ( int n, double poly_cof[], double xval )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_DER_VAL evaluates the derivative of a real polynomial in standard form.
//
//  Discussion:
//
//    A polynomial in standard form, with coefficients POLY_COF(*),
//    may be written:
//
//    P(X) = POLY_COF[0]
//         + POLY_COF[1] * X
//         ...
//         + POLY_COF[N-1] * X**(N-1)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.
//
//    Input, double POLY_COF[N], the polynomial coefficients.  POLY_COF[0]
//    is the constant term, and POLY_COF[N-1] is the coefficient of
//    X**(N-1).
//
//    Input, double XVAL, a value where the derivative of the polynomial
//    is to be evaluated.
//
//    Output, double R8POLY_DER_VAL, the value of the derivative of the polynomial
//    at XVAL.
//
{
  int i;
  double value;

  value = ( double ) ( n - 1 ) * poly_cof[n-1];

  for ( i = n-2; 1 <= i; i-- )
  {
    value = value * xval + ( double ) i * poly_cof[i];
  }

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

int r8poly_order ( int na, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_ORDER returns the order of a polynomial.
//
//  Discussion:
//
//    The order of a polynomial is the degree plus 1.
//
//    The order of a constant polynomial is 1.  The order of the
//    zero polynomial is debatable, but this routine returns the
//    order as 1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NA, the order of the polynomial (ignoring zero coefficients).
//
//    Input, double A[NA], the coefficients of the polynomial.
//
//    Output, int R8POLY_ORDER, the degree of the polynomial.
//
{
  int value;

  value = na;

  while ( 1 < value )
  {
    if ( a[value-1] != 0.0 )
    {
      return value;
    }

    value = value - 1;

  }

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

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

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_PRINT prints out a real polynomial in standard form.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.
//
//    Input, double POLY_COF[N], the polynomial coefficients.
//    POLY_COF[0] is the constant term and
//    POLY_COF[N-1] is the coefficient of X**(N-1).
//
//    Input, char *TITLE, a title to be printed first.
//    TITLE may be blank.
//
{
  int i;
  int n2;
//
//  Print the title, it if's not a blank.
//
  if ( s_len_trim ( title ) != 0 )
  {
    cout << "\n";
    cout << title << "\n";
  }
//
//  Determine the "real" degree of the polynomial.
//  N2 is the index of the highest order nonzero coefficient.
//
  n2 = r8poly_order ( n, poly_cof );

  cout << "\n";
  cout << "  p(x) = ";

  if ( poly_cof[n2-1] < 0.0 )
  {
    cout << "- ";
  }
  else
  {
    cout << "  ";
  }

  cout << setw(8) << fabs ( poly_cof[n2-1] );

  if ( 2 <= n2 )
  {
    cout << " * x";
  }

  if ( 3 <= n2 )
  {
    cout << " ^ " << n2-1;
  }
  cout << "\n";

  for ( i = n2-1; 1 <= i; i-- )
  {

    if ( poly_cof[i-1] < 0.0 )
    {
      cout << "       - ";
    }
    else if ( poly_cof[i-1] == 0.0 )
    {
      break;
    }
    else if ( 0.0 <= poly_cof[i-1] )
    {
      cout << "       + ";
    }

    cout << setw(10) << fabs ( poly_cof[i-1] );

    if ( 2 <= i )
    {
      cout <<" * x";
    }

    if ( 3 <= i )
    {
      cout << " ^ " << i-1;
    }

    cout << "\n";

  }

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

void r8poly_shift ( double scale, double shift, int n, double poly_cof[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_SHIFT adjusts the coefficients of a polynomial for a new argument.
//
//  Discussion:
//
//    Assuming P(X) is a polynomial in the argument X, of the form:
//
//      P(X) =
//          C(N-1) * X**(N-1)
//        + ...
//        + C(1) * X
//        + C(0),
//
//    and that Z is related to X by the formula:
//
//      Z = SCALE * X + SHIFT
//
//    then this routine computes coefficients C for the polynomial Q(Z):
//
//      Q(Z) =
//          C(N-1) * Z**(N-1)
//        + ...
//        + C(1) * Z
//        + C(0)
//
//    so that:
//
//      Q(Z(X)) = P(X)
//
//  Example:
//
//    P(X) = 2 * X**2 - X + 6
//
//    Z = 2.0 * X + 3.0
//
//    Q(Z) = 0.5 *         Z**2 -  3.5 * Z + 12
//
//    Q(Z(X)) = 0.5 * ( 4.0 * X**2 + 12.0 * X +  9 )
//            - 3.5 * (               2.0 * X +  3 )
//                                            + 12
//
//            = 2.0         * X**2 -  1.0 * X +  6
//
//            = P(X)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double SHIFT, SCALE, the shift and scale applied to X,
//    so that Z = SCALE * X + SHIFT.
//
//    Input, int N, the order of the polynomial.
//
//    Input/output, double POLY_COF[N].
//    On input, the coefficient array in terms of the X variable.
//    On output, the coefficient array in terms of the Z variable.
//
{
  int i;
  int j;

  for ( i = 1; i <= n; i++ )
  {
    for ( j = i+1; j <= n; j++ )
    {
      poly_cof[j-1] = poly_cof[j-1] / scale;
    }
  }

  for ( i = 1; i <= n; i++ )
  {
    for ( j = n-1; i <= j; j-- )
    {
      poly_cof[j-1] = poly_cof[j-1] - shift * poly_cof[j];
    }
  }

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

double r8poly_val_horner ( int n, double poly_cof[], double xval )

//****************************************************************************80
//
//  Purpose:
//
//    R8POLY_VAL_HORNER evaluates a real polynomial in standard form.
//
//  Discussion:
//
//    A polynomial in standard form, with coefficients POLY_COF(*),
//    may be written:
//
//    P(X) = POLY_COF[0]
//         + POLY_COF[1] * X
//         ...
//         + POLY_COF[N-1] * X**(N-1)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.
//
//    Input, double POLY_COF[N], the polynomial coefficients.  POLY_COF[0]
//    is the constant term, and POLY_COF[N-1] is the coefficient of
//    X**(N-1).
//
//    Input, double XVAL, a value where the polynomial is to be evaluated.
//
//    Output, double R8POLY_VAL_HORNER, the value of the polynomial at XVAL.
//
{
  int i;
  double value;

  value = poly_cof[n-1];

  for ( i = n-2; 0 <= i; i-- )
  {
    value = value * xval + poly_cof[i];
  }

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

bool r8vec_distinct ( int n, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_DISTINCT is true if the entries in an R8VEC are distinct.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input, double X[N], the vector to be checked.
//
//    Output, bool R8VEC_DISTINCT is true if all N elements of X are distinct.
//
{
  int i;
  int j;

  for ( i = 1; i <= n-1; i++ )
  {
    for ( j = 1; j <= i - 1; j++ )
    {
      if ( x[i] == x[j] )
      {
        return false;
      }
    }
  }

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

void r8vec_indicator ( int n, double a[] )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_INDICATOR sets an R8VEC to the indicator vector {1,2,3...}.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of elements of A.
//
//    Output, double A[N], the array to be initialized.
//
{
  int i;

  for ( i = 0; i <= n-1; i++ ) 
  {
    a[i] = ( double ) ( i + 1 );
  }

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

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

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_PRINT prints an R8VEC.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    16 December 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, double A[N], the vector to be printed.
//
//    Input, char *TITLE, a title to be printed first.
//    TITLE may be blank.
//
{
  int i;

  if ( s_len_trim ( title ) != 0 )
  {
    cout << "\n";
    cout << title << "\n";
  }

  cout << "\n";
  for ( i = 0; i <= n-1; i++ )
  {
    cout << setw(6)  << i+1  << "  " 
         << setw(10) << a[i] << "\n";
  }

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

void roots_to_dif ( int nroots, double roots[], int *ntab, double xtab[], 
  double diftab[] )

//****************************************************************************80
//
//  Purpose:
//
//    ROOTS_TO_DIF sets a divided difference table for a polynomial from its roots.
//
//  Discussion:
//
//    This turns out to be a simple task, because of two facts:
//
//    * The divided difference polynomial of one smaller degree which
//      passes through the values ( ROOT(I), 0 ) is the zero polynomial,
//      and hence has a zero divided difference table.
//
//    * We want a polynomial of one degree higher, but we don't want it
//      to pass through an addditional point.  Instead, we specify that
//      the polynomial is MONIC.  This means that the divided difference
//      table is almost the same as for the zero polynomial, except that
//      there is one more pair of entries, an arbitrary X value, and
//      a Y value of 1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 August 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NROOTS, is the number of roots.
//
//    Input, double ROOTS[NROOTS], the roots of the polynomial.
//
//    Output, int *NTAB, the size of the divided difference table.
//    This is NROOTS+1
//
//    Output, double XTAB[NTAB], the abscissas of the table.
//
//    Output, double DIFTAB[NTAB], the divided difference table.
//
{
  int i;

  *ntab = nroots + 1;
//
//  Build the appropriate difference table for the polynomial
//  through ( ROOTS(I), 0 ) of degree NTAB-2.
//
  for ( i = 0; i < *ntab-1; i++ )
  {
    diftab[i] = 0.0;
  }
  for ( i = 0; i < *ntab-1; i++ )
  {
    xtab[i] = roots[i];
  }
//
//  Append the extra data to make a monic polynomial of degree NTAB-1
//  which is zero at the NTAB-1 roots.
//
  xtab[*ntab-1] = 0.0;
  diftab[*ntab-1] = 1.0;

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

void roots_to_r8poly ( int nroots, double roots[], int *nc, double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    ROOTS_TO_R8POLY converts polynomial roots to polynomial coefficients.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int NROOTS, the number of roots specified.
//
//    Input, double ROOTS[NROOTS], the roots.
//
//    Output, integer *NC, the order of the polynomial, which will be NROOTS+1.
//
//    Output, double C[*NC], the coefficients of the polynomial.
//
{
  int i;
  int j;
  double *xtab;

  *nc = nroots + 1;
//
//  Initialize C to (0, 0, ..., 0, 1).
//  Essentially, we are setting up a divided difference table.
//
  xtab = new double[nroots+1];
  for ( i = 0; i < *nc-1; i++ )
  {
    xtab[i] = roots[i];
  }
  xtab[*nc-1] = 0.0;

  for ( i = 0; i < *nc-1; i++ )
  {
    c[i] = 0.0;
  }
  c[*nc-1] = 1.0;
//
//  Convert to standard polynomial form by shifting the abscissas
//  of the divided difference table to 0.
//
  dif_shift_zero ( *nc, xtab, c );

  delete [] xtab;

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

int s_len_trim ( char* s )

//****************************************************************************80********
//
//  Purpose:
//
//    S_LEN_TRIM returns the length of a string to the last nonblank.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    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 != ' ' )
    {
      return n;
    }
    t--;
    n--;
  }

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

void timestamp ( void )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    02 October 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

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

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

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

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}