【HDU1724】Ellipse-自适应Simpson积分法

测试地址：Ellipse
题目大意：给定四个实数$a,b,l,r$，求椭圆$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$与半平面$x\ge l$和$x\le r$相交部分的面积。
做法：前段时间突击复习中考去了，结果还是考崩了2333…回来做一做比较好做的题找找手感~
这一道题需要使用自适应Simpson积分法来求解函数的定积分。
我们发现椭圆上下对称，所以题目所求的面积就是函数$f(x)=\sqrt{b^2(1-\frac{x^2}{a^2})}$在区间$[l,r]$上的积分的$2$倍。而自适应Simpson积分法是一种求函数定积分近似值的方法。它通过把曲边近似成过$(l,f(l)),(r,f(r)),(\frac{l+r}{2},f(\frac{l+r}{2}))$的抛物线，然后就可以简便地计算出这条抛物线在区间上的定积分，可以推导出这个定积分为：

$$\frac{r-l}{6}\Big(f(l)+f(r)+4f\Big(\frac{l+r}{2}\Big)\Big)$$
但是单纯这样算的话，精度误差是非常可观的，要保证整个值的误差在一定范围（设为 $eps$）内，则需要使用自适应Simpson积分法。思想上和二分法差不多，当算一个区间 $[a,b]$上的积分值时，先直接用公式计算区间 $[a,b]$上的积分值，然后再用公式计算区间 $[a,\frac{a+b}{2}]$和 $[\frac{a+b}{2},b]$上的积分值之和，如果这两个数值误差不超过一个精度范围，则认为没有细分下去的必要，否则递归求解两个半区间的积分值。需要注意的是，为了保证整体积分值的精度误差不超过 $eps$，所以当递归到第 $i$层时，允许的误差范围应该为$2^{-i}\times eps$。这样计算下去，最后的时间复杂度应该近似 $O(\log \frac{r-l}{eps})$，可能会更大，不过已经比较优秀了。
那么这一题就直接套用自适应Simpson积分法算所求的积分值即可，比较简单。
以下是本人代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
int T;
double A,B,l,r;

double f(double x)
{
  return sqrt(B*B*(1-x*x/(A*A)));
}

double calc(double a,double b)
{
  double mid=(a+b)/2;
  return (b-a)/6*(f(a)+f(b)+4*f(mid));
}

double simpson(double a,double b,double eps)
{
  double mid=(a+b)/2,s1=calc(a,mid),s2=calc(mid,b),s3=calc(a,b);
  if (fabs(s3-s1-s2)<=15*eps) return s1+s2+(s3-s1-s2)/15;
  else return simpson(a,mid,eps/2)+simpson(mid,b,eps/2);
}

int main()
{
  scanf("%d",&T);
  while(T--)
  {
    scanf("%lf%lf%lf%lf",&A,&B,&l,&r);
    printf("%.3lf\n",2*simpson(l,r,1e-6));
  }

  return 0;
}