Ellipse

题目描述

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 at this sample picture:

An ellipse in a plane centered on 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 to turn to you, a talent of programmer. Your task is telling me the result of calculation. (defined $\pi=3.14159265$, the area of an ellipse $A=\pi ab$)

题意概述

给定$a, b, l, r$，求椭圆${x^2 \over a^2}+{y^2 \over b^2}=1$在直线$x=l$与$x=r$之间的面积。
数据范围：$-a \le l \le r \le a$。

算法分析

这是一个不规则图形，无法直接计算，但我们可以用积分来求面积。设$f(n)$表示直线$x=n$与椭圆的两个交点之间的距离。
根据辛普森公式
$$
\int_a^b f(x) \, {\rm d}x \approx {b-a \over 6} \left(f(a)+4f\left({a+b \over 2}\right)+f(b)\right)
$$
但这样的精度并不是很高。考虑二分，令
$$g(l, r)=\int_l^r f(x) \, {\rm d}x, \; h(l, r)={r-l \over 6} \left(f(l)+4f\left({l+r \over 2}\right)+f(r)\right)$$
当$h(l, r)$与$h\left(l, {l+r \over 2}\right)+h\left({l+r \over 2}, r\right)$的差足够小时，令
$$g(l, r)=h(l, r)$$
否则
$$g(l, r)=g\left(l, {l+r \over 2}\right)+g\left({l+r \over 2}, r\right)$$
这样就可以达到精度要求。

代码

/*
 * Today is what happened to yesterday.
 */

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

static double const EPS = 1e-8;
int a, b;

double get_y(double x) { return 2 * b * sqrt((1 - x * x / a / a)); }

double get_s(double l, double r) {
  return (get_y(l) + 4 * get_y((l + r) / 2) + get_y(r)) * (r - l) / 6;
}

double calc(double l, double r) {
  double mid = (l + r) / 2;
  if (fabs(get_s(l, r) - get_s(l, mid) - get_s(mid, r)) < EPS)
    return get_s(l, r);
  return calc(l, mid) + calc(mid, r);
}

int main() {
  int T;
  scanf("%d", &T);
  for (; T--;) {
    int l, r;
    scanf("%d%d%d%d", &a, &b, &l, &r);
    printf("%.3lf\n", calc(l, r));
  }
  return 0;
}