/*****************************************************************************/
/* page53 表3.2
/* ガウス・ルジャンドル積分公式のプログラム
/* f(x)=1.d0//sqrt(1.d0-x**2)
/* x=[-1,1]
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#define N 200

double pnx(int n, double x); /*ルジャンドル多項式の関数*/	
int zeroten(int n, double x[]);	/*Pn(x)のゼロ点を求める関数*/

int main(void){
	double x[N],w[N],y[N];
	double f,s,s0,z,pi=4.0*atan(1.0);
	int i,k,ni;

	s0 = pi;
	for(ni=50; ni<=N; ni=ni+50){
		zeroten(ni,x);
		for(i=0; i<ni; i++){
			f = 0.0;
			for(k=0; k<ni; k++){
				z = pnx(k,x[i]);
				f = f + (2.0*k+1.0)*z*z;
			}
			w[i] = 2.0/f;
			y[i] = 1.0/sqrt(1.0-x[i]*x[i]);
		}
		s=0.0;
		for(i=0; i<ni; i++){
			s =s+w[i]*y[i];
		}
		printf("|I-π|＝%20.16f\n", fabs(s0-s));
	}
	return 0;
}

int zeroten(int n, double x[]){
    int i,m;
    double a,b,c,fa,fb,fc;
    const double eps=1.e-12,pi=4.0*atan(1.0);
    
    m = n/2; 
    x[m] = 0.0;    
    for(i=0; i<m; i++){
		a = sin(pi*(double)(n-2*i-2)/(double)(2*n));
    	b = sin(pi*(double)(n-2*i  )/(double)(2*n));
    	fa = pnx(n,a);
    	fb = pnx(n,b);
    	if( fa*fb>0.0 ){
    		printf("Error a,b,i==%d\n",i);
    		return 0;
    	}
    	fc = 1.0;
    	while( fabs(fc)>eps && fabs(a-b)>eps ){
    		c = (a + b)/2.0;
            fc = pnx(n,c);
            if( fc*fa<0.0 ) b = c;
            else            a = c;
        }
        x[n-i-1] = c;
    	x[i] = -c;
    	
    }
	 return 1;
}

double pnx(int n, double x){
	double p0, p1, p2;
	int j;
	p0 = 1.0 ; 
	p1 = x;
	if( n == 0 )  return p0;
    if( n == 1 )  return p1;
    for(j=1; j<n; j++){
        p2 = ((double)(2*j+1)*x*p1 - (double)j*p0)/(double)(j+1);
        p0 = p1 ;
        p1 = p2;
    }
    return p2;

}