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

Problem Description

Math is important!! Many students failed in 2+2’s mathematical test, so let’s AC this problem to mourn for our lost youth..
Look this sample picture:

A ellipses in the plane and center in point O. the L,R lines will be vertical through the X-axis. The problem is calculating the blue intersection area. But calculating the intersection area is dull, so I have turn to you, a talent of programmer. Your task is tell me the result of calculations.(defined PI=3.14159265 , The area of an ellipse A=PI*a*b )

---

Input

Input may contain multiple test cases. The first line is a positive integer N, denoting the number of test cases below. One case One line. The line will consist of a pair of integers a and b, denoting the ellipse equation , A pair of integers l and r, mean the L is (l, 0) and R is (r, 0). (-a <= l <= r <= a).

---

Output

For each case, output one line containing a float, the area of the intersection, accurate to three decimals after the decimal point.

---

Sample Input

2
2 1 -2 2
2 1 0 2

---

Sample Output

6.283
3.142

---

Solution

这是一道自适应辛普森积分的入门题，也是我写的第一道辛普森积分。

题目就是求两条线所夹的椭圆的面积。这是个标准的椭圆，其中心在原点。我好像看到了题目连椭圆面积公式都给了，然而并没有任何用，直接上辛普森！

根据椭圆的对称性，我们求$x$轴上方的面积然后乘以2就是答案。根据椭圆方程${x^2\over a^2}+{y^2\over b^2}=1$，构造函数$y=b\sqrt{1-({x\over a})^2}$，其中$y>0$，然后对此函数进行自适应辛普森积分就可以了。

写完程序过了样例后立马交，TLE，1000ms。结果发现是Eps取的太小了，取了1e-10，随着递归下去，Eps不断除以2，导致不断递归划分，做的次数很多就超时了。这题只要求保留3位小数，所以精度没有太大问题，将Eps改为1e-5立马就AC了，109ms。

---

Code

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdlib>

using namespace std;

int nG;
double a, b;

double f(double x){
    return b * sqrt(1.0 - (x * x) / (a * a));
}

double simpson(double L, double R){
    double mid = (L + R) / 2.0;
    return (f(L) + 4.0 * f(mid) + f(R)) * (R - L) / 6.0;
}

double integral(double L, double R, double Eps){
    double mid = (L + R) / 2.0;
    double ST = simpson(L, R), SL = simpson(L, mid), SR = simpson(mid, R);
    if(fabs(SL + SR - ST) <= 15.0 * Eps)  return SL + SR + (SL + SR - ST) / 15.0;
    //这里直接返回ST会WA
    return integral(L, mid, Eps/2.0) + integral(mid, R, Eps/2.0);
}

int main(){

    scanf("%d", &nG);

    double L, R, Eps = 1e-5;//Eps不能取得太小
    while(nG --){
      scanf("%lf%lf%lf%lf", &a, &b, &L, &R);
      printf("%.3lf\n", 2.0 * integral(L, R, Eps));
    }

    return 0;
}

---