/*************************************************************************/
/* Program file name: lagrange.c                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 5, 2006                                           */
/*                                                                       */
/* (1) This C program is created for the book, "An Introduction to       */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

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

/* An example of extracting an approximate function
   value via the Lagrange interpolation scheme. */

main() {
  int i;
  double x = 0.9, f, df;
  double xi[] = {0, 0.5, 1, 1.5, 2};
  double fi[] = {1.0, 0.938470, 0.765198, 0.511828, 0.223891};
  void aitken();

  aitken(nmax, xi, fi, x, &f, &df);
  printf("%10.6lf %10.6lf %10.6lf\n", x, f, df);
}

/* Function to carry out the Aitken recursions. */

void aitken(n, xi, fi, x, f, df) 
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;
}