// 龍貝格求積法.cpp : 定義控制臺應用程序的入口點。
//

#include "stdafx.h"

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

using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
	//float a=0,b=2,e=0.000001,result;
	//float a=3.14159265/2,b=3.14159265*3/4,e=0.000001,result;
	float a=2,b=3,e=0.000001,result;

	float f(float x);
	float romberg(float (*f)(float),float a,float b,float e);

	cout<<"【龍貝格求積法計算定積分】"<<endl<<endl;	
	result=romberg(f,a,b,e);
	cout<<"運算結果是："<<result<<endl;
	cout<<endl<<"按任意鍵繼續...";
	getch();
	return 0;
}
//龍貝格算法函數
float romberg(float (*f)(float),float a,float b,float e)
{
	float t1,t2,s1,s2,c1,c2,r1=0,r2=e;
	float s,h,x,k=1;
	h=b-a,t1=h*(f(a)+f(b))/2;
	while(r2-r1>=e || r1-r2>=e)
	{
		for(s=0,x=a+h/2;x<b;x+=h)
		s+=f(x);
		t2=(t1+h*s)/2;
		s2=t2+(t2-t1)/3;
		if(k==1)
		{
			k++,h/=2;
			t1=t2,s1=s2;
			continue;
		}c2=s2+(s2-s1)/15;
		if(k==2)
		{
			k++,h/=2;
			t1=t2,s1=s2,c1=c2;
			continue;
		}r2=c2+(c2-c1)/63;
		r1=r2,c1=c2;
		t1=t2,s1=s2;
		k++,h/=2;
	}
	return r2;
}

//積分函數
float f(float x)
{
	//return x*x*exp(-x*x);
	//return 1/tan(x);
	return 1/(x*x-1);

}