ACM-ICPC World Finals 2016 I解法

I: Road Times

問題

n 頂点 m 辺の有向グラフが与えられる。
各辺は距離[km]が与えられる。
あなたが頂点uからvへ移動するときには、必ず距離が最短に成るように移動する。
各辺は速度が定数で制限されており、それは時速30km以上60km以下の実数の定数で、あなたはその辺を通るときはちょうどその速度で移動する(要はx[km]の辺はx~2x[min.]の実数の定数時間が掛かる)
その速度の定数は分からないが、情報としてr個のヒント(出発点、到達点、時間[min.])が与えられる。
q個の質問:「頂点sからdに移動するときに掛かる、最短時間と最長時間を出力」に答えよ。

制約

グラフは隣接行列で与えられる
1 <= n <= 30
m <= 100
1 <= r <= 100
1 <= q <= 100
r + q 個の出発点・到達点の全てのペアには唯一の最短経路が存在し、r個の時間のヒントは正当で矛盾がない。

解法

線形計画法が解けますか?
問題文を意訳するとこれに尽きる。

1 (最短経路):

最短経路も自明ではないので一応考えると、サイズの小さい行列で与えられるのでFloyd-Warshallが良い。
最短経路の辺が列挙できるように出発点から向かう頂点もDPで求める。

FW[i][j] := iからjへ移動するときの最短距離
nxt[i][j] := iからjへ移動するとき、iの次に進むべき頂点

FWを隣接行列で初期化。辺が無いところは無限、FW[i][i]は使わないので0でも何でもいい。
nxt[i][j] = jで初期化

REP (k, n)
    REP (i, n)
        REP (j, n) 
            if FW[i][j] > FW[i][k] + FW[k][j] then
                FW[i][j] = FW[i][k] + FW[k][j]
                nxt[i][j] = nxt[i][k]

FW法は中継点kが外側のループである必要がある。

2 (線形計画):

長さc_e[km]の辺eを移動するときに掛かる時間x_e[min]は c_e <= x_e <= 2 c_e 。
頂点s_j, d_j間の最短経路の辺集合をE_jとして、その時の時間がt_jがわかっているので sum_{e in E_j} x_e = t_j 。
これらが線形計画問題の制約条件。
目的関数は同様に、クエリの最短経路の辺集合をEとした時、
min and max : sum_{e in E} x_e

線形計画法は単体法一択だが、コードが書けるならAC、書けないなら終わり。
私は書けないのでindy256先生のコードをc++, doubleに移植した。
Simplex algorithm - Algorithms and Data Structures

制約条件が変わらないので、単体法の計算を一部使い回すこともできるがTLEの心配はなさそう。

コード

#include<cassert>
#include<cstdio>
#include<iostream>
#include<vector>
#include<algorithm>
#include<string>
#include<cstring>
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; }


// max c * x s.t. A*x <= b, x >= 0
struct Simplex {
    const double EPS = 1e-8;
    typedef double Double;
    typedef vector<Double> Arr;
    typedef vector<Arr> Mat;
    bool LT(Double a, Double y) {
	return y-a > EPS;
    }
    bool EQ(Double a, Double y) {
	return abs(a-y) < EPS;
    }
    bool LE(Double a, Double y) {
	return y-a > -EPS;
    }

    bool infinity;
    bool none;
    Double ans;
    vector<Double> x;

