/* divdifc.c  29 April 1999 */

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

#include "divdifc.h"

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

void adddif ( float diftab[], int ntab, float xtab[], float xval, 
  float ytab[], float yval ) {

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

/*
  Purpose:

    ADDDIF adds a pair of data values to a divided difference table.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float DIFTAB[NTAB].

    On input, entries 0 through NTAB-2 of DIFTAB contain divided 
    differences corresponding to the data in entries 0 through NTAB-1
    of XTAB and YTAB.

    On output, entries 0 through NTAB-1 of DIFTAB contain divided
    differences corresponding to the data in entries 0 through NTAB-1
    of XTAB and YTAB. 

    Input, int NTAB, the number of data values to be used, including
    the new value to be added.  NTAB may be 1, in which case it is
    assumed that the input DIFTAB, XTAB and YTAB are empty.

    Input/output, float XTAB[NTAB].
    On input, entries 0 through NTAB-2 contain data abscissas.
    On output, entries 0 through NTAB-1 contain data abscissas,
    and in particular, XTAB[0] = XVAL.

    Input, float XVAL, the data abscissa to be added to the table.

    Input/output, float YTAB[NTAB].
    On input, entries 0 through NTAB-2 contain data values.
    On output, entries 0 through NTAB-1 contain data values, and
    in particular, YTAB[0] = YVAL.

    Input, float YVAL, the data value to be added to the table.
*/
  int i;
/*
  Move the original data up one index.
*/
  for ( i = ntab-1; i >= 1; i-- ) {
    diftab[i] = diftab[i-1];
    xtab[i] = xtab[i-1];
    ytab[i] = ytab[i-1];
  }
/*
  Recompute the data.
*/
  xtab[0] = xval;
  ytab[0] = yval;
 
  diftab[0] = yval;
  for ( i = 1; i < ntab; i++ ) {
    diftab[i] = ( diftab[i] - diftab[i-1] ) / ( xtab[i] - xtab[0] );
  }
 
  return;
}
/**********************************************************************/

