#include<iostream.h>
#include<math.h>

double f(double x)
{
	return 1.0/(1+x);
}

void Romberg(double a, double b, double e)
{
	double h = b - a;
	int k = 0, j;
	double T[10][10];
	T[0][0] = ( f(a) + f(b) ) * h / 2;
	cout<<'\t'<<"k\t"<<"T\t"<<"S\t"<<"C\t"<<"R\n";
	cout<<'\t'<<k<<'\t'<<T[0][0]<<'\n';
	for(k = 1; ; k++)
	{
		double temp = 0;
		for(int i = 1; i<=pow(2,k-1); i++)
		{
			temp += f(a + (i-1.0/2)*h );
		}
		T[0][k] = 1.0/2 * ( T[0][k-1] + h*temp );
		cout<<'\t'<<k<<'\t'<<T[0][k];
		for(j=1; j<=3&&j<=k; j++)
		{
			T[j][k-j] = ( pow(4,j)*T[j-1][k-j+1] - T[j-1][k-j] ) / ( pow(4,j) - 1 );
			cout<<'\t'<<T[j][k-j];
		}
		cout<<'\n';
		if( fabs(T[j-1][k-j+1]-T[j-2][k-j+2] )<e )
		{
			cout<<"結果為:"<<T[j-1][k-j+1]<<'\n';
			break;
		}
		else
		{
			h = h/2;
		}
	}	
}

void main(void)
{
	double a, b, e;
	cout<<"輸入上限:";
	cin>>b;
	cout<<"輸入下限:";
	cin>>a;
	cout<<"輸入精度:";
	cin>>e;
	Romberg(a, b, e);
}