hdoj-1724 Ellipse(自适应Simpson积分)

浙大月赛的时候卡在了一道积分题目上，被积函数无法求出原函数。

连接：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3898

问了学长方法，学习了刘汝佳的白书《算法竞赛入门经典训练指南》p169附近介绍的方法。

自适应Simpson积分也是在被积函数上取微元然后积分，普通的积分取点越多，精度越高，

但是时间也越长，在设定一个精度eps之后，自适应Simpson积分可以根据情况在容易近似成梯形的地方少划分，

在不容易拟合和梯形的地方多取点，这是一个递归的过程，时间和精度都是比较好的。

先拿hdoj上这道题练练手吧。

题意：给出一个椭圆的参数，求x=[l,r]之间椭圆被直线划分出来的面积

注意l，r正负时的情况

（请教数学老师后得知还有一种精度更高的Romberg算法）

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;
double PI=3.14159265;
double a,b;
double e=1e-8;

//被积函数的表达式，这里是x轴上方的椭圆函数
double F(double x)
{
    return sqrt(fabs(a*a-x*x))*b/a;
}

double simpson(double a,double b)
{
    double c=a+(b-a)/2;
    return (F(a)+4*F(c)+F(b))*(b-a)/6;
}

double asr(double a,double b, double eps,double A)
{
    double c=a+(b-a)/2;
    double L=simpson(a,c);
    double R=simpson(c,b);
    if(fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15.0;
    return asr(a,c,eps/2,L)+asr(c,b,eps/2,R);
}

double asr(double a,double b,double eps)
{
    return asr(a,b,eps,simpson(a,b));
}

void work()
{
    double ans,l,r;
    scanf("%lf%lf%lf%lf",&a,&b,&l,&r);
    ans=0;
    if(l<0) ans=asr(0,fabs(l),e);
    else ans=-asr(0,fabs(l),e);
    if(r<0) ans-=asr(0,fabs(r),e);
    else ans+=asr(0,fabs(r),e);
    ans*=2.0;
    printf("%0.3lf\n",ans);

}

int main()
{
    int T;
    double a,b;
    scanf("%d",&T);
    while(T--)
        work();
    return 0;
}