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

/********** 表2.1のラグランジュ補間の誤差ε_nを出力するプログラム ***********

	・実行すると，ラグランジュ補間の不等間隔
		補間区間の数N(textのNに対応)のデータ入力を求めます．

	・コマンドラインの第一引数にNを与えることでも
		即時に誤差を出力することができます．
		例）	>lagrange_n 12

	・L_2ノルムの2乗の計算には積分の台形公式(式(3.4))を用います．

****************************************************************************/

#define nn 100			//小区間中のきざみ数(L2nの精度に関わる)
#define Pi acos(-1.)		//πの値

double f(double x){		//f(x)の計算式
	double y;
	y = 1./(1.+25.*pow(x,2));
	return y;
}

/*** x=xaにおけるラグランジュ多項式の計算結果を返す関数 ***/
double Lagrange(int n, double xa, double *x, double *y){
	int i, j;
	double temp, f=0.;

	for(j=0; j<n+1; j++){
		temp = 1.;
		for(i=0; i<n+1; i++)
			if(i != j)
				temp *= ((xa-x[i])/(x[j]-x[i]));
		f += temp*y[j];
	}

	return f;
}

int main(int argc, char **argv){
	int i, j, N;
	double x, y, *xi, *yi, dx, h, L2n;

		/********* 入力処理(始) ***********/

	if(!(argc==1 || argc==2)){
		printf("引数のエラーです\n");
		return 1;
	}
	if(argc == 1){
		printf("補間点の数N+1に対応するNとして(表2.1のNとして)					Nを入力してください．\n N = ");
		scanf("%d",&N);
	}else{
		argv++;
		N = atoi(*argv);
	}
	if(N > 1000){
		printf("Nの値が大きすぎます．\n");
		return 2;
	}

		/********** 入力処理(終) ***********/

	/*** double型の配列xi[N+1]，yi[N+1]と同等のメモリ領域を確保 ***/
	xi  = malloc(sizeof(double)*(N+1));	
	yi  = malloc(sizeof(double)*(N+1));

	if(xi==NULL || yi==NULL){
		printf("can't malloc!!\n");
		return 4;
	}

	dx = 2./N;	//小区間の幅(補間点間の幅) dx
	h = dx/nn;	//小区間のきざみ幅 h(台形公式のh)

/******************************************************************************
			不等間隔にラグランジュ補間を適用
******************************************************************************/

	for(i=0; i<N+1; i++){
		xi[i] = cos(Pi*(N-i)/N);
		yi[i] = f(xi[i]);
	}

	L2n = 0.;

	/************************************************************
	下のコメント"//"を消去すると，一定間隔でf^N(x)を出力します．
	Excelなどでグラフ化するなら，>lagrange_n 12 > data.csv
	などとコマンドライン上でファイルにリダイレクトしてください．
	*************************************************************/
	for(j=0, x=-1; j<N*nn+1; j++, x+=h){

		/***** ラグランジュ補間を適用 *****/
		y = Lagrange(N, x, xi, yi);

		if(i==0 || i==N*nn)
			L2n += .5*h*pow((y-f(x)),2);
		else
			L2n += h*pow((y-f(x)),2);


//		if(fmod(j,nn)==0)
//			printf("%6.3lf,%11.3e\n",x, y);
	}

	printf("L2n=%11.3e\n",L2n);
	free(xi);			//メモリ領域の開放
	free(yi);
	return 0;
}