//
// interpolation.c
//
// f(x) = 1 / (1 + x^2) の区間[-1,1]上での等間隔標本点Lagrange補間（反復補間法）
// 
// このファイルをテキストファイルにコピペして，拡張子を.cとして保存してください．
// 反復補間法の途中の差分商も表示されます．
// f(j, i) = f(x(1), ..., f(x(j-1)), f(x(i)); x) です．
//
#include <stdio.h>
#include <math.h>
#define NMAX 100
//
double interpolation(double x, double xj[NMAX+1], double fj[NMAX+1], int n_intpol)
{
  double yj[NMAX+1], gj[NMAX+1], ymin, y_tmp, g_tmp;
  double f_diff[NMAX+1][NMAX+1], fn, eps;
  int i, j, i1, l_conv;
  //
  eps = 1.0e-14;
  //
  for (j=1; j<=n_intpol; ++j)
    {
      yj[j] = x - xj[j];
      gj[j] = fj[j];
    }
  //
  l_conv = 0;
  for (i=1; i<=n_intpol; ++i)
    {
      ymin = fabs(yj[i]);
      for (j=i+1; j<=n_intpol; ++j)
	if (ymin > fabs(yj[j])) 
	  {
	    ymin = fabs(yj[j]);
	    i1 = j;
	  }
      if (i1!=i) 
	{
	  y_tmp = yj[i];
	  yj[i] = yj[i1];
	  yj[i1] = y_tmp;
	  g_tmp = gj[i];
	  gj[i] = gj[i1];
	  gj[i1] = g_tmp;
	}
      //      printf("|x(%3d) - x| = %22.15e\n", i, fabs(yj[i]));
      //
      f_diff[1][i] = gj[i];
      printf("f(%3d, %3d) = %22.15e\n", 1, i, f_diff[1][i]);
      for (j=2; j<=i; ++j)
	{
	  f_diff[j][i] = (yj[i] * f_diff[j-1][j-1] - yj[j-1] * f_diff[j-1][i]) / (yj[i] - yj[j-1]);
	  printf("f(%3d, %3d) = %22.15e\n", j, i, f_diff[j][i]);
	}
      printf("\n");
      if (i>1)
	if (fabs(f_diff[i][i] - f_diff[i-1][i-1]) < eps) 
	  {
	    l_conv = 1;
	    fn = f_diff[i][i];
	    break;
	  }
    }
  //
  if (l_conv == 0)
    {
      printf("Lagrange interpolation not convergent.\n");
      fn = f_diff[n_intpol][n_intpol];
    }
  //
  return fn;
}
//
void data_intpol(double (*func)(double), double xmin, double xmax, int n_intpol, 
		 double xj[NMAX+1], double fj[NMAX+1])
{
  double dx;
  int j;
  //
  dx = (xmax - xmin) / (n_intpol - 1);
  for (j=1; j<=n_intpol; ++j)
    {
      xj[j] = xmin + (j - 1) * dx;
      fj[j] = func(xj[j]);
    }
  //
  /* for (j=1; j<=n_intpol; ++j) */
  /*   printf("x(%3d) = %22.15e f(%3d) = %22.15e\n", j, xj[j], j, fj[j]); */
}
//
double func(double x)
{
  return 1.0 / ( 1.0 + x*x );
}
//
main()
{
  int n_intpol, k;
  double xj[NMAX+1], fj[NMAX+1], xmin, xmax;
  double x, fn, fn_exact, error;
  double x_array[2] = {.12345,.98765};
  //
  xmin = - 1.0;
  xmax =   1.0;
  n_intpol = 55;
  data_intpol(func, xmin, xmax, n_intpol, xj, fj);
  //
  for (k=0; k<=1; ++k)
    {
      x = x_array[k];
      printf("x = %22.15e\n", x);
      fn = interpolation(x, xj, fj, n_intpol);
      fn_exact = func(x);
      error = fabs(fn - fn_exact);
      printf("f(x)(interpolation) = %22.15e\n", fn);
      printf("f(x)(exact)         = %22.15e\n", fn_exact);
      printf("error               = %10.3e\n", error);
      printf("\n");
    }
  //
  return 0;
}