题解 P6671 【[清华集训2016] 定向越野】

· · 题解

发篇题解纪念一下我一天的 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;
}