/*************************    Program 2.1    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 5

main()
/* Main program for the Lagrange interpolation with the Aitken method.
   NMAX is the number of given points.  Copyright (c) Tao Pang 1997. */
{
int i;
double h = 0.5, x = 0.9;
double f, df;
double xi[NMAX], fi[] = {1.0, 0.938470, 0.765198, 0.511828, 0.223891};
void aitken();

for (i = 0; i < NMAX; i++)
  {
  xi[i] = h*i;
  }
aitken(NMAX, xi, fi, x, &f, &df);
printf("%10.6lf %10.6lf %10.6lf\n", x, f, df);
}

void aitken(n, xi, fi, x, f, df)
/* Subroutine for the interpolation with the Aitken method.
   Copyright (c) Tao Pang 1997. */
int n;
double x;
double xi[], fi[];
double *f, *df;
{
int i, j;
double x1, x2, f1, f2;
double ft[NMAX];

for (i = 0; i <= (n-1); ++i)
  {
  ft[i] = fi[i];
  }
for (i = 0; i <= (n-2); ++i)
  {
  for (j = 0; j <= (n-i-2); ++j)
    {
    x1 = xi[j];
    x2 = xi[j+i+1];
    f1 = ft[j];
    f2 = ft[j+1];
    ft[j] = (x-x1)/(x2-x1)*f2+(x-x2)/(x1-x2)*f1;
    }
  }
*f  = ft[0];
*df = (fabs(*f-f1)+fabs(*f-f2))/2;
}