/*
 * runge.c --- 等間隔標本点の補間多項式はポシャるという Runge の現象
 * 参考: 森正武, 数値解析, 共立出版 (1973, 第2版 2002).
 * gcc runge.c ; ./a.out > runge.data
 * gnuplot> f(x)=1/(1+25*x*x)
 * gnuplot> plot [-1:1] [-1:1] "runge.data" with linespoints, f(x)
 * gnuplot> plot [-1:1] [-1:10] "runge.data" with linespoints, f(x)
 *
 * ここでは Lagrange 補間多項式として計算している。
 */

#include <stdio.h>
#include <stdlib.h>

/* [-1,1] での等間隔標本点がうまく行かないことで有名な関数 */
double f(double x)
{
  return 1.0 / (1.0 + 25.0 * x * x);
}

/* Lagrange 補間係数 */
double L(double x, int k, int N, double xv[])
{
  int j;
  double t = 1;
  for (j = -N; j <= N; j++)
    if (j != k)
      t *= (x - xv[j+N]) / (xv[k+N] - xv[j+N]);
  return t;
}

/* Lagrange 補間公式 */
double fn(double x, int N, double xv[], double fv[])
{
  int k;
  double s = 0;
  for (k = -N; k <= N; k++)
    s += fv[k+N] * L(x, k, N, xv);
  return s;
}

int main(void)
{
  int j, N, nn;
  double h;
  double *xv, *fv;
  N = 10;
  xv = malloc(sizeof(double) * (2 * N + 1)); // エラーチェックさぼり
  fv = malloc(sizeof(double) * (2 * N + 1)); // 同上
  h = 1.0 / N;
  for (j = -N; j <= N; j++) {
    xv[j + N] = j * h;
    fv[j + N] = f(xv[j + N]);
  }  
  nn = 200;
  h = 2.0 / nn;
  printf("%g %g\n", -1.0, fn(-1.0, N, xv, fv));
  for (j = 1; j <= nn; j++)
    printf("%g %g\n", -1.0 + j * h, fn(-1.0 + j * h, N, xv, fv));
}