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

/************** ルジャンドル補間の誤差を出力するプログラム *****************

	・下記Nをルジャンドル多項式の次数として
	
		y=1/(1+25x^2)
		
		のL_2ノルムの2乗を計算します．

	・Nを変更する場合はコンパイルしなおしてください．

	・多項式の0点を求めるのに，4.2節の二分法を用います．

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

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

#define N 20			//ルジャンドル多項式の次数
#define delta pow(10,-10)	//δ(二分法条件式(4.10))
#define epsilon pow(10,-10)	//ε(二分法条件式(4.10))
#define kmax 1000		//収束判定条件との比較回数の上限
#define Pi acos(-1.)		//πの値
#define nn 1000			//L2nの精度に関わる[-1,1]のきざみ数
#define h 2./nn			//台形公式(式(3.4))のh

double f(double x){
	return (1./(1.+25.*pow(x,2)));
}

double P_k(int k, double x){		//式(2.87)
	int n;
	double p0, p1, p2;
	p0 = 1;
	p1 = x;
	if( k==0 ){
		return p0;
	}else{
		for(n=1; n<k; n++){
			p2 = ((2*(double)n+1)*x*p1-(double)n*p0)/(n+1);
			p0 = p1;
			p1 = p2;
		}
	}
	return p1;
}

double fPN(double x, double ck[]){		//式(2.94)
	int k;
	double f=0.;
	for(k=0; k<N; k++)
		f += ck[k]*P_k(k,x);
	return f;
}

/***** "区間"を引数にとり，根を返す関数 *****/
double NIBUN(double a, double b){
	int k;
	double c;
	for(k=1; k<kmax+1; k++){
		c = (a+b)/2;
		if(fabs( P_k(N,c) )<delta || fabs(a-b)<epsilon)
			return c;	//根を返す
		if(P_k(N,c)*P_k(N,a)<0){
//				a=a;
				b=c;
		}else{
				a=c;
//				b=b;
		}
//		printf("| %4d回 | %13.5e | %9.6lf | %13.5e |\n"			//			, k, fabs(a-b), c, P_k(N,c));
	}
	printf("error!!\n");
	exit(1);			//エラーによる強制終了
}

int main(void){
	int i, m, k, kk;
	double a, b;
	double lambdak, lambdakk, wi, L2n, x;
	double xi[N+1], ck[N];

	if( fmod(N,2)==0 )		//次数が偶数のとき
			m = (int)(N/2);
	else{				//次数が奇数のとき
			m = (int)((N-1)/2);
	}

	/***** N次ルジャンドル多項式の根すべてをxi[]に登録 *****/
	for(i=1; i<m+1; i++){

		//p.70の区間[a,b]
		a = sin((Pi*(N-2*i))/(2*N));
		b = sin((Pi*(N+2-2*i))/(2*N));
		xi[i] = NIBUN(a,b);	//根を代入
		xi[i+m] = -xi[i];	//-xiも根
//		printf("\n x[%2d]=%8.7lf\n", i, xi[i]);
	}

	if( fmod(N,2)==1 )	xi[N] = 0;	//奇数のときは0も根

	/***** ck[]を求める *****/
	for(k=0; k<N; k++){
		lambdak = 2/(2*(double)k+1);
		ck[k] = 0.;
		for(i=1; i<N+1; i++){
			wi = 0.;
			for(kk=0; kk<N; kk++){
				lambdakk = 2/(2*(double)kk+1);
				wi += pow(P_k(kk,xi[i]),2)/lambdakk;
								//式(2.96)部分
			}
			wi = 1/wi;				//式(2.96)
			ck[k] += wi*P_k(k,xi[i])*f(xi[i]);	//式(2.94)部分
		}
		ck[k] /= lambdak;				//式(2.94)
	}

/***** L_2ノルムの2乗(式(2.59))を積分の台形公式(式(3.4))を用いて求める *****/
	L2n = 0.;

	/************************************************************
	下のコメント"//"を消去すると，一定間隔でf^P_N(x)を出力します．
	Excelなどでグラフ化するなら，>legendre > data.csv　などと
	コマンドライン上でファイルにリダイレクトしてください．
	*************************************************************/
	for(i=0; i<nn+1; i++){
		x = -1.+h*(double)i;
//		if(fmod(i,(int)(nn/10))==0)					//			printf("%6.3lf,%11.3e\n", x, fPN(x,ck));

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

	}

	printf("N=%4d, L2n=%16.10e\n", N, L2n);	
	return 0;
}