ICPC 2016 国内予選 G -- ワープ航法

問題文:
http://icpcsec.storage.googleapis.com/icpc2016-domestic/problems/all_ja.html#section_G

テストデータ:
http://icpc.iisf.or.jp/past-icpc/domestic2016/judgedata/G/

「二点を通る直線で切るパターンを全部試す」はよくありそうで、「円で切り分けられた領域ごとに考える（ために二円の交点を列挙する）」のもよくありそうだけど、このふたつの組み合わせは見たことがなかった気がする。

求めるのが最小値なので、領域を切り分けるといいつつ、各領域における包含のパターンだけが列挙できればよくて、計算するときはワープ位置は平面上の好きな点にとれると考えてもいいのも面白かった。

テストデータが通ることは確認している。

#include <cstdio>
#include <iostream>
#include <set>
#include <cmath>
#include <vector>

using namespace std;

typedef long long Int;

template<typename T> void chmax(T& t, T f) { if (t < f) t = f; }
template<typename T> void chmin(T& t, T f) { if (t > f) t = f; }
int in() { int x; scanf("%d", &x); return x; }
double fin() { double x; scanf("%lf", &x); return x; }

const double EPS = 1e-6;

struct Pt {
  double x, y;
  Pt(double x = 0, double y = 0) : x(x), y(y) { }
  Pt operator + (const Pt& o) const { return Pt(x + o.x, y + o.y); }
  Pt operator - (const Pt& o) const { return Pt(x - o.x, y - o.y); }
  Pt operator * (const double k) const { return Pt(k*x, k*y); }
  Pt operator / (const double k) const { return Pt(x/k, y/k); }
  double det(const Pt& o) const { return x*o.y - o.x*y; }
  double abs2() const { return x*x + y*y; }
  double abs() const { return sqrt(abs2()); }
  Pt rot(const double t) const { return Pt(x*cos(t) - y*sin(t), x*sin(t) + y*cos(t)); }
  Pt rot90() const { return Pt(-y, x); }
};

pair<Pt, Pt> pCC(Pt a, double r, Pt b, double s)
{
  double d = (b - a).abs();
  double x = (d * d + r * r - s * s) / (d * 2);
  Pt e = (b - a) / d, w = e.rot90() * sqrt(max(r * r - x * x, 0.0));
  return make_pair(a + e * x - w, a + e * x + w);
}

int N, M, A[64], B[64];
double V[64];
Pt P[32];

void acc(double& a, double& b, double& c, double t, double v) {
  a += 1.0 / v / v;
  b += -2.0 * t / v / v;
  c += t * t / v / v;
}

double minimum(double a, double b, double c) {
  if (fabs(a) < EPS) {
    return c;
  } else {
    double x = -b / (2 * a);
    return a * x * x + b * x + c;
  }
}

double calc(const set<Int>& pats, const int ps) {
  double res = 1e12;
  for (const Int pat : pats) {
    double Ax = 0.0, Bx = 0.0, Cx = 0.0;
    double Ay = 0.0, By = 0.0, Cy = 0.0;
    for (int k = 0; k < M; ++k) {
      if (!(ps & (1 << A[k]))) {
        continue;
      }
      if (pat & (1LL << k)) {
        acc(Ax, Bx, Cx, P[A[k]].x, V[k]);
        acc(Ay, By, Cy, P[A[k]].y, V[k]);
      } else {
        Cx += pow(P[A[k]].x - P[B[k]].x, 2) / V[k] / V[k];
        Cy += pow(P[A[k]].y - P[B[k]].y, 2) / V[k] / V[k];
      }
    }
    chmin(res, minimum(Ax, Bx, Cx) + minimum(Ay, By, Cy));
  }
  return res;
}

void solve() {
  set<Int> pats;
  vector<Pt> cands;
  cands.emplace_back(5000, 5000);
  for (int i = 0; i < M; ++i) {
    cands.push_back(P[A[i]] + (P[B[i]] - P[A[i]]).rot(1.0) * (1.0 - EPS));
  }
  for (int i = 0; i < M; ++i) {
    for (int j = i + 1; j < M; ++j) {
      auto tcc = pCC(P[A[i]], (P[A[i]] - P[B[i]]).abs(), P[A[j]], (P[A[j]] - P[B[j]]).abs());
      Pt cc[2] = {tcc.first, tcc.second};
      for (int k = 0; k < 2; ++k) {
        Pt dir = cc[1 - k] - cc[k];
        dir = dir / dir.abs();
        for (int d = 0; d < 4; ++d) {
          cands.push_back(cc[k] + dir * EPS);
          dir = dir.rot90();
        }
      }
    }
  }
  for (const Pt& p : cands) {
    Int pat = 0;
    for (int i = 0; i < M; ++i) {
      if ((p - P[A[i]]).abs() < (P[A[i]] - P[B[i]]).abs()) {
        pat |= 1LL << i;
      }
    }
    pats.insert(pat);
  }

  double res = 1e12;
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      if (i == j) {
        continue;
      }
      Pt v1 = (P[i] + P[j]) / 2;
      Pt v2 = v1 + (P[j] - P[i]).rot(EPS);
      int p1 = 0, p2 = 0;
      for (int k = 0; k < N; ++k) {
        if ((v2 - v1).det(P[k] - v1) < 0) {
          p1 |= 1 << k;
        } else {
          p2 |= 1 << k;
        }
      }
      chmin(res, calc(pats, p1) + calc(pats, p2));
    }
  }
  chmax(res, 0.0);

  printf("%.9f\n", sqrt(res / M));
}

int main() {
  while (true) {
    N = in();
    M = in();
    if (N == 0 && M == 0) {
      break;
    }

    for (int i = 0; i < N; ++i) {
      P[i].x = in();
      P[i].y = in();
    }
    for (int i = 0; i < M; ++i) {
      A[i] = in() - 1;
      B[i] = in() - 1;
      V[i] = fin();
    }

    solve();
  }

  return 0;
}