#include<stdio.h>
#include<stdlib.h>
#include <math.h>
#define f(x) (4.0/(1+(x)*(x)))//函數定義
int MAX_N=0;//最大的迭代次數
double Romberg(double a,double b,double eps,int *k);
void main(){
	int k=0;
	double r,a,b,eps;
	printf("Please input the maximum iteration number MAX_N:");
	scanf("%d",&MAX_N);
	printf("Please input the interval [a,b]:");
	scanf("%lf,%lf",&a,&b);
	printf("Please input the precision:");
	scanf("%lf",&eps);
	printf("The list of calculational process:\n");
	printf("*******************************\n");
	printf("Num T           S            C            R\n");
	r=Romberg(a,b,eps,&k);
	printf("*******************************\n");
	printf("The result of Romberg is:%.10lf. And k is:%d.\n",r,k);
}
double Romberg(double a,double b,double eps,int *k)
{
	double h=1,summ,*p;
	int m=1;
	int i,j;
	double p1[4]={0},p2[4]={0};
	p1[0]=h*(f(a)+f(b))*1.0/2;
	printf("1:");
	for (j=0;j<4;j++)
		printf("%.10lf ",p1[j]);
	printf("\n");
	for (i=2;i<=MAX_N;i++)
	{
		m=2*m;
		h=h/2.0;
		p=malloc(sizeof(double)*(m+1));
		if (p==NULL)
	 	{
			printf("Memory is NULL!\n");
			exit(1);
		}
		for (j=0;j<=m;j++)
			p[j]=a+h*j;
		summ=0;
		for (j=0;j<m;j++)
			summ+=f((p[j]+p[j+1])/2.0);
		p2[0]=0.5*p1[0]+h/2.0*summ;        //T
		p2[1]=(4*p2[0]-p1[0])/3.0;		   //S
		if (i>=3)
			p2[2]=(16*p2[1]-p1[1])/15.0;   //C
		if (i>=4)
			p2[3]=(64*p2[2]-p1[2])/63.0;   //R	
		if ((i>=5)&&(fabs(p1[3]-p2[3])<=eps))
		{
			free(p);
			break;
		}
		printf("%d:",i);
		for (j=0;j<4;j++)
		{
			printf("%.10lf ",p2[j]);
			p1[j]=p2[j];
		}
		printf("\n");
		free(p);
	}
	printf("%d:",i);
	for (j=0;j<4;j++)
		printf("%.10lf ",p2[j]);
	printf("\n");
	if (i<=MAX_N)
		*k=i;
	else
		*k=-1;
	return p2[3];
}