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

/******* 表2.1のラグランジュ補間の誤差ε_1,2,3,4を出力するプログラム ********

	・実行すると，ラグランジュ多項式の次数n(textのεの添字に対応)と
		補間区間の数N(textのNに対応)のデータ入力を求めます．

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

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

	・補間点で重複計算がありますが，
			誤差はないので無視します．

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

#define nn 100			//小区間中のきざみ数(L2nの精度に関わる)

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, imax, jmax, N, n;
	double x, y, *xi, *yi, dx, h, L2n=0.;
	double Min;		/*[-1,1]を分割した小区間の最小値(変数)*/
	double Max;		/*[-1,1]を分割した小区間の最大値(変数)*/

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

	if(!(argc==1 || argc==3)){
		printf("引数のエラーです\n");
		return 1;
	}
	if(argc == 1){
		printf("ラグランジュ多項式の次数を入力してください．					ε_1なら1，ε_2なら2などとしてください．\n n = ");
		scanf("%d",&n);
		printf("補間点の数N+1に対応するNとして(表2.1のNとして)					Nを入力してください．\n N = ");
		scanf("%d",&N);
	}else{
		argv++;
		n = atoi(*argv);
		argv++;
		N = atoi(*argv);
	}
	if(n > 1000){
		printf("次数nの値が大きすぎます．\n");
		return 2;
	}
	if(fmod(N, n) != 0){
		printf("入力エラーです．\n");
		printf(" N = %dとするなら\n n=",N);
		for(i=1; i<N+1; i++){
			if(fmod(N, i) == 0)
				printf("%d ",i);
		}
		printf("の中から選んで入力してください．\n");
		return 3;
	}

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

	/*** 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)

/******************************************************************************
			小区間ごとにラグランジュ補間を適用
******************************************************************************/
	imax = (int)(N/n);			//小区間の数

	for(i=0; i<imax; i++){
		Min = -1.+(double)i*n*dx;	//小区間の最小値
		Max = Min+n*dx;			//小区間の最大値

		for(j=0; j<n+1; j++){
			xi[j] = Min+(double)j*dx;	//適用小区間内のxi
			yi[j] = f(xi[j]);		//適用小区間内のyi
		}

		jmax = (int)(n*nn);			//int切捨て

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

			/***** 小区間にn次ラグランジュ補間を適用 *****/
			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,(int)(N/10))==0)
//				printf("%6.3lf,%11.3e\n",x, y);
		}
	}

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