補間法（ラグランジュ）

/****************************/
/* 補間法（ラグランジュ）   */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

double Lagrange(int, double *, double *, double);

int main()
{
	double x[3], x1, y[3], y1;
	int i1, n;
					// 線形補間
	n    = 1;
	x[0] = 0.0;
	x[1] = 2.0;
	y[0] = 0.0;
	y[1] = 4.0;

	printf("線形補間\n");
	x1 = 0.0;
	for (i1 = 0; i1 <= 8; i1++) {
		y1 = Lagrange(n, x, y, x1);
		printf("   x %f y %f\n", x1, y1);
		x1 += 0.25;
	}

					// ２次補間
	n    = 2;
	x[0] = 0.0;
	x[1] = 1.0;
	x[2] = 2.0;
	y[0] = 0.0;
	y[1] = 1.0;
	y[2] = 4.0;

	printf("２次補間\n");
	x1 = 0.0;
	for (i1 = 0; i1 <= 8; i1++) {
		y1 = Lagrange(n, x, y, x1);
		printf("   x %f y %f\n", x1, y1);
		x1 += 0.25;
	}

	return 0;
}

/**********************************/
/* ラグランジュ補間法             */
/*      n : 次数（n=1は線形補間） */
/*      x, y : (n+1)個の点の座標  */
/*      x1 : 補間したい点のｘ座標 */
/*      return : 補間結果         */
/**********************************/
double Lagrange(int n, double *x, double *y, double x1)
{
	double s1, s2, y1 = 0.0;
	int i1, i2;

	for (i1 = 0; i1 <= n; i1++) {
		s1 = 1.0;
		s2 = 1.0;
		for (i2 = 0; i2 <= n; i2++) {
			if (i2 != i1) {
				s1 *= (x1 - x[i2]);
				s2 *= (x[i1] - x[i2]);
			}
		}
		y1 += y[i1] * s1 / s2;
	}

	return y1;
}