ACM-ICPC 2016 Asia Tsukuba Regional 国内予選, G 解答

G: ワープ航法

問題

平面上に2つのワープ場を設置したい。ワープ場にたどり着いた航空便は平面の任意の位置に瞬時に移動する事ができる。今、2次元平面上に n箇所の街がある。m本の航空便があり、航空便はある街a_iから他の街b_iへ、速さv_iで直線距離で最短のルートを選択し、移動するものである。
任意のワープポイントの置き方に対し、「各便について最短時間の2乗」の総和を最小化せよ。

制約

テストケース <= 35
n <= 20
m <= 40
座標の絶対値 <= 1000

解法

簡単な問題から段階的に難しくしていく。

Level 1: ワープ場1箇所、全ての航空便が必ずワープをする

コストの関数を式変形するとx, y独立に最短コストの座標は重心のような簡単な式で表すことができる。
航空便の集合をΛとして、ワープ場(x, y)の時のコストは

と定義されている。コストを最小にする(x, y)は

Level 2: ワープ場1箇所、ワープを使っても使わなくてもよい

ワープを使う航空便とは、目的地よりもワープ場のほうが近い場合である。
全ての航空便について、出発の街から目的地までの距離を半径として円を描くと、平面がO(m^2)に分割される。
分割された各領域の中では、どこに置いてもワープを使う航空便の集合は固定であるので「ある領域のワープを使うべき航空便の集合」はその領域の代表点が1点あればよい。これは、円の内点や、2つの円の交点を適当に出せば良い。

ワープを使うべき集合についてはLevel 1でコストを求め、それ以外の航空便は直線距離でコストを出す。各領域で求めたコストの内、最小のものがLevel 2の答え。
航空便の中心から距離を半径として円を書くと、青いドットの領域では赤い円の航空便のみワープ場を使う。

Level 3: ワープ場2箇所、ワープを使っても使わなっくても良い

本問である。
m個の航空便を2つ二分割すると、それぞれを、Level 2によって求めることができるが、分け方がO(2^m)あり厳しい。
2つのワープ場があるとすると、「たとえ、ワープをするとしたら出発する街から近い方を使う」はずである。よって、街を二分割すればよいが、これでO(2^n)に落ちたがまだ厳しい。

2つのワープ場のどちらに近いかは、ワープ場を結ぶ線分の垂直二等分線で分断される平面で決まる。つまり、任意の直線で平面を分割し、それぞれに1つづつ最適なワープ場をLevel 2で求める。
分断する直線の位置や角度は関係なく、2つに別れる街の集合だけ必要となるが、これは、「ある2つの街をギリギリ分断する直線」を全て列挙することで、全ての分割を求めることができる。
分割はO(n^2)通り。
Level 2と合わせて全体で、O(n^2 m^3)。これはジャッジ解通りの計算量だが各入力データに対し30秒程度の実行時間が掛かる。国内予選のルールなら気にする必要はないが、流石に想定解法で10秒以上掛かるのは実装が悪いと思われる。
ギリギリで平面を分断する直線の雰囲気の画像。

コード

#include<cstdio>
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<cstring>
#include<complex>
#include<cmath>
using namespace std;

typedef long long LL;
typedef vector<int> VI;

#define REP(i,n) for(int i=0, i##_len=(n); i<i##_len; ++i)
#define EACH(i,c) for(__typeof((c).begin()) i=(c).begin(),i##_end=(c).end();i!=i##_end;++i)
#define eprintf(s...)  fprintf(stderr, s)

template<class T> inline void amin(T &x, const T &y) { if (y<x) x=y; }
template<class T> inline void amax(T &x, const T &y) { if (x<y) x=y; }

typedef double Double;
typedef complex<Double> Point;
const Double PI = acos(-1);
const Double EPS = 1e-7;

const Point I(0, 1);

// abs(P a) ::= sqrt(a.x^2 + a.y^2);
// norm(P a) ::= a.x^2 + a.y^2;
// Point polar(rho, theta=0)

int sign(Double x) {
    if (x < -EPS) return -1;
    if (EPS < x) return 1;
    return 0;
}

namespace std {
    bool operator<(const Point &x, const Point &y) {
	return real(x) != real(y)? real(x) < real(y): imag(x) < imag(y);
    }
}

Double cross(const Point &x, const Point &y) {
    return imag(conj(x)*y);
}
Double dot(const Point &x, const Point &y) {
    return real(conj(x)*y);
}
Point normal(const Point &a) { return a / abs(a); }

