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

/****************** 表3.5の誤差を出力するプログラム *************************

	・実行すると，誤差の絶対値がマクロのepsよりも
		小さくなったときにおけるL，hとその誤差 |I-rPi|を出力します．

	・実行結果の最後部から一部を下に示します．
		(マクロの値は，Lstep:1，eps:pow(10,-8)，cntl:3)
				・
				・
			L=7.000000, h=0.700,								 |I-rPi|=6.33905e-09
			L=7.000000, h=0.667,								 |I-rPi|=8.04261e-10
			L=7.000000, h=0.636,								 |I-rPi|=9.22641e-11
			L=6.000000, h=0.706,								 |I-rPi|=8.85542e-09
			L=6.000000, h=0.667,								 |I-rPi|=8.04260e-10
			L=6.000000, h=0.632,								 |I-rPi|=6.36878e-11
			L=5.000000, h=0.667,								 |I-rPi|=8.13534e-10
			L=5.000000, h=0.625,								 |I-rPi|=2.90233e-11
			L=5.000000, h=0.588,								 |I-rPi|=9.65417e-12

	・これより，厳密ではありませんが，L>=5，h<=0.7の条件が得られます．

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

#define Pi acos(-1.)		//πの値
#define rPi sqrt(Pi)		//ルートπの値
#define iniL 10			//Lの初期値(textに対応)
#define Lstep 1			//Lのステップ(初期値からの減算)
#define eps pow(10,-8)		//許容誤差
#define nmax 300		//区間[-L, L]の分割数の最大値
#define cntl 3			//L=constにおけるhの出力数(上記の結果を参照)

double f(double x){
	double y;
	y = exp(-pow(x,2));
	return y;
}

int main(void){
	int n, i, cnt;
	double I, L, h, xn;

	for(L=iniL; 0<L; L-=Lstep){
		cnt = 0;
		for(n=1; n<nmax+1; n++){
			h = 2*L/n;
			I = 0.;
			for(i=0; i<n+1; i++){
				xn = -L+i*h;
				if(i==0 || i==n)
					I += .5*h*f(xn);
				else
					I += h*f(xn);
			}
			if(fabs(I-rPi) < eps){
				cnt++;
				printf("L=%3lf, h=%5.3lf,								 |I-rPi|=%11.5e\n",									L, h, fabs(I-rPi));
			}
			if(cnt == cntl)
				goto label1;
		}
		label1:
	}
	return 0;
}