/*************************    Program 2.A    ****************************/
/*                                                                      */
/************************************************************************/
/* 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 21

main()
/* Main program for the Lagrange interpolation with the
   upward and downward correction method.
   Copyright (c) Tao Pang 1997. */
{
int n;
double xi[] = {0,0.5,1,1.5,2};
double fi[] = {1,0.938470,0.765198,0.511828,0.223891};
double x,f,df;
void updown();

n = 5;
x = 0.9;
updown(n,xi,fi,x,&f,&df);

printf("%16.8lf %16.8lf %16.8lf\n", x,f,df);
}

void updown (n,xi,fi,x,f,df)
/* Subroutine performing the Lagrange interpolation with the
   upward and downward correction method.  F: interpolated
   value.  DF: error estimated.  Copyright (c) Tao Pang 1997. */
int n;
double x,*f,*df;
double xi[],fi[];
{
  int i,i0,j,j0,k,it;
  double dx0,dx1,dt;
  double dp[NMAX][NMAX],dm[NMAX][NMAX];

  if (n > NMAX)
  {
    printf("Dimension too large\n");
    exit(1);
  }
  i0  = 0;
  dx0 = fabs(xi[n-1]-xi[0]);
  for (i = 0; i < n; ++i)
  {
    dp[i][i] = fi[i];
    dm[i][i] = fi[i];
    dx1 = fabs(x-xi[i]);
    if (dx1 < dx0)
    {
      i0  = i;
      dx0 = dx1;
    }
  }
  j0 = i0;

/* Evaluate correction matrices */

  for (i = 0; i < n-1; ++i)
  {
    for (j = 0; j < n-i-1; ++j)
    {
      k = j+i+1;
      dt = (dp[j][k-1]-dm[j+1][k])/(xi[k]-xi[j]);
      dp[j][k] = dt*(xi[k]-x);
      dm[j][k] = dt*(xi[j]-x);
    }
  }

/* Update the approximation */

  *f = fi[i0];
  it = 0;
  if (x < xi[i0]) it = 1;
  for (i = 0; i < n-1; ++i)
  {
    if ((it == 1) || (j0 == n-1))
    {
      i0  = i0-1;
      *df = dp[i0][j0];
      *f  = *f+*df;
      it = 0;
      if (j0 == n-1) it = 1;
    }
    else if ((it == 0) || (i0 == 0))
    {
      j0 = j0+1;
      *df = dm[i0][j0];
      *f  = *f+*df;
      it = 1;
      if (i0 == 0) it = 0;
    }
  }
  *df = fabs(*df);
}