Listing 3 Program to demonstrate divided-difference algorithm

/* Divided Difference Interpolation */
/* Copyright 1994 by Philip Gage */

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

#define POINTS 10 /* Max number of data points */

/* Function prototypes */
void interpolate(int,int,int,double*,double*,double*);
double evaluate (int, int, double*, double);
void print_function (int, int, double*);
void print_table(int,int,int,double*,double*,double*);

/* Generate polynomial or exponential function
    polynomial   1=polynomial, 0=exponential
    degree       Degree of desired solution
    n            Number of data points
    x,y          Data point values
    coefficient  Output coefficient array */

void interpolate (int polynomial, int degree, int n,

double *coefficients)

{
  double xpower, data[POINTS], diff[POINTS];
  int i, j, k;

  /* Copy input y values to data array */
  for (i = 0; i < n; i++)
    data[i] = y[i];

  /* Initialize coefficients to identity element */
  for (i = 0; i <= degree; i++)
    if (polynomial)
      coefficients[i] = 0.0;
    else
      coefficients[i] = 1.0;

  */ Solve each high order term */
  for (i = degree; i >= 0; i--) {

    /* Copy data values to difference array */
    for (j = 0; j < n; j++)
      diff[j] = data[j];

    /* Compute divided differences or quotients */
    for (j= 1; j <= i; j++)
      for (k = 0; k < n-j; k++)
        if (polynomial)
          diff[k] = (diff[k+1] - diff[k]) /
                   (x[k+j] - x[k]);
        else
          diff[k] = pow(diff[k+1] / diff[k],
                      1.0/(x[k+j] - x[k]));

    /* Average current difference row */
    coefficients[i] = 0.0;
    for (j = 0; j < n-i; j++)
      coefficients[i] += diff[j];
    coefficients[i] /= n-i;

    /* Remove high order term from data */
    for (j = 0; j < n; j++) {
      /* Compute x[j] raised to i power */
      xpower = 1.0;
      for (k = 0; k< i; k++)
        xpower *= x[j];
      if (polynomial)
        data[j] -= coefficients[i] * xpower;
      else
        data[j] /= pow(coefficients[i],xpower);
    }
  }
}

/* Evaluate function at xvalue using Homer's rule */
double evaluate (int polynomial, int degree,
              double *coefficients, double xvalue)
{
  double yvalue;
  int i;

  if (polynomial)
    yvalue = 0.0;
  else
    yvalue = 1.0;
  for (i = degree; i >= 0; i--)
    if (polynomial)
      yvalue = yvalue * xvalue + coefficients[i];
    else
      yvalue = pow(yvalue,xvalue) * coefficients[i];
  return yvalue;
}

/* Print symbolic function terms using ^ for power */

void print_function (int polynomial, int degree,
                 double *coefficients)
{
  int i, identity, printing = 0;
  char operator1, operator2;

  /* Setup for polynomial or exponential */
  if (polynomial) {
    identity = 0.0;
    operator1 = '+';
    operator2 = '*';
  }
  else {
    identity = 1.0;
    operator1 = '*';
    operator2 = '^';
  }

  printf("\nF(X) = ");
  for (i = degree; i >= 0; i--) {
    if (coefficients[i] != identity) {
      if (printing)
        printf(" %c ",operator1);
      printf("%g",coefficients[i]);
      if (1 > 0)
        printf("%cX",operator2);
      if (i > 1)
        printf("^%d",i);
      printing = 1;
    }
  }
  printf("\n");
}

/* Print table of results */
void print_table (int polynomial, int degree, int n,
        double *x, double *y, double *coefficients)
{
  int i;
  double result, error, percent;

  printf("%15s%15s%15s%15s%15s %\n",
        "X","y","F(X)","Error","Error");
  for (i = 0; i < n; i++) {
    result = evaluate(polynomial,degree,
                   coefficients,x[i]);
    error = result - y[i];
    if (y[i] != 0.0)
      percent = 100.0 * error / y[i];
    printf("%l5g%l5g%15g%15g%15g %\n\n",
          x[i],y[i],result,error,percent);
  }
}

void main()
{
  int i, n, degree, polynomial, choice;
  double x[POINTS], y[POINTS], coefficients[POINTS];

  /* Read number of data points from 2 to POINTS */
  do {
    printf("Number of data points: ");
    scanf("%d",&n);
  } while (n<2 || n>POINTS);

  /* Read each data point from user */
  for (i = 0; i < n; i++) {
    printf("x%d: ",i);
    scanf("%lf",&x[i]>;
    printf("y%d: ",i);
    scanf("%lf",&y[i]);
  }

  /* Generate solutions until Quit selected */
  while (1) {
    printf("l Polynomial interpolation\n");
    printf("2 Exponential interpolation\n");
    printf("3 Quit program\n");
    printf("Choice: ");
    scanf("%1d",&choice);
    if (choice < 1 || choice > 2) break;
    polynomial = (choice == 1);

    /* Read degree of solution from 1 to n-1 */
    do {
      printf("Degree: ");
      scanf("%d",&degree);
    } while (degree < 0 || degree >= n);

    /* Compute function and print results */
    interpolate(polynomial,degree,n,x,y,coefficients);
    print_function(polynomial, degree, coefficients);
    print_table(polynomial,degree,n,x,y,coefficients);
  }
}

/* End of File */