    Simplex(const Mat &A, const Arr &b, const Arr &c) : infinity(false), none(false) {
	int m = A.size(), n = A[0].size() + 1;
	VI index(n+m);
	REP (i, n+m) index[i] = i;
	Mat a(m+2, Arr(n+1, 0));
	int L = m;
	REP (i, m) {
	    REP (j, n-1) a[i][j] = -A[i][j];
	    a[i][n-1] = 1;
	    a[i][n] = b[i];
	    if (a[L][n] > a[i][n]) L = i;
	}
	REP (j, n-1) a[m][j] = c[j];
	a[m+1][n-1] = -1;
	for (int E=n-1;;) {
	    if (L < m) {
		swap(index[E], index[L+n]);
		a[L][E] = 1 / a[L][E];
		REP (j, n+1) if (j != E) a[L][j] *= -a[L][E];
		REP (i, m+2) if (i != L) {
		    REP (j, n+1) if (j != E) {
			a[i][j] += a[i][E] * a[L][j];
		    }
		    a[i][E] = a[i][E] * a[L][E];
		}
	    }
	    E = -1;
	    REP (j, n) if (E < 0 || index[E] > index[j])
		if (LT(0, a[m+1][j]) || (EQ(a[m+1][j], 0) && LT(0, a[m][j])))
		    E = j;
	    if (E < 0) break;
	    L = -1;
	    REP (i, m) if (LT(a[i][E], 0))
		if (L < 0 || LE(a[L][n] / a[L][E], a[i][n] / a[i][E]))
		    L = i;
	    if (L < 0) {
		infinity = true;
		return;
	    }
	}
	if (LT(a[m+1][n], 0)) {
	    none = true;
	    return;
	}
	x.assign(n-1, 0);
	REP (i, m) if (index[n+i] < n-1) x[index[n+i]] = a[i][n];
	ans = a[m][n];
    }
};

int N, M, R, Q, Z;
int FW[31][31], nxt[31][31], E[31][31], e_id[31][31];
int H[31][31];


int main() {

    // input
    memset(e_id, -1, sizeof e_id);
    scanf("%d", &N);
    REP (i, N) REP (j, N) {
	scanf("%d", E[i]+j);

	if (i == j) FW[i][j] = 0;
	else if (E[i][j] == -1) FW[i][j] = 1<<29;
	else {
	    FW[i][j] = E[i][j];
	    e_id[i][j] = M++;
	}
	nxt[i][j] = j;
    }

    REP (k, N) REP (i, N) REP (j, N) if (FW[i][k] + FW[k][j] < FW[i][j]) {
	FW[i][j] = FW[i][k] + FW[k][j];
	nxt[i][j] = nxt[i][k];
    }

    // hints
    Simplex::Mat A;
    Simplex::Arr b;

    scanf("%d", &R);
    Z = M;
    REP ($, R) {
	int s, d, t;
	scanf("%d%d%d", &s, &d, &t);
	H[s][d] = t;

	Simplex::Arr a(Z, 0);
	double k = t;

	while (s != d) {
	    int v = nxt[s][d];
	    int id = e_id[s][v];
	    a[id] = 1;
	    k -= E[s][v];
	    s = v;
	}
	A.push_back(a);
	b.push_back(k);
	REP (i, Z) a[i] = -a[i];
	k = -k;
	A.push_back(a);
	b.push_back(k);
    }
    REP (i, N) REP (j, N) if (e_id[i][j] != -1) {
	A.push_back(Simplex::Arr(Z, 0));
	b.push_back(E[i][j]);
	A.back()[e_id[i][j]] = 1;
    }

   // queries
    scanf("%d", &Q);
    REP ($, Q) {
	int s_, d;
	scanf("%d%d", &s_, &d);
	int s = s_;

	Simplex::Arr c(Z);
	double offset = 0;
	while (s != d) {
	    int v = nxt[s][d];
	    int id = e_id[s][v];
	    c[id] = -1;
	    offset += E[s][v];
	    s = v;
	}

	double lo, hi;

	if (H[s_][d] != 0) {
	    lo = hi = H[s_][d];
	} else {
	    REP (t, 2) {
		Simplex X(A, b, c);

		if (t == 0) {
		    if (X.none) lo = H[s_][d];
		    else lo = -X.ans + offset;
		    REP (i, Z) c[i] = -c[i];
		} else {
		    if (X.none) hi = H[s_][d];
		    else hi = X.ans + offset;
		}
	    }
	}

	printf("%d %d %.9f %.9f\n", s_, d, lo, hi);
    }

    return 0;
}

日記

単体法は一般の入力に対して多項式時間で抑えられていない。
一方で指数時間掛かる入力を作るのも困難である、ましてやコンテスタントはどんなコードを出してくるか事前に知る由もない。
など考えたが、実際早いので良し。
また、

  1. この問題設定のインスタンスにおける単体法の計算量
  2. この問題の線形計画法以外の多項式時間の解法

については全く別の話である。