[SkcRound]凸包板题

题目

题目背景
$\mathsf{D}\mathsf{D}(\mathsf{X}\mathsf{Y}\mathsf{X})$ 是一根标杆，永远地标志着战力天花板。

题意概要
给出若干个圆，随机一个角度 $\theta$，取倾斜角为 $\theta$ 的直线 $\ell$ 去 “卡壳”，问两条直线的期望距离。输出结果乘 $\pi$，精度误差不超过 $10^{-8}$ 。

数据范围与提示
$n⩽5\times 10^5$，但时限 $3\mathrm{s}$ 。

思路

设圆心坐标为 { ( x i , y i ) } \{(x_i,y_i)\} {(xi​,yi​)}，半径为 { r i } \{r_i\} {ri​}，显然就是求
$${\int }_{-\pi }^{\pi }\max (x_i\cos \theta +y_i\sin \theta +r_i)\cdot \text{d}\theta$$

直接考虑每个圆作为最大值的区间，难以使用排序后求凸包算法，毕竟该式并非直线。但我的视野太狭窄了；本来就存在 分治 求凸包的方法，只是对于二维直线的处理有更简单的方式罢了。最终我只拿了 $3\mathrm{p}\mathrm{t}\mathrm{s}$，约为 $\mathsf{D}\mathsf{D}(\mathsf{X}\mathsf{Y}\mathsf{X})$ 得分的 $6\%$ 。

类似于凸包，我们需要证明：区间数量是 $\mathcal{O}(n)$ 的。这里的 “区间” 是指，使得上式取$\max$ 的圆（即最大值点）不变的 $\theta$ 的区间。不妨限定 $\theta \in [-\pi ,\pi ]$ 。

只需要注意到一个简单的性质——这个性质是直线也具有的——两个不同的圆的 最大值点不互相跨越。也就是说，不存在单增序列 { θ 1 , θ 2 , θ 3 , θ 4 } \{\theta_1,\theta_2,\theta_3,\theta_4\} {θ1​,θ2​,θ3​,θ4​} 使得 ${\theta }_1,{\theta }_3$ 处的 $\max$ 值在圆 $a$ 处取得，而 ${\theta }_2,{\theta }_4$ 的 $\max$ 值在 $b$ 处取得。因为 $[-\pi ,\pi ]$ 中，只有 $\theta$ 恰好为两个圆的外公切线倾斜角时，才会让 $\max$ 取值点变化；变化两次，中间夹住的区间就是某个圆的 $\max$ 贡献范围，就不会互相跨越了。

只要有这个性质，就可以建立树形结构。形式化的证明：任意时刻，存在某个圆，其 $\max$ 贡献范围是一个区间（而非多个）。因为圆 $a$ 的两个区间之间，设存在圆 $b$ 的区间，则圆 $b$ 的所有区间都在此处（由上述性质）。要么 $b$ 只有一个区间，得证；要么 $b$ 在此处有至少两个区间，对这两个区间递归分析——圆只有 $n$ 个，所以递归有出口，得证。

然后归纳法。将 $\max$ 贡献范围是一个区间的圆删去，得到 $(n-1)$ 个圆的情形，由归纳法，最多 $2(n-1)$ 个区间；然后将这一个区间加入，除去自己这一个区间，至多使得区间数量 $+1$，即两侧本来是一个大区间、现在裂开了。于是最多 $2n$ 个区间，归纳得证。

