P4682 [ZJOI2007] 粒子运动 题解

观察到数据范围很小，枚举产生最短距离的两个粒子，如果知道它们的轨迹，轨迹只有 \(k\) 段，那么双指针这两个粒子的轨迹就可以计算答案。

考虑计算一个粒子的轨迹，首先可以解出来它会在点 \(P\) 碰到圆，然后需要计算的是改变后的速度，考虑速度向量和 \(\overrightarrow {PO}\) 的夹角，反射之后的速度一定是逆时针旋转二倍夹角之后的速度。

熟知 \(a\times b = |a||b|\sin {\alpha}, a·b=|a||b|\cos \alpha\)，可以根据二倍角公式计算出 \(\sin 2\alpha, \cos 2\alpha\)，根据这两个值可以用向量旋转公式计算出旋转之后的向量。

对于两个粒子的两段轨迹，想要计算出粒子在这两段轨迹运动时的最小距离，设其中一个粒子的坐标为 \((x+v_xt, y+v_yt)\) 另一个粒子同理，根据两点距离公式得到距离是一个二次函数，找到这个二次函数的顶点 \(x\) 坐标就是最小距离的 \(t\)，需要注意二次函数定义域的问题。

// Problem: P4682 [ZJOI2007] 粒子运动
// Author: Moyou
// Copyright (c) 2025 Moyou All rights reserved.

#include <algorithm>
#include <cstring>
#include <iostream>
#include <queue>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef long double ld;
typedef pair<ld, ld> PII;
const int N = 100 + 10;

ld X, Y, R;
int n, k;
struct qwq {
    ld x, y, vx, vy, l, r;
} p[N][N], t1[N], t2[N];
ld inner(ld a, ld b, ld c, ld d) {
    return a * c + b * d;
}
ld outer(ld a, ld b, ld c, ld d) {
    return -a * d + b * c;
}
ld dis(ld a, ld b) {
    return sqrtl(a * a + b * b);
}
PII rotate(ld x, ld y, ld s, ld c) {
    return {x * c - y * s, x * s + y * c};
}
ld S(ld x) {
    return x * x;
}
ld calc(ld a, ld b, ld c, ld x) {
    return a * x * x + b * x + c;
}
signed main() {
    ios::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    cin >> X >> Y >> R >> n >> k;
    for(int i = 1; i <= n; i ++) {
        cin >> p[i][0].x >> p[i][0].y >> p[i][0].vx >> p[i][0].vy;
        p[i][0].x -= X, p[i][0].y -= Y, p[i][0].l = 0;
        for(int j = 0; j < k; j ++) {
            if(j) p[i][j].l = p[i][j - 1].r;
            ld a = p[i][j].vx * p[i][j].vx + p[i][j].vy * p[i][j].vy,
               b = 2 * (p[i][j].x * p[i][j].vx + p[i][j].y * p[i][j].vy),
               c = p[i][j].x * p[i][j].x + p[i][j].y * p[i][j].y - R * R,
               t1 = (-b + sqrtl(b * b - 4 * a * c)) / a / 2,
               t2 = (-b - sqrtl(b * b - 4 * a * c)) / a / 2;
            if(t1 < t2) swap(t1, t2);
            p[i][j].r = p[i][j].l + t1;
            ld xx = p[i][j].x + t1 * p[i][j].vx, yy = p[i][j].y + t1 * p[i][j].vy;
            p[i][j + 1].x = xx, p[i][j + 1].y = yy;
            
            // 计算碰撞后速度
            ld l1 = dis(xx, yy), l2 = dis(-p[i][j].vx, -p[i][j].vy);
            ld sina = outer(-xx, -yy, -p[i][j].vx, -p[i][j].vy) / l1 / l2,
               cosa = inner(-xx, -yy, -p[i][j].vx, -p[i][j].vy) / l1 / l2;
            ld sin2a = 2 * sina * cosa, cos2a = cosa * cosa - sina * sina;
            auto tmp = rotate(-p[i][j].vx, -p[i][j].vy, sin2a, cos2a);
            p[i][j + 1].vx = tmp.x, p[i][j + 1].vy = tmp.y;
        }
    }
    ld ans = 1e18;
    for(int i = 1; i <= n; i ++) {
        for(int j = i + 1; j <= n; j ++) {
            for(int t = 0; t < k; t ++) t1[t] = p[i][t], t2[t] = p[j][t];
            for(int p1 = 0, p2 = 0; p1 < k && p2 < k; ) {
                ld lim = min(t1[p1].r, t2[p2].r) - t1[p1].l;
                ld x = t1[p1].x - t2[p2].x, y = t1[p1].y - t2[p2].y, vx = t1[p1].vx - t2[p2].vx, vy = t1[p1].vy - t2[p2].vy;
                ld a = S(vx) + S(vy),
                   b = 2 * (x * vx + y * vy),
                   c = S(x) + S(y),
                   X = -b / (a * 2);
                ans = min(ans, calc(a, b, c, max((ld)0, min(lim, X))));
                if(t1[p1].r < t2[p2].r) t2[p2].x += t2[p2].vx * lim, t2[p2].y += t2[p2].vy * lim, t2[p2].l = t1[p1 ++].r;
                else t1[p1].x += t1[p1].vx * lim, t1[p1].y += t1[p1].vy * lim, t1[p1].l = t2[p2 ++].r;
            }
        }
    }
    ans = max(ans, (ld)0);
    printf("%.3Lf\n", sqrtl(ans));

    return 0;
}