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

double e=0.00000005;

double f1(double x)
{
	return (-2)/(x*x-1);	
}

double f2(double x)
{
	return 4/(x*x+1);	
}

double f3(double x)
{
	return pow(3,x);	
}

double f4(double x)
{
	return x*exp(x);	
}

void Tixing(int a,int b,double M,double (*f)(double x))
{
	double w,s;
	w=fabs((-1)*(b-a)*(b-a)*(b-a)*M/12/e);
	s=sqrt(w);
	int n;
	n=int(ceil(s))+1;
	double *x=new double[n+1];
	double h=double(b-a)/n;
	x[0]=a;x[n]=b;
	for(int i=1;i<n;i++)
	{
		x[i]=x[i-1]+h;
		
	}
	double *y=new double[n+1];
	for(int j=0;j<n+1;j++)
	{
		y[j]=f(x[j]);
		
	}
	
	double T=h/2*(y[0]+y[n]);
	for(i=1;i<n;i++)
	{
		T=T+h*y[i];

	}
	
	cout<<"   "<<"用復(fù)化梯形公式求得："<<"T="<<setprecision(10)<<T<<endl;
	delete[]x;
	delete[]y;
}

void Simpson(int a,int b,double M,double (*f)(double x))
{
	double w,s;
	w=sqrt(fabs((-1)*(b-a)*(b-a)*(b-a)*(b-a)*(b-a)*M/180/16/e));
	s=sqrt(w);
	int n,m;
	n=int(ceil(s)+1);
	m=2*n;
	//cout<<n<<endl;
	double h=double(b-a)/n;
	double *x=new double[m+1];
	x[0]=a;x[m]=b;
	for(int i=1;i<m;i++)
	{
		x[i]=x[i-1]+h/2;
		//cout<<x[i]<<"\t";
	}
	double *y=new double[m+1];
	for(int j=0;j<m+1;j++)
	{
		y[j]=f(x[j]);
		//cout<<y[j]<<"\t";
	}
	double S=h/6*(y[0]+y[m]);
	for(i=1;i<m;i++)
	{
		if(i%2==1)
			S=S+4*h/6*y[i];
		else
			S=S+2*h/6*y[i];
		//cout<<S<<"\t";
	}
	cout<<"   "<<"用復(fù)化Simpson公式求得："<<"S="<<setprecision(10)<<S<<endl;
	delete[]x;
	delete[]y;
}

void Romberg(int a,int b,double (*f)(double x))
{
	int n=1;
	double T[10];
	T[0]=(b-a)*(f(a)+f(b))/2;
	for (int i=1;i<10;i++)
	{
		double h=(b-a)/pow(2,i);
		T[i]=h*(f(a)+f(b))/2;
		for (int j=1;j<pow(2,i);j++)
		{
			T[i]+=h*f(a+h*j);
		}

	}

	double S[9];
	S[0]=4*T[1]/3-T[0]/3;
	for (i=1;i<9;i++)
	{
		S[i]=4*T[i+1]/3-T[i]/3;
	}
	
	double C[8];
	C[0]=16*S[1]/15-S[0]/15;
	for (i=1;i<8;i++)
	{
		C[i]=16*S[i+1]/15-S[i]/15;
	}
	double R[7];
	R[0]=64*C[1]/63-C[0]/63;
	for (i=1;i<7;i++)
	{
		R[i]=64*C[i+1]/63-C[i]/63;
		if (fabs(R[i]-R[i-1])<e)
			break;
	}
	cout<<"   "<<"用Romberg公式求得"<<"R="<<setprecision(10)<<R[i]<<endl;
}

void main()
{
	cout<<"相對(duì)誤差限為1*10^(-7)/2,分別用復(fù)化梯形公式、復(fù)化Simpson公式、Romberg公式計(jì)算："<<endl;
	cout<<endl;
	cout<<"用三種方法解(-2)/(x*x-1)從2到3的積分：（即求ln(2/3)，精確解為-0.405465108）"<<endl;
	Tixing(2,3,1,f1);
	Simpson(2,3,5,f1);
	Romberg(2,3,f1);
	cout<<endl;
	cout<<"用三種方法解4/(x*x+1)從0到1的積分：（即求圓周率π，精確解為3.141592654）"<<endl;
	Tixing(0,1,5,f2);
	Simpson(0,1,1,f2);
	Romberg(0,1,f2);
	cout<<endl;
	cout<<"用三種方法解3的x次方從0到1的積分：（即求2/ln3，精確解為1.820478453）"<<endl;
	Tixing(0,1,2.5,f3);
	Simpson(0,1,2,f3);
	Romberg(0,1,f3);
	cout<<endl;
	cout<<"用三種方法解x乘e的x次方從1到2的積分：（即求e的平方，精確解為7.389056099）"<<endl;
	Tixing(1,2,18,f4);
	Simpson(1,2,20,f4);
	Romberg(1,2,f4);
}