Solution -「SPOJ-VCIRCLES」Area of Circles

(mathcal{Description})

  Link.

  求平面上 (n) 个圆的并的面积。

  (nle50)，可能被圆覆盖的横纵坐标区域在 ([-10^4,10^4])。

(mathcal{Solution})

  做了那么多计算几何之后写了这道不那么计算几何的计算几何题的题解。

  若想直接处理面积，就需要处理圆的各种相交关系，但是仅计算在某个 (x_0) 处，被圆的并覆盖的线段长度是很好处理地：每次 (mathcal O(nlog n)) 排序后扫一遍即可。我们记 (f(x)) 表示 (x_0=x) 时，被圆的并覆盖的线段长度。所以问题转为求 (int_{-infty}^{+infty}f(x)dx)。

  运用自适应 Simpson 积分，有 Simpson 公式：

[int_a^bf(x) ext dxapproxfrac{b-a}6left(f(a)+4f(frac{a+b}2)+f(b) ight) ]

  其本质是用二次函数 (g) 近似地替换一个复杂的函数 (f)，即设：

[g(x)=Ax^2+Bx+Capprox f(x)Longrightarrowint_a^bf(x) ext dxapproxint_b^ag(x) ext dx ]

  利用 (g(x)) 化简积分，注意 (g) 仅是一种形式，我们的目的使用 (f) 的一些点值近似表示原积分。推导：

[egin{aligned}int_a^bf(x) ext dx&approxint_a^b(Ax^2+Bx+C) ext dx\&=frac{A}3(b^3-a^3)+frac{B}2(b^2-a^2)+C(b-a)\&=frac{2A(b^3-a^3)+3B(b^2-a^2)+6C(b-a)}6\&=frac{2A(b-a)(a^2+ab+b^2)+3B(b-a)(b+a)+6C(b-a)}6\&=frac{b-a}6left(2Aa^2+2Aab+2Ab^2+3Bb+3Ba+6C ight)\&=frac{b-a}6left[(Aa^2+Ba+C)+(Ab^2+Bb+C)+4left(A(frac{a+b}2)^2+Bcdotfrac{a+b}{2}+C ight) ight]\&=frac{b-a}6left(g(a)+4g(frac{a+b}2)+g(b) ight)\&approxfrac{b-a}6left(f(a)+4f(frac{a+b}2)+f(b) ight)end{aligned} ]

  记 (h(a,b)=frac{b-a}6left(f(a)+4f(frac{a+b}2)+f(b) ight))，自适应 Simpson 积分的表达式如下：

[int_a^bf(x) ext dxapproxegin{cases}h(a,b)&h(a,b)approx^*h(a,frac{a+b}2)+h(frac{a+b}2,b)\int_a^{frac{a+b}2}f(x) ext dx+int_{frac{a+b}2}^bf(x) ext dx& ext{otherwise}end{cases} ]

  其中 (approx^*) 根据题目精度要求等自行判断。若精度不合适，则递归求解提升精确度。

  注意自适应 Simpson 积分只能处理平滑函数。例如依照上述算法，(int_0^{2pi}|sin x| ext dx=0)（大雾）。

  回到本题，结合求 (mathcal O(nlog n)) 求 (f(x_0)) 的方法近似求出 (int_{-infty}^{+infty}f(x)dx) 即可。复杂度与圆的位置和要求精度有关（aka. 我不知道 qwq）。

(mathcal {Code})

  讲道理若某两块并在 (x) 轴上的投影相离应该分开求积分，但那样写我 WA 掉了 qwq。

/* Clearink */

#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>

typedef std::pair<double, double> PDD;

const int MAXN = 50;
const double EPS = 1e-9;
int n;
bool ban[MAXN + 5];

inline double sqr ( const double a ) { return a * a; }
inline double dabs ( const double a ) { return a < 0 ? -a : a; }
inline double dmin ( const double a, const double b ) { return a < b ? a : b; }
inline double dmax ( const double a, const double b ) { return a < b ? b : a; }

struct Circle { double x, y, r; } O[MAXN + 5];

inline double func ( const double x ) {
	static std::vector<PDD> sec; sec.clear ();
	for ( int i = 1; i <= n; ++i ) {
		const Circle& C ( O[i] );
		if ( dabs ( x - C.x ) < C.r ) {
			double yd = sqrt ( sqr ( C.r ) - sqr ( dabs ( C.x - x ) ) );
			sec.push_back ( { C.y - yd, C.y + yd } );
		}
	}
	std::sort ( sec.begin (), sec.end () );
	double las = -1e9, ret = 0;
	for ( PDD s: sec ) {
		las = dmax ( las, s.first );
		ret += dmax ( 0, s.second - las );
		las = dmax ( las, s.second );
	}
	return ret;
}

inline double scalc ( const double lx, const double rx,
	const double fl, const double fm, const double fr ) {
	return ( fl + 4 * fm + fr ) / 6 * ( rx - lx );
}

inline double simpson ( const double lx, const double rx,
	const double fl, const double fm, const double fr ) {
	double mx = 0.5 * ( lx + rx );
	double fll = func ( 0.5 * ( lx + mx ) ), frr = func ( 0.5 * ( mx + rx ) );
	double s = scalc ( lx, rx, fl, fm, fr );
	double sl = scalc ( lx, mx, fl, fll, fm ), sr = scalc ( mx, rx, fm, frr, fr );
	if ( dabs ( s - sl - sr ) < EPS ) return sl + sr;
	return simpson ( lx, mx, fl, fll, fm ) + simpson ( mx, rx, fm, frr, fr );
}

int main () {
	scanf ( "%d", &n );
	double mn = 1e9, mx = -1e9;
	for ( int i = 1; i <= n; ++i ) {
		scanf ( "%lf %lf %lf", &O[i].x, &O[i].y, &O[i].r );
		mn = dmin ( mn, O[i].x - O[i].r );
		mx = dmax ( mx, O[i].x + O[i].r );
	}
	printf ( "%.5f
", simpson ( mn, mx,
		func ( mn ), func ( 0.5 * ( mn + mx ) ), func ( mx ) ) );
	return 0;
}