/*用龍貝格積分計算P140的例6.8*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 20               /*最大二份次數*/
#define EPS 0.00001        /*精度控制*/
double f(double x)
{
	if(x==0)
		return 1;
	else
		return(sin(x)/x);
}
void Romberg(float a,float b)
{
	double T[N][4],h,p;
	h=fabs(b-a);
	T[0][0]=h*(f(a)+f(b))/2;
	p=EPS+1;
	int k=1;
	while(k<N&&p>EPS)
	{
		double sum=0;
		for(int i=1;i<=(pow(2,k-1));i++)
			sum=sum+f(a+(i-0.5)*h);
		T[k][0]=(T[k-1][0]+h*sum)/2;
		for(int m=1;m<=k;m++)
		{
			if(m>3)
				break;
			T[k][m]=(pow(4,m)*T[k][m-1]-T[k-1][m-1])/(pow(4,m)-1);
		}
		if(k<=3)
			p=fabs(T[k][k-1]-T[k-1][k-1]);
		else
			p=fabs(T[k][3]-T[k-1][3]);
		k++;
		h=h/2;
	}
	if(k>N)
		printf("方法失敗!\n");
	else
	{
		printf("\tT\t\tS\t\tC\t\tR\n");
		for(int i=0;i<k;i++)
		{
			printf("k=%d\t",i);
			for(int j=0;j<4;j++)
			{
				if(i>=j)
					printf("%-10.9f\t",T[i][j]);
			}
			printf("\n");
		}
	}
}
void main()
{
	float a,b;
	printf("請輸入積分區間:\n");
	scanf("%f%f",&a,&b);
	Romberg(a,b);
}