int ccw(Point a, Point b, Point c) {
    b-=a; c-=a;
    if (cross(b, c) > EPS) return 1; // counter clockwise
    if (cross(b, c) < -EPS) return -1; // clockwise
    if (dot(b, c) < -EPS) return 2; // c--a--b on line
    if (norm(b) < norm(c)) return -2; // a--b--c on line
    return 0;
}

struct Circle {
    Point p;
    Double r;
    Circle(Point p_, Double r_) : p(p_), r(r_) {
    }
};

vector<Point> intersection(Circle c0, Circle c1) {
    vector<Point> ret;
    Double l = abs(c1.p - c0.p);
    if (l < EPS) {
	ret.push_back(c0.p);
	return ret;
    }
    if (l > c0.r + c1.r) return ret;
    if (c0.r > l + c1.r) return ret; // c1 in c0
    if (c1.r > l + c0.r) return ret; // c0 in c1


    Double A = acos((pow(c0.r, 2) + pow(l, 2) - pow(c1.r, 2)) / (2 * c0.r * l));
    Double aa = arg(c1.p - c0.p);
    for (int d=-1; d<2; d+=2) {
	for (int e=-1; e<2; e+=2) {
	    ret.push_back(c0.p + polar(c0.r + d * EPS, aa + A + e * EPS));
	    ret.push_back(c0.p + polar(c0.r + d * EPS, aa - A + e * EPS));
	}
    }

    return ret;
}


int N, M;
Point Z[22];
int A[44], B[44];
Double V[44];

Point P[2];
Double ans;
Point ansP[2];

Double half_cost(int S) {
    vector<Point> G;
    REP (i, M) if (S & (1<<A[i])) {
	Double r = abs(Z[A[i]] - Z[B[i]]);
	G.push_back(Z[A[i]]);

	REP (j, i) if (S & (1<<A[j])) {
	    Double r_j = abs(Z[A[j]] - Z[B[j]]);
	    vector<Point> points = intersection(Circle(Z[A[i]], r), Circle(Z[A[j]], r_j));
	    G.insert(G.end(), points.begin(), points.end());
	}
    }

    vector<LL> mask(G.size());

    Double ans = 1e100;

    REP (k, G.size()) {
	Double sum_v2 = 0;
	Double sum_x1 = 0;
	Double sum_y1 = 0;
	REP (i, M) if (S & (1<<A[i])) {
	    if (abs(Z[A[i]] - Z[B[i]]) > abs(Z[A[i]] - G[k])) {
		mask[k] |= 1LL<<i;
		Double v2 = 1.0 / pow(V[i], 2);
		sum_v2 += v2;
		sum_x1 += Z[A[i]].real() * v2;
		sum_y1 += Z[A[i]].imag() * v2;
	    }
	}

	if (mask[k] == 0) {
	    G[k] = Point(0, 0);
	} else {
	    G[k] = Point(sum_x1 / sum_v2, sum_y1 / sum_v2);
	}

	Double guess = 0;

	REP (i, M) if (S & (1<<A[i])) {
	    Double l = min(abs(Z[A[i]] - Z[B[i]]), abs(Z[A[i]] - G[k]));
	    guess += pow(l / V[i], 2);
	}

	if (ans > guess) {
	    ans = guess;
	}
    }

    return ans;
}

int main() {
    while (scanf("%d%d", &N, &M), N|M) {
	REP (i, N) {
	    double x, y;
	    scanf("%lf%lf", &x, &y);
	    Z[i] = Point(x, y);
	}
	REP (i, M) {
	    double v;
	    scanf("%d%d%lf", A+i, B+i, &v);
	    A[i]--;
	    B[i]--;
	    V[i] = v;
	}

	Double ans = half_cost((1<<N) - 1);


	REP (i, N) REP (j, i) {
	    Point vec = Z[j] - Z[i];
	    vec = vec / abs(vec) * polar<Double>(1.0, PI / 2.0);
	    REP (t, 2) {
		Point zi = Z[i] + vec * EPS;
		Point zj = Z[j] - vec * EPS;

		int S1 = 0, S2 = 0;
		REP (k, N) {
		    if (ccw(zi, zj, Z[k]) == 1) S1 |= 1 << k;
		    else S2 |= 1 << k;
		}

		Double cst1 = half_cost(S1);
		Double cst2 = half_cost(S2);
		if (ans > cst1 + cst2) {
		    ans = cst1 + cst2;
		}

		vec = -vec;
	    }
	}

	ans = sqrt(ans / M);
	printf("%.9lf\n", (double)ans);
    }

    return 0;
}

日記

難しすぎる。本番中通したチームはいたのだろうか。
平面を分割するのにO(n^2)で列挙するものは過去のICPCでも3次元O(n^3)のものが出題されていたはず。(地区予選だったか？)