龙贝格求积算法

函数f(x)=sinx/x,区间[0,1]
算法流程图如图：

/*
-------------Analysis----------------------
*龙贝格算法求积分
*1.函数f(x)=sinx/x,区间[0,1]
*2.输入a,b,e;其中a,b为区间[a,b]，e为精度
*3.T为复化梯形公式，S为辛甫生公式，C为复化科特斯公式
*4.h为步长 x为求积节点,k为加速次数
*5.输出R2
-----------------------------------
*/
#include<iostream>
#include<math.h>
using namespace std;
//get f(x) according to a specified x
double f(double x)
{
	if(0==x)
		return 1;
	return sin(x)/x;
}
int main()
{
	double a,b;
	double e;
	cout<<"Please give values of section [a,b] and precision e"<<endl;
	cout<<"------------------------------------------"<<endl;
	cout<<"Please input pre-section a:"<<endl;
	cin>>a;
	cout<<"Please input suf-section b:"<<endl;
	cin>>b;
	cout<<"Please input precision e:"<<endl;
	cin>>e;
	cout<<"Calculating.........plese wait..."<<endl;
	
	unsigned k;
	double h;
	double x;
	double S;
	double T1,T2;
	double S2,S1;
	double C2,C1;
	double R2,R1;

	h=b-a;
	T1=(h/2)*(f(a)+f(b));
	k=1;
	S=0;
	x=a+h/2;
	S=S+f(x);
	x=x+h;
	while(x<b)
	{
		S=S+f(x);
		x=x+h;
	}
	T2=T1/2+(h/2)*S;	
	S2=T2+(1/3)*(T2-T1);
	while(1==k)
	{
		k=k+1;
		h=h/2;
		T1=T2;
		S1=S2;
		S=0;
		x=a+h/2;
		S=S+f(x);
		x=x+h;
		while(x<b)
		{
			S=S+f(x);
			x=x+h;
		}
		T2=T1/2+(h/2)*S;		
		S2=T2+(1/3)*(T2-T1);
	}

	C2=S2+(1/15)*(S2-S1);
	while(2==k)
	{
		C1=C2;
		k=k+1;
		h=h/2;
		T1=T2;
		S1=S2;
		
		S=0;
		x=a+h/2;

		//one cycle
		S=S+f(x);
		x=x+h;
		while(x<b)
		{
			S=S+f(x);
			x=x+h;
		}
		T2=T1/2+(h/2)*S;
		
		S2=T2+(1/3)*(T2-T1);
		while(1==k)
		{
			k=k+1;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
			x=a+h/2;
			S=S+f(x);
			x=x+h;
			while(x<b)
			{
				S=S+f(x);
				x=x+h;
			}
			T2=T1/2+(h/2)*S;	
			S2=T2+(1/3)*(T2-T1);
		}
		C2=S2+(1/15)*(S2-S1);
	}
	R2=C2+(1/63)*(C2-C1);
	while(3==k)
	{
		R1=R2;
		C2=C1;
		k=k+1;
		h=h/2;
		T1=T2;
		S1=S2;	
		S=0;
		x=a+h/2;
		S=S+f(x);
		x=x+h;
		while(x<b)
		{
			S=S+f(x);
			x=x+h;
		}
		T2=T1/2+(h/2)*S;
		S2=T2+(1/3)*(T2-T1);
		while(1==k)
		{
			k=k+1;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
			x=a+h/2;
			S=S+f(x);
			x=x+h;
			while(x<b)
			{
				S=S+f(x);
				x=x+h;
			}
			T2=T1/2+(h/2)*S;	
			S2=T2+(1/3)*(T2-T1);
		}
		C2=S2+(1/15)*(S2-S1);
		while(2==k)
		{
			C1=C2;
			k=k+1;
			h=h/2;
			T1=T2;
			S1=S2;			
			S=0;
			x=a+h/2;
			S=S+f(x);
			x=x+h;
			while(x<b)
			{
				S=S+f(x);
				x=x+h;
			}
			T2=T1/2+(h/2)*S;
			
			double S2;
			S2=T2+(1/3)*(T2-T1);
			while(1==k)
			{
				k=k+1;
				h=h/2;
				T1=T2;
				S1=S2;
				S=0;
				x=a+h/2;
				S=S+f(x);
				x=x+h;
				while(x<b)
				{
					S=S+f(x);
					x=x+h;
				}
				T2=T1/2+(h/2)*S;
				
				S2=T2+(1/3)*(T2-T1);
			}
		
			C2=S2+(1/15)*(S2-S1);	
		}
		R2=C2+(1/63)*(C2-C1);
	}
	while(fabs(R2-R1)>e)
	{
		R1=R2;
		C2=C1;
		k=k+1;
		h=h/2;
		T1=T2;
		S1=S2;
		S=0;
		x=a+h/2;
		S=S+f(x);
		x=x+h;
		while(x<b)
		{
			S=S+f(x);
			x=x+h;
		}
		T2=T1/2+(h/2)*S;
		
		S2=T2+(1/3)*(T2-T1);
		while(1==k)
		{
			k=k+1;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
			x=a+h/2;
			S=S+f(x);
			x=x+h;
			while(x<b)
			{
				S=S+f(x);
				x=x+h;
			}
			T2=T1/2+(h/2)*S;	
			S2=T2+(1/3)*(T2-T1);
		}
		C2=S2+(1/15)*(S2-S1);
		while(2==k)
		{
			C1=C2;
			k=k+1;
			h=h/2;
			T1=T2;
			S1=S2;
			S=0;
			x=a+h/2;
			S=S+f(x);
			x=x+h;
			while(x<b)
			{
				S=S+f(x);
				x=x+h;
			}
			T2=T1/2+(h/2)*S;
			double S2;
			S2=T2+(1/3)*(T2-T1);
			while(1==k)
			{
				k=k+1;
				h=h/2;
				T1=T2;
				S1=S2;
				S=0;
				x=a+h/2;
				S=S+f(x);
				x=x+h;
				while(x<b)
				{
					S=S+f(x);
					x=x+h;
				}
				T2=T1/2+(h/2)*S;
				S2=T2+(1/3)*(T2-T1);
			}
			C2=S2+(1/15)*(S2-S1);
		}
		R2=C2+(1/63)*(C2-C1);
	}
	cout<<"R2:"<<R2<<endl;
	cout<<"Hello Boker.."<<endl;
	system("pause");
	return 0;
}

运行结果：

虽然能运行出来，但比较繁琐，希望能有好的数据结构。