外层套个分治，时间复杂度 $\mathcal{O}(n\log n)$ 。计算答案可以用开头的式子；另有一个神奇结论是，答案为 凸包的周长。不妨设这个凸包只由圆弧构成（直线部分视为曲率半径非常大的圆弧），设列向量 $\left[\begin{array}{c}x(\theta )\\ y(\theta )\end{array}\right]$ 为 $\theta$ 对应的最优解（凸包上的点），则答案为
$$\begin{array}{rl}& {\int }_{-\pi }^{\pi }\left[\begin{array}{c}x(\theta )\\ y(\theta )\end{array}\right]\cdot \left[\begin{array}{c}\cos \theta \\ \sin \theta \end{array}\right]\cdot \text{d}\theta \\ =& {\int }_{-\pi }^{\pi }x(\theta )\cdot \text{d}(\sin \theta )-y(\theta )\cdot \text{d}(\cos \theta )\\ =& [x(\theta )\sin \theta -y(\theta )\cos \theta {]}_{-\pi }^{\pi }+{\int }_{-\pi }^{\pi }\cos \theta \cdot \text{d}y(\theta )-\sin \theta \cdot \text{d}x(\theta )\\ =& {\int }_{-\pi }^{\pi }\left[\begin{array}{c}\text{d}x(\theta )\\ \text{d}y(\theta )\end{array}\right]\cdot \left[\begin{array}{c}-\sin \theta \\ \cos \theta \end{array}\right]\end{array}$$

由定义可知 $\left[\begin{array}{c}\text{d}x(\theta )\\ \text{d}y(\theta )\end{array}\right]$ 是垂直于 $\theta$ 方向的。所以最后这个式子，实际上是向量长度；这就等价于 $\oint |\overrightarrow{l}|=C$ 即周长。多么神奇！

代码细节

稍微聊聊，分治之后合并时，怎么求出两个圆各自的 $\max$ 贡献区间。代码实现能力强的，可以立刻退出这篇没用的博客了。

设圆心分别为 $m_1,m_2$，半径为 $r_1,r_2$ 。根据我们的分析，只需要解出两个临界点
$$\frac{\overrightarrow{m_1}\cdot \overrightarrow{n}}{|\overrightarrow{n}|}+r_1=\frac{\overrightarrow{m_2}\cdot \overrightarrow{n}}{|\overrightarrow{n}|}+r_2$$

我比较懒，直接令 $|\overrightarrow{n}|=1$，先求出 $\overrightarrow{n}$ 在 $(\overrightarrow{m_1}-\overrightarrow{m_2})$ 上的投影向量 $\overrightarrow{n_0}$，再做垂线与单位圆相交。

实现时，$\overrightarrow{n_0}=\frac{r_2-r_1}{|\overrightarrow{m_1}-\overrightarrow{m_2}{|}^2}(\overrightarrow{m_1}-\overrightarrow{m_2})$，求向量长度的平方是简单的。偏移量 $|v|=\sqrt{1-|\overrightarrow{n_0}{|}^2}$ 所以 $\frac{|v|}{|\overrightarrow{n_0}|}=\sqrt{\frac{1}{|\overrightarrow{n_0}{|}^2}-1}$，很容易得到。

吐槽一句：$\mathtt{c}\mathtt{m}\mathtt{a}\mathtt{t}\mathtt{h}$ 很脆弱。不仅 $\mathtt{s}\mathtt{q}\mathtt{r}\mathtt{t}(0)$ 会得到 $\mathrm{i}\mathrm{n}\mathrm{f}$，而且 $\mathrm{i}\mathrm{n}\mathrm{f}$ 乘 $0$ 又会得到 $\mathrm{n}\mathrm{a}\mathrm{n}$，完全搞不懂 😂

完整代码

我还是按照最初的计算式求的。

坑：当相邻区间对应的最优圆是同一个时，应当合并。我忘了这茬……但仍然过了……