void antdif ( float diftab[], float anttab[], int ntab, float xtab[] ) {

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

/*
  Purpose:

    ANTDIF integrates a polynomial in divided difference form.

  Discussion:

    ANTDIF uses the divided difference representation (XTAB, DIFTAB)
    of a polynomial to compute the divided difference representation
    (XTAB, ANTTAB) of the antiderivative of the polynomial.

  Definition:

    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.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float DIFTAB[NTAB].

    On input, (XTAB,DIFTAB) is a divided difference table, which may
    be thought of as representing a polynomial through data points
    (XTAB,YTAB).

    On output, (XTAB,DIFTAB) represent the same polynomial, but the
    numeric values in DIFTAB will have changed, because XTAB is changed.

    Output, float ANTTAB[NTAB+1], the divided difference table computed for
    the antiderivative of the polynomial represented by (XTAB,DIFTAB).
    ANTTAB[0] is the arbitrary integration constant, and is set to zero,
    but may be reset to any desired value.  The entries of ANTTAB have
    been computed with X data based at zero.

    Input, int NTAB, the dimension of the array DIFTAB.

    Input/output, float XTAB[NTAB+1].  (Note the extra entry!)

    On input, entries 0 through NTAB-1 of XTAB contain the X locations where
    data was taken.

    On output, all the entries of XTAB have been set to zero.  This
    corresponds to the changes in DIFTAB, and to the way ANTTAB was
    computed.
*/

/*
  Recompute XTAB and DIFTAB, so that XTAB(1:NTAB) is entirely zero.
*/
  zerdif ( diftab, ntab, xtab );
/*
  Append a final zero to XTAB.
*/
  xtab[ntab] = 0.0;
/*
  Get the antiderivative of the standard form polynomial.
*/
  poly_ant_cof ( ntab-1, diftab, anttab );

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

void dif_to_poly ( float diftab[], int ntab, float xtab[] ) {

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

/*
  Purpose:

    DIF_TO_POLY converts a divided difference polynomial to standard form.

  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.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float DIFTAB[NTAB].

    On input, DIFTAB contains the coefficients of the divided difference
    polynomials, corresponding to the XTAB array.

    On output, DIFTAB contains the standard form polyomial coefficients.
    DIFTAB[0] is the constant term, and DIFTAB[NTAB-1] is the coefficient
    of X**(NTAB-1).

    Input, int NTAB, the number of coefficients, and abscissas.

    Input, float XTAB[NTAB], the X values used in the divided difference
    representation of the polynomial.
*/
  int i;
  int j;
/*
  Recompute the divided difference coefficients.
*/
  for ( j = 1; j <= ntab - 1; j++ ) {
    for ( i = ntab-1; i >= j; i-- ) {
      diftab[i-1] = diftab[i-1] - xtab[i-j] * diftab[i];
    }
  }

  return;
}

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

void derdif ( float diftab[], float dertab[], int ntab, float xtab[] ) {

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

/*
  Purpose:

    DERDIF computes the derivative of a polynomial in divided difference form.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float DIFTAB[NTAB], the divided difference
    table for the polynomial whose derivative is desired.

    On output, DIFTAB has been updated so that it is zero-based, but
    still represents the same polynomial.

    Output, float DERTAB[NTAB-1], contains the zero based divided
    difference table for the derivative of the polynomial represented by DIFTAB.

    Input, int NTAB, the dimension of the arrays DIFTAB, DERTAB, and XTAB.

    Input/output, float XTAB[NTAB].

    On input, XTAB contains the abscissas for the divided difference table
    stored in DIFTAB.

    On output, XTAB has been set to all zeroes.  
*/
  int i;
/*
  Move DIFTAB to zero based form.
*/
  zerdif ( diftab, ntab, xtab );

  for ( i = 0; i < ntab-1; i++ ) {
    dertab[i] = ( float ) ( i + 1 ) * diftab[i+1];
  }

  return;
}

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

void movdif ( float diftab[], int ntab, float xtab[], float xval ) {

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

/*
  Purpose:

    MOVDIF replaces one abscissa of a divided difference table with a new one.

  Discussion:

    MOVDIF 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
    is should be identical.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float DIFTAB[NTAB], the divided difference coefficients
    corresponding to the XTAB array.  On output, this array has been
    adjusted.

    Input, int NTAB, the number of divided difference coefficients, and
    the number of entries in XTAB.

    Input/output, float XTAB[NTAB], the X values used in the representation of
    the divided difference polynomial.

    After a call to MOVDIF, 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, float 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; i >= 0; i-- ) {
    diftab[i] = diftab[i] + ( xval - xtab[i] ) * diftab[i+1];
  }
/*
  Shift the X values up one position.
*/
  for ( i = ntab-1; i > 0; i-- ) {
    xtab[i] = xtab[i-1];
  }

  xtab[0] = xval;

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

void outdif ( float diftab[], int ntab, float xtab[], float ytab[] ) {

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

/*

  Purpose:

    OUTDIF sets up a divided difference table and prints out intermediate data.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Output, real DIFTAB[NTAB], the divided difference coefficients 
    corresponding to the input (XTAB,YTAB).

    Input, integer NTAB, the number of pairs of points
    (XTAB[I],YTAB[I]) which are to be used as data.  

    Input, real XTAB[NTAB], the X values at which data was taken.
    These values must be distinct.

    Input, real YTAB[NTAB], the corresponding Y values.
*/
  int i;
  int j;

  printf ( "\n" );
  printf ( "Divided difference table:\n" );
  printf ( "\n" );
  for ( i = 0; i < ntab; i++ ) {
    printf ( "%f ", xtab[i] );
  }
  printf ( "\n" );
  printf ( "\n" );
  printf ( "0 " );
  for ( i = 0; i < ntab; i++ ) {
    printf ( "%f ", ytab[i] );
  }
  printf ( "\n" );
/*
  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 ) {
        printf ( "\n" );
        printf ( "OUTDIF - Fatal error!\n" );
        printf ( "  Two entries of XTAB are equal!\n" );
        printf ( "  XTAB[%d] = %f\n", i, xtab[i] );
        printf ( "  XTAB[%d] = %f\n", j, xtab[j] );
        exit ( 1 );
      }
    }
  }
/*
  Compute the divided differences.
*/
  for ( i = 1; i <= ntab-1; i++ ) {
    printf ( "%d ", i );
    for ( j = ntab-1; j >= i; j-- ) {
      diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] );
    }
    for ( j = i; j < ntab; j++ ) {
      printf ( "%f ", diftab[j] );
    }
    printf ( "\n" );
  }
 
  return;
}
/**********************************************************************/

