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

/***** 表3.4のガウス・ルジャンドル積分公式の数値解を求めるプログラム *****

	・まず，変数変換を行います．
		S=8/3integral^1_-1 (1-ξ^2)^(3/2)dξ

	・下記Nをルジャンドル多項式の次数として誤差を求めます．

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

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

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

#define N 3			//ルジャンドル多項式の次数
#define delta pow(10,-10)	//δ(二分法条件式(4.10))
#define epsilon pow(10,-10)	//ε(二分法条件式(4.10))
#define kmax 1000		//収束判定条件との比較回数の上限
#define Pi  acos(-1.)		//πの値

double f(double x){
	double y;
	y = 8./3.*pow((1-pow(x,2)), 1.5);
	return y;
}

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 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, wi;
	double xi[N+1], yi, I;
	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も根

	I = 0.;
	for(i=1; i<N+1; i++){
		wi = 0.;
		for(k=0; k<N; k++){
			lambdak = 2./(2.*(double)k+1.);
			wi += pow(P_k(k,xi[i]),2.)/lambdak;
							//式(2.96)部分
		}
		wi = 1./wi;				//式(2.96)
		yi = f(xi[i]);
		I += wi*yi;				//式(3.20)
	}

	printf("N=%4d, S=%19.12e\n", N, I);	
	return 0;
}