#include<iostream.h>//Romberg求積公式
#include<cmath>
double **BuildMatrix(int m,int n);
void FreeMatrix(double **A,int m);
double f(double x);
double sum(int k);
double ABS(double x);
double power(double x,int n);
double a,b,h;	
void main()
{
	cout<<"請輸入積分區間："<<endl;
	cout<<"積分下限a=";
	cin>>a;
	cout<<"積分上限b=";
	cin>>b;
	h=b-a;
	int n=2;
	while(1)
	{
		double **T=BuildMatrix(n,n);
		T[0][0]=h/2*(f(a)+f(b));
		for(int k=1;k<n;k++)
			T[k][0]=0.5*T[k-1][0]+h/power(2,k)*sum(k);
		for(int l=1;l<n;l++)
			T[0][l]=0;
		for(int i=1;i<n;i++)
		{
			for(int j=1;j<n;j++)
			{
				if(i<j)
					T[i][j]=0;
				else
					T[i][j]=(power(4,j)*T[i][j-1]-T[i-1][j-1])/(power(4,j)-1);
			}
		}
		if(ABS(T[n-1][n-1]-T[n-2][n-2])<=1e-5)
		{
			cout<<"所求定積分值為:"<<T[n-1][n-1];
			FreeMatrix(T,n);
			break;
		}
		else
		{
			FreeMatrix(T,n);
			n++;
		}
	}
}
double **BuildMatrix(int m,int n)
{
	double **A;
	A=new double *[m];
	for(int i=0;i<m;i++)
		A[i]=new double [n];
	return A;
}
void FreeMatrix(double **A,int m)
{
	for(int i=0;i<m;i++)
		delete[] A[i];
	delete[] A;
}
double sum(int k)
{
	double s=0;
	for(int i=1;i<power(2,k);i+=2)
		s+=f((b-a)*i/power(2,k));
	return s;
}
double ABS(double x)
{
	if(x>=0)
		x=x;
	else
		x=-x;
	return x;
}
double f(double x)
{
	double y;
	if(x==0)
		y=1;
	else
		y=sin(x)/x;
	return y;
}
double power(double x,int n)
{
	double val=1;
	while(n--)
		val*=x;
	return val;
}