void poly_ant_cof ( int n, float poly_cof[], float poly_cof2[] ) {

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

/*
  Purpose:

    POLY_ANT_COF integrates a polynomial in standard form.

  Definition:

    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.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the degree of the original polynomial.

    Input, float 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, float POLY_COF2(0:N+1), the coefficients of the antiderivative
    polynomial, in standard form.  The constant term is set it to zero.
*/
  int i;
/*
  Set the constant term.
*/
  poly_cof2[0] = 0.0;
/*
  Integrate the polynomial.
*/
  for ( i = 1; i <= n+1; i++ ) {
    poly_cof2[i] = poly_cof[i-1] / ( float ) i;
  }

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

void poly_ant_val ( int n, float poly_cof[], float xval, float *yval ) {

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

/*
  Purpose:

    POLY_ANT_VAL evaluates the antiderivative of a polynomial in standard form.

  Discussion:

    The constant term of the antiderivative is taken to be zero.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the degree of the polynomial.

    Input, float POLY_COF[N+1], the polynomial coefficients.  POLY_COF[0]
    is the constant term, and POLY_COF[N] is the coefficient of
    X**N.

    Input, float XVAL, the point where the antiderivative is to be
    evaluated.

    Output, float *YVAL, the value of the antiderivative of the polynomial
    at XVAL.
*/

  int i;

  *yval = poly_cof[n] * xval / ( float ) ( n + 1 );

  for ( i = n-1; i >= 0; i-- ) {
    *yval = ( *yval + poly_cof[i] / ( float ) ( i + 1 ) ) * xval;
  }

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

void poly_der_cof ( int n, float poly_cof[], float poly_cof2[] ) {

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

/*
  Purpose:

    POLY_DER_COF computes the coefficients of the derivative of a polynomial.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the degree of the original polynomial.

    Input, float POLY_COF[N+1], the coefficients of the polynomial to
    be differentiated.  POLY_COF[0] is the constant term, and
    POLY_COF[N] is the coefficient of X**N.

    Output, float POLY_COF2[N], the coefficients of the derivative of
    the polynomial.
*/
  int i;

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

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

void poly_der_val ( int n, float poly_cof[], float xval, float *yval ) {

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

/*
  Purpose:

    POLY_VAL evaluates the derivative of a polynomial in standard form.

  Definition:

    A polynomial in standard form, with coefficients POLY_COF(*),
    may be written:

  P(X) = POLY_COF[0]
       + POLY_COF[1] * X
       ...
       + POLY_COF[N] * X**N

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the degree of the polynomial.

    Input, float POLY_COF[N+1], the polynomial coefficients.  POLY_COF[0]
    is the constant term, and POLY_COF[N] is the coefficient of
    X**N.

    Input, float XVAL, a value where the derivative of the polynomial
    is to be evaluated.

    Output, float *YVAL, the value of the derivative of the polynomial
    at XVAL.
*/
  int i;

  *yval = ( float) n * poly_cof[n];
  for ( i = n-1; i >= 1; i-- ) {
    *yval = (*yval) * xval + ( float ) i * poly_cof[i];
  }

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

void poly_print ( int n, float poly_cof[] ) {

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

/*
  Purpose:

    POLY_PRINT prints out a polynomial.

  Modified:

    29 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the degree of the polynomial.

    Input, float POLY_COF[N+1], the polynomial coefficients.
    POLY_COF[0] is the constant term and
    POLY_COF[N] is the coefficient of X**N.
*/
  int i;
  float mag;
  int n2;

  n2 = 0;
  for ( i = 0; i <= n; i++ ) {
    if ( poly_cof[i] != 0.0 ) {
      n2 = i;
    }
  }
 
  printf ( "\n" );
  printf ( "  p(x) = " );

  if ( poly_cof[n2] < 0.0 ) {
    printf ( "- " );
  }

  mag = fabs ( poly_cof[n2] );
  printf ( "%f", mag );

  if ( n2 >= 1 ) {
    printf ( " * x" );
  }

  if ( n2 >=2 ) {
    printf ( " ^ %d", n2 );
  }
  printf ( "\n" );

  for ( i = n2-1; i >= 0; i-- ) {

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

    mag = fabs ( poly_cof[i] );
    printf ( "%f", mag );

    if ( i >= 1 ) {
      printf (" * x" );
    }

    if ( i >= 2 ) {
      printf ( " ^ %d", i );
    }

    printf ( "\n" );

  }

  return;
}

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

void poly_val ( int n, float poly_cof[], float xval, float *yval ) {

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

/*
  Purpose:

    POLY_VAL evaluates a polynomial in standard form.

  Definition:

    A polynomial in standard form, with coefficients POLY_COF(*),
    may be written:

  P(X) = POLY_COF[0]
       + POLY_COF[1] * X
       ...
       + POLY_COF[N] * X**N

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, int N, the degree of the polynomial.

    Input, float POLY_COF[N+1], the polynomial coefficients.  POLY_COF[0]
    is the constant term, and POLY_COF[N] is the coefficient of
    X**N.

    Input, float XVAL, a value where the polynomial is to be evaluated.

    Output, float YVAL, the value of the polynomial at XVAL.
*/

  int i;

  *yval = poly_cof[n];

  for ( i = n-1; i >= 0; i-- ) {
    *yval = *yval * xval + poly_cof[i];
  }

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

void pridif ( float diftab[], int ntab, float xtab[] ) {

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

/*
  Purpose:

    PRIDIF prints the polynomial represented by a divided difference table.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, float DIFTAB[NTAB], the divided difference table
    for the polynomial.

    Input, int NTAB, the dimension of the arrays DIFTAB and XTAB.

    Input, float XTAB[NTAB], the X values for the polynomial.
*/
  int i;

  printf ( "\n" );
  printf ( "\n" );
  printf ( "Divided difference polynomial:\n" );
  printf ( "\n" );
  printf ( "  p(x) =                    %f\n", diftab[0] );

  for ( i = 1; i < ntab; i++ ) {
    printf ( " + ( x - %f) * ( %f\n", xtab[i-1], diftab[i] );
  }

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

  return;
}

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

void setdif ( float diftab[], int ntab, float xtab[], float ytab[] ) {

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

/*

  Purpose:

    SETDIF 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.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Output, real DIFTAB[NTAB], the divided difference coefficients 
    corresponding to the input (XTAB,YTAB).

    Input, integer NTAB, the number of pairs of points
    (XTAB[I],YTAB[I]) which are to be used as data.  

    Input, real XTAB[NTAB], the X values at which data was taken.
    These values must be distinct.

    Input, real YTAB[NTAB], the corresponding Y values.
*/
  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 ) {
        printf ( "\n" );
        printf ( "SETDIF - Fatal error!\n" );
        printf ( "  Two entries of XTAB are equal!\n" );
        printf ( "  XTAB[%d] = %f\n", i, xtab[i] );
        printf ( "  XTAB[%d] = %f\n", j, xtab[j] );
        exit ( 1 );
      }
    }
  }
/*
  Compute the divided differences.
*/
  for ( i = 1; i <= ntab-1; i++ ) {
    for ( j = ntab-1; j >= i; j-- ) {
      diftab[j] = ( diftab[j] - diftab[j-1] ) / ( xtab[j] - xtab[j-i] );
    }
  }
 
  return;
}
/**********************************************************************/

void valdif ( float diftab[], int ntab, float xtab[], float xval, 
  float *yval ) {

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

/*
  Purpose:

    VALDIF evaluates a divided difference polynomial at a point.

  Discussion:

    SETDIF must be called first to set up the divided difference table.

  Modified:

    12 April 1999

  Author:

    John Burkardt

  Parameters:

    Input, float DIFTAB[NTAB], the divided difference table.

    Input, integer NTAB, the number of divided difference
    coefficients in DIFTAB, and the number of points XTAB.

    Input, float XTAB[NTAB], the X values upon which the
    divided difference polynomial is based.

    Input, float XVAL, a value of X at which the polynomial
    is to be evaluated.

    Output, float *YVAL, the value of the polynomial at XVAL.
*/
  int i;

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

void zerdif ( float diftab[], int ntab, float xtab[] ) {

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

/*
  Purpose:

    ZERDIF 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.

  Modified:

    13 April 1999

  Author:

    John Burkardt

  Parameters:

    Input/output, float DIFTAB[NTAB], the divided difference table
    for the polynomial.  On output, DIFTAB is also
    the coefficient array for the standard representation
    of the polynomial.

    Input, int NTAB, the length of the DIFTAB and XTAB arrays.

    Input/output, float XTAB[NTAB], the X values that correspond to the
    divided difference table.  On output, XTAB contains only zeroes.
*/
  int i;
  float xval;

  xval = 0.0;

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

  return;
}