#include <cstdio> // JZM yydJUNK!!!
#include <iostream> // XJX yyds!!!
#include <algorithm> // XYX yydLONELY!!!
#include <cstring> // (the STRONG LONG for loneliness)
#include <cctype> // DDG yydDOG & ZXY yydBUS !!!
#include <initializer_list>
#include <vector>
# define _USE_MATH_DEFINES
#include <cmath>
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define rep0(i,a,b) for(int i=(a); i!=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long llong;
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MAXN = 500005;
const double EPS = 1e-10;
struct Circle{ int x, y, r; };
typedef std::pair<double,double> PDD;
# define v_len_2(x,y) ((x)*(x)+(y)*(y))
/// @return the angle (radian) of inclination
PDD getTangent(Circle &a, Circle &b){
	llong vx = a.x-b.x, vy = a.y-b.y; // M1-M2
	if(!vx && !vy) return PDD{-M_PI,M_PI};
	if(a.r == b.r){ // special case
		double theta = atan2(-vx,vy); // (y,-x)
		if(theta > 0) return PDD{theta,theta-M_PI};
		else return PDD{theta,theta+M_PI};
	}
	double ratio = double(b.r-a.r)/v_len_2(vx,vy);
	double nx = ratio*vx, ny = ratio*vy; // fake
	if(v_len_2(nx,ny) > 1-EPS) // no solution
		return PDD{-M_PI,M_PI}; // cover all
	ratio = sqrt(1/v_len_2(nx,ny)-1); // reuse
	double dx = ratio*ny, dy = -ratio*nx; // (y,-x)
	PDD res = {atan2(ny+dy,nx+dx),atan2(ny-dy,nx-dx)};
	return res;
}

Circle circ[MAXN];
typedef std::pair<double,int> PDI;
typedef std::vector<PDI> vecset;
void calc(double l, double r, int a, int b, vecset &v){
	if(l >= r) return ; // empty range
	PDD qx = getTangent(circ[a],circ[b]); // common tangent
	if(qx.first > qx.second) std::swap(qx.first,qx.second);
	int core; ///< which one is inside
	if(qx.first < 0 && 0 < qx.second)
		core = (circ[a].x+circ[a].r > circ[b].x+circ[b].r) ? a : b;
	else core = (circ[a].x+circ[a].r < circ[b].x+circ[b].r) ? a : b;
	if(qx.first <= l && r <= qx.second) // in
		return void(v.push_back(PDI{r,core}));
	if(r <= qx.first || qx.second <= l) // out
		return void(v.push_back(PDI{r,a^b^core}));
	if(r <= qx.second){ // two parts
		v.push_back(PDI{qx.first,a^b^core});
		v.push_back(PDI{r,core});
	}
	else if(qx.first <= l){ // two parts
		v.push_back(PDI{qx.second,core});
		v.push_back(PDI{r,a^b^core});
	}
	else{ // three parts
		v.push_back(PDI{qx.first,a^b^core});
		v.push_back(PDI{qx.second,core});
		v.push_back(PDI{r,a^b^core});
	}
}
vecset solve(int l, int r){
	if(l == r) return vecset{PDI{M_PI,l}};
	vecset &&lc = solve(l,(l+r)>>1);
	vecset &&rc = solve((l+r+2)>>1,r);
	int lenl = int(lc.size()), lenr = int(rc.size());
	double lst = -M_PI; vecset res; res.clear();
	for(int i=0,j=0; i!=lenl&&j!=lenr; )
		if(lc[i].first < rc[j].first){
			calc(lst,lc[i].first,lc[i].second,rc[j].second,res);
			lst = lc[i].first, ++ i; // merge sort
		}
		else{
			calc(lst,rc[j].first,lc[i].second,rc[j].second,res);
			lst = rc[j].first, ++ j; // keep moving
		}
	return res;
}

double contribute(double l, double r, int id){
	double kx = sin(r)-sin(l), ky = cos(r)-cos(l);
	return kx*circ[id].x-ky*circ[id].y+(r-l)*circ[id].r;
}
int main(){
	int n = readint();
	rep(i,1,n){
		circ[i].x = readint(), circ[i].y = readint();
		circ[i].r = readint();
	}
	vecset &&convex = solve(1,n);
	double lst = -M_PI, ans = 0;
	for(int len=int(convex.size()),i=0; i!=len; ++i)
		ans += contribute(lst,convex[i].first,
			convex[i].second), lst = convex[i].first;
	printf("%.12f\n",ans);
	return 0;
}