自动搬运

查看原文

来自洛谷，原作者为

	81179332_ 这个家伙很菜，什么也没有留下

搬运于2025-08-24 22:22:57，当前版本为作者最后更新于2020-08-04 20:53:57，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

发篇题解纪念一下我一天的 debug

把起点和终点看成两个半径为 $0$ 的圆，最短路只会经过圆弧和圆之间的公切线

对于两个圆一共有两组四条公切线，我们需要求出切点坐标

对于交叉的两条公切线，我们将垂直于切线的半径画出来，连接圆心，发现有两对相似三角形，可直接求出圆心连线和半径的夹角，套用向量旋转公式即可

对于不交叉的两条公切线，我们从一个圆心向另一个圆的半径作垂线，出现了一个直角三角形，它两边已知，则可求出圆心连线和半径的夹角，进而求出切点坐标

我们需要判断公切线是否经过了其他的圆，由于线段的端点不会在圆的内部，所以一条线段与圆相交当且仅当圆心到线段的距离不大于半径，并且圆心与线段端点的连线与线段的夹角为锐角

对于圆弧，我们将同一个圆上的点极角排序，直接 $l=\theta r$ 就好啦

最后我们跑一个最短路算法

一下是我犯过的部分不那么智障的错误

公切线一共有 $n^2$ 条，每次 check 是 $O(n)$ 的，对于最短路，我们有 $n^2$ 个点，$n^2$ 条边，所以总时间复杂度为 $O(n^3+\text{最短路复杂度})$

这看起来就很不友好，计算几何常数还是挺大的，虽然在这种图 SPFA 应该不会被卡，但是它还是比 Dijkstra 慢了很多，所以如果你要使用 SPFA，你可能需要一点高超的卡常技巧

注意，如果你的圆使用了一个结构体来存储，并且你将圆上的点都压入了结构体的 vector 中，切记在函数传参时传入地址即 const circle &c 而不是只是一个 circle c，否则复杂度是 $O(n^4)$ 的

将不合法的公切线端点从圆的 vector 中删去，否则会引发各类玄学问题

为了通过 UOJ 的存在半径为 $0$ 的圆的 hack 数据，请将相切不视为相交

求圆弧长度的时候，半径为 $0$ 特判掉，否则 $0/0$ 会出现 $nan$

对于半径为 $0$ 的圆，它和其它圆的公切线只计算“不相交”的两条，如果计算了四条会出现 $nan$

最后一个小伙伴犯的错误：尽量不要使用解析几何，并请规范写法，否则会炸精度

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cctype>
#include<queue>
using namespace std;
#define C const
#define eps 0.0000000001
#define fi first
#define se second
int read()
{
	int x = 0;int f = 1;char ch = getchar();
	for(;!isdigit(ch);ch = getchar()) if(ch == '-') f = -1;
	for(;isdigit(ch);ch = getchar()) x = x * 10 + (ch ^ 48);
	return x * f;
}
#define Point Vector
C int N = 510;
struct Vector
{
	double x,y;
 	inline friend Vector operator + (C Vector &a,C Vector &b) { return {a.x + b.x,a.y + b.y}; }
	inline friend Vector operator - (C Vector &a,C Vector &b) { return {a.x - b.x,a.y - b.y}; }
	inline friend double operator * (C Vector &a,C Vector &b) { return a.x * b.x + a.y * b.y; }
	inline friend double operator ^ (C Vector &a,C Vector &b) { return a.x * b.y - b.x * a.y; }
	inline double len() C{ return sqrt((*this) * (*this)); }
	inline double dis(C Point &a,C Point &b) C{ return fabs((*this - a) ^ (*this - b)) / (b - a).len(); }
}p[N * N * 4];int tot,n;
struct edge { int to,nxt;double w; }e[N * N * 12];int head[N * N * 4],cnt;
inline void add(int u,int v,double w) { e[++cnt] = {v,head[u],w},head[u] = cnt; }
inline void add_edge(int u,int v,double w) { add(u,v,w),add(v,u,w); }
inline double angle(C Vector &a,C Vector &b) { return acos(a * b / a.len() / b.len()); }
inline Vector rotate(C Vector &a,double b)
{
	double c = cos(b),s = sin(b);
	return {a.x * c - a.y * s,a.x * s + a.y * c};
}

