Hdu 1724 Ellipse (自适应辛普森积分法)

辛普森积分公式：http://zh.wikipedia.org/zh-tw/%E8%BE%9B%E6%99%AE%E6%A3%AE%E7%A9%8D%E5%88%86%E6%B3%95

代码参考：http://blog.csdn.net/oceanlight/article/details/12259983

以下原理摘自：http://wenku.baidu.com/link?url=5Mkew3kQ34hji4mkY_zouYLyuDtq7-zdjk8z7BDogLKSoCS3hS4fGH_0MYTm0u0Pc3YUAEH9h6LLCNEyYIyqrBiyluHYhiEttxMSFNC4qH7

自适应积分法是一种比较经济而且快速的求积分的方法。他能自动地在被积函数变化剧烈的区域增多节点，
而在被积函数变化平缓的地方减少节点。因此它是一种不均匀区间的积分方法。
按照子区间上的积分方式它可以分为自适应辛普森积分法和自适应梯形积分法。
通常是采用自适应辛普森积分法作为子区间的积分方式。


自适应积分法的基本步骤如下：
(1) 将积分区间[a,b]分成两个相等的1级子区间[a,a+0.5h]和[a+0.5h,a+h],且h=b-a;
(2) 在上述两个1级子区间上用辛普森积分得到积分I_{a,a+0.5h}^{(1)}和I_{a+0.5h,a+h}^{(1)};
(3) 将子区间[a,a+0.5h]分成两个相等的2级子区间[a,a+0.5^2*h]和[a+0.5^2*h,a+0.5h];
(4)  采用辛普森积分计算得到:I_{a,a+0.5h}^{(2)}=I_{a,a+0.5^2*h}^{(1)}+I_{a+0.5^2*h,a+0.5h}^{(1)};
(5) 比较I_{a,a+0.5h}^{(2)}和I_{a,a+0.5h}^{(1)}, 如果|I_{a,a+0.5h}^{(1)} - I_{a,a+0.5h}^{(2)}|<10*0.5*epsi，
其中epsi为整体积分所需要的精度，则认为子区间[a,a+0.5h]上的积分I_{a,a+0.5h}^{(1)}
已达到所需精度, 不需要再细分;
否则就需要再细分, 对每个2级子区间做同样的判断.1级子区间[a+0.5h,a+h]的操作过程完全与上面相同

#include <cstdio>
#include <cmath>

const double eps = 1e-6;  //积分精度

double a,b , l, r;
// simpson公式用到的函数
double F (double x){//也就是被积函数,本题是椭圆
     return sqrt(b*b*(1-x*x/(a*a)));
}

// 三点simpson法. 这里要求F是一个全局函数
double simpson (double a,double b){
	double c =  a+(b-a)/2;
	return (F(a) + 4*F(c) + F(b))*(b-a)/6;
}
// 自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点simpson值A
double asr(double a , double b ,double eps ,double A){
	double c = a+ (b-a)/2;
	double L = simpson(a,c) ,R = simpson(c,b);
	if (fabs(A-L-R)<=15*eps) return L + R +(A-L-R)/15;
	return asr(a,c,eps/2,L) + asr(c,b,eps/2,R);
}
// 自适应Simpson公式（主过程）  
double asr(double a, double b, double eps) {  
	return asr(a, b, eps, simpson(a, b));  
}

int main ()
{
	int T;
	scanf("%d",&T);
	while (T--)
	{
		scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
		double ans = asr(l,r,eps);
		ans *= 2;//积分只计算了x轴上半部分
		printf("%.3f\n",ans);
	}
	return 0;
}