BZOJ 1094. [ZJOI2007]粒子运动

#include <bits/stdc++.h>

#define db long double
const db eps = 1e-8;
const db pi = acos(-1.0);
const int N = 110;

inline int sign(db k) { return k < -eps ? -1 : k > eps; }
inline int cmp(db k1, db k2) { return sign(k1 - k2); }

struct P {
  db x, y;
  P() {}
  P(db x, db y): x(x), y(y) {}
  P operator + (const P &p) const { return {x + p.x, y + p.y}; }
  P operator - (const P &p) const { return {x - p.x, y - p.y}; }
  P operator * (db k) const { return {x * k, y * k}; }
  P operator / (db k) const { return {x / k, y / k}; }
  bool operator < (const P &p) const {
    int c = cmp(x, p.x);
    return c ? c == -1 : cmp(y, p.y) == -1;
  }
  bool operator == (const P &p) const {
    return !cmp(x, p.x) && !cmp(y, p.y);
  }
  db distTo(const P &p) const { return (*this - p).abs(); }
  db alpha() { return atan2(y, x); }
  void read() { double a, b; scanf("%lf%lf", &a, &b); x = a, y = b; }
  //void print() { printf("%.10f %.10f\n", x, y); }
  db abs() { return sqrt(abs2()); }
  db abs2() { return x * x + y * y; }
  P rot(const db &k) { return P(x * cos(k) - y * sin(k), x * sin(k) + y * cos(k)); }
  P rot90() { return P(-y, x); }
  P unit() { return *this / abs(); }
  P normal() { return rot90() / abs(); }
  int quad() const { return sign(y) == 1 || (sign(y) == 0 && sign(x) >= 0); }
  P getdel() {if (sign(x) == -1 || (sign(x) == 0 && sign(y) == -1)) return (*this) * (-1); else return (*this);}
  db dot(const P &p) { return x * p.x + y * p.y; }
  db det(const P &p) { return x * p.y - y * p.x; }
} o;

struct L {
  P p, v;
  void read() { 
    p.read(); v.read(); 
  }
} A[N][N];
int n, k;
db r, cost[N][N];

// 投影
P proj(const P &p1, const P &p2, const P &q) {
  P dir = p2 - p1;
  return p1 + dir * (dir.dot(q - p1) / dir.abs2());
}

// 反射
P reflect(const P &p1, const P &p2, const P &q) {
  return proj(p1, p2, q) * 2 - q;
}

db solve(int x, int y, int k1, int k2, db t1, db t2) {
  P p = (A[x][k1].p - A[x][k1].v * cost[x][k1]) - (A[y][k2].p - A[y][k2].v * cost[y][k2]), q = A[x][k1].v - A[y][k2].v;
  db a = q.abs2(), b = 2 * p.dot(q), c = p.abs2();
  if (!a) return std::min(sqrt(t1 * b + c), sqrt(t2 * b + c));
  db t = -b / (2.0 * a);
  if (t > t1 && t < t2) return sqrt(a * t * t + b * t + c);
  return std::min(a * t1 * t1 + b * t1 + c, a * t2 * t2 + b * t2 + c);
}

int main() {
  o.read();
  double rr;
  scanf("%lf", &rr);
  r = (db)rr;
  scanf("%d%d", &n, &k);
  for (int i = 1; i <= n; i++) {
    A[i][0].read();
    A[i][0].p = A[i][0].p - o;
  }
  o.x = o.y = 0;
  for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= k; j++) {
      P p = A[i][j - 1].p, q = A[i][j - 1].v;
      db a = q.abs2(), b = 2 * p.dot(q), c = p.abs2() - r * r;
      db t = (-b + sqrt(b * b - 4 * a * c)) / (a * 2);
      cost[i][j] = cost[i][j - 1] + t;
      A[i][j].p = A[i][j - 1].p + A[i][j - 1].v * t;
      P re = reflect(o, A[i][j].p, A[i][j - 1].p);
      A[i][j].v = (re - A[i][j].p) / t;
    }
  }
  db ans = 1e10;
  for (int i = 1; i <= n; i++) 
    for (int j = i + 1; j <= n; j++) {
      int k1 = 0, k2 = 0;
      while (k1 < k && k2 < k) {
        ans = std::min(ans, solve(i, j, k1, k2, std::max(cost[i][k1], cost[j][k2]), std::min(cost[i][k1 + 1], cost[j][k2 + 1])));
        if (cost[i][k1 + 1] < cost[j][k2 + 1]) k1++;
        else k2++;
      }
    }
  printf("%.3Lf\n", ans);
  return 0;
}