struct circle
{
	Point c;double r;vector<int> id;
	inline void Pop(int q)
	{
		int now = id.back();id.pop_back();
		if(now == q) return;
		id.pop_back(),id.push_back(now);
	}
}c[N];
inline bool cross(C circle &x,C Point &a,C Point &b)
{
	if(x.r <= x.c.dis(a,b)) return 0;
	return (x.c - a) * (b - a) > 0 && (x.c - b) * (a - b) > 0;
}
inline void solve1(circle &a,circle &b)
{
	Vector v = b.c - a.c;
	double k = acos((a.r - b.r) / v.len());
	Vector p1 = rotate(v,k),p2 = rotate(v,-k);
	double l = p1.len();
	p1 = a.c + (Vector){p1.x / l * a.r,p1.y / l * a.r};
	p2 = a.c + (Vector){p2.x / l * a.r,p2.y / l * a.r};
	p[++tot] = p1,a.id.push_back(tot);
	p[++tot] = p2,a.id.push_back(tot);
}
inline void solve2(circle &a,circle &b)
{
	Vector v = b.c - a.c;
	double k = acos((a.r + b.r) / v.len());
	Vector p1 = rotate(v,k),p2 = rotate(v,-k);
	double l = p1.len();
	p1 = a.c + (Vector){p1.x / l * a.r,p1.y / l * a.r};
	p2 = a.c + (Vector){p2.x / l * a.r,p2.y / l * a.r};
	p[++tot] = p1,a.id.push_back(tot);
	p[++tot] = p2,a.id.push_back(tot);
}
inline bool check(C Point &a,C Point &b,int x,int y)
{
	for(int i = 1;i <= n;i++) if(i != x && i != y)
		if(cross(c[i],a,b)) return 0;
	return 1;
}
inline void add(int i,int j)
{
	solve1(c[i],c[j]),solve1(c[j],c[i]);
	if(check(p[tot],p[tot - 3],i,j))
		add_edge(tot - 3,tot,(p[tot] - p[tot - 3]).len());
	else c[i].Pop(tot - 3),c[j].Pop(tot);
	if(check(p[tot - 2],p[tot - 1],i,j))
		add_edge(tot - 2,tot - 1,(p[tot - 2] - p[tot - 1]).len());
	else c[i].Pop(tot - 2),c[j].Pop(tot - 1);
	if(fabs(c[i].r) < eps || fabs(c[j].r) < eps) return;
	solve2(c[i],c[j]),solve2(c[j],c[i]);
	if(check(p[tot - 3],p[tot - 1],i,j))
		add_edge(tot - 3,tot - 1,(p[tot - 3] - p[tot - 1]).len());
	else c[i].Pop(tot - 3),c[j].Pop(tot - 1);
	if(check(p[tot - 2],p[tot],i,j))
		add_edge(tot - 2,tot,(p[tot - 2] - p[tot]).len());
	else c[i].Pop(tot - 2),c[j].Pop(tot);
}
inline double get_len(C circle &c,C Point &a,C Point &b)
{
	if(fabs(c.r) < eps) return 0;
	return angle(a - c.c,b - c.c) * c.r;
}
double dis[N * N * 4];priority_queue<pair<double,int> > q;bool vis[N * N * 4];
double Dijkstra(int S,int T)
{
	for(int i = 1;i <= tot;i++) dis[i] = 1e100;
	q.push({0,S}),dis[S] = 0;
	while(!q.empty())
	{
		int u = q.top().se;q.pop();
		if(vis[u]) continue;vis[u] = 1;
		for(int i = head[u],v;i;i = e[i].nxt)
			if(dis[v = e[i].to] > dis[u] + e[i].w)
				dis[v] = dis[u] + e[i].w,q.push({-dis[v],v});
	}return dis[T];
}
int main()
{
	scanf("%lf %lf %lf %lf",&c[0].c.x,&c[0].c.y,&c[N - 1].c.x,&c[N - 1].c.y);
	n = read();c[n + 1] = c[N - 1];
	for(int i = 1;i <= n;i++) scanf("%lf %lf %lf",&c[i].c.x,&c[i].c.y,&c[i].r);
	for(int i = 0;i <= n;i++)
		for(int j = i + 1;j <= n + 1;j++)
			add(i,j);
	for(int i = 0;i <= n + 1;i++)
	{
		#define v c[i].id
		sort(v.begin(),v.end(),[&](int a,int b)
		{
			Vector x = p[a] - c[i].c,y = p[b] - c[i].c;
			return atan2(x.y,x.x) < atan2(y.y,y.x);
		});int s = v.size();
		for(int j = 0;j < s;j++)
		{
			int now = v[j],nxt = v[(j + 1) % s];
			add_edge(now,nxt,get_len(c[i],p[now],p[nxt]));
		}
		#undef v
	}
	printf("%.1lf",Dijkstra(c[0].id[0],c[n + 1].id[0]));
	return 0;
}