結果
| 問題 |
No.1776 Love Triangle 2 (Hard)
|
| ユーザー |
 hitonanode
|
| 提出日時 | 2021-12-05 10:53:04 |
| 言語 | C++17(clang) (17.0.6 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 3,784 ms / 10,000 ms |
| コード長 | 13,395 bytes |
| コンパイル時間 | 5,438 ms |
| コンパイル使用メモリ | 157,720 KB |
| 実行使用メモリ | 6,572 KB |
| 最終ジャッジ日時 | 2024-07-07 07:10:46 |
| 合計ジャッジ時間 | 188,759 ms |
|
ジャッジサーバーID （参考情報） |
judge1 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 176 |
ソースコード
// O(n^4 log n) weighted-linear-matroid-parity-based solution
#include <algorithm>
#include <cassert>
#include <iostream>
#include <random>
#include <utility>
#include <vector>
using namespace std;
#include <atcoder/mincostflow>
#include <atcoder/modint>
using mint = atcoder::static_modint<(1 << 30) + 3>;
mt19937 mt(530629);
uniform_int_distribution<int> rndgen(0, mint::mod());
// Upper Hessenberg reduction of square matrices
// Complexity: O(n^3)
// Reference:
// http://www.phys.uri.edu/nigh/NumRec/bookfpdf/f11-5.pdf
template <class Tp> void hessenberg_reduction(std::vector<std::vector<Tp>> &M) {
assert(M.size() == M[0].size());
const int N = M.size();
for (int r = 0; r < N - 2; r++) {
int piv = -1;
for (int j = r + 1; j < N; ++j) if (M[j][r] != 0) {
piv = j;
break;
}
if (piv < 0) continue;
for (int i = 0; i < N; i++) std::swap(M[r + 1][i], M[piv][i]);
for (int i = 0; i < N; i++) std::swap(M[i][r + 1], M[i][piv]);
const auto rinv = Tp(1) / M[r + 1][r];
for (int i = r + 2; i < N; i++) {
const auto n = M[i][r] * rinv;
for (int j = 0; j < N; j++) M[i][j] -= M[r + 1][j] * n;
for (int j = 0; j < N; j++) M[j][r + 1] += M[j][i] * n;
}
}
}
// Characteristic polynomial of matrix M (|xI - M|)
// Complexity: O(n^3)
// R. Rehman, I. C. Ipsen, "La Budde's Method for Computing Characteristic Polynomials," 2011.
template <class Tp> std::vector<Tp> characteristic_poly(std::vector<std::vector<Tp>> &M) {
hessenberg_reduction(M);
const int N = M.size();
std::vector<std::vector<Tp>> p(N + 1); // p[i + 1] = (Characteristic polynomial of i-th leading principal minor)
p[0] = {1};
for (int i = 0; i < N; i++) {
p[i + 1].assign(i + 2, 0);
for (int j = 0; j < i + 1; j++) p[i + 1][j + 1] += p[i][j];
for (int j = 0; j < i + 1; j++) p[i + 1][j] -= p[i][j] * M[i][i];
Tp betas = 1;
for (int j = i - 1; j >= 0; j--) {
betas *= M[j + 1][j];
Tp hb = -M[j][i] * betas;
for (int k = 0; k < j + 1; k++) p[i + 1][k] += hb * p[j][k];
}
}
return p[N];
}
// (X, N), (Y, N + 1), (Z, N + 2) shortest S-paths, without using vertices in `used_vs`
// return -1 if not found
int shortest_len(int N, int X, int Y, int Z, const vector<vector<int>> &conn, const vector<int> &used_vs) {
vector<int> is_alive(N + 5, 1);
for (auto i : used_vs) is_alive[i] = 0;
vector<int> alive_vs;
for (int i = 0; i < int(is_alive.size()); ++i) {
if (is_alive[i]) alive_vs.push_back(i);
}
vector<int> vid(N + 5, -1);
for (int j = 0; j < int(alive_vs.size()); ++j) vid[alive_vs[j]] = j;
const vector<int> terminal_vs{X, Y, Z, N, N + 1, N + 2, N + 3, N + 4};
auto Label = [&](int i) -> int {
if (i == X or i == N) return 1;
if (i == Y or i == N + 1) return 2;
if (i == Z or i == N + 2) return 3;
if (i == N + 3) return 4;
if (i == N + 4) return 5;
return 0;
};
// 論文通りに実装すると Z は 2(N + 5) * M 行列となるが，all zero な行が 8 つ発生して厄介
// なので，これらを潰して (2N + 2) * M 行列（feasible ならばこれは行フルランク）を構築．
// 潰した行を飛ばした index を取得する関数
auto truncate = [&](int i) -> int {
int red = 0;
for (auto l : terminal_vs) {
if (l * 2 < i) ++red;
}
return vid[i / 2] * 2 + (i % 2) - red;
};
int r = alive_vs.size() * 2 - 8;
// M(x) = mat0 + mat1 * x
// mat = [mat1, mat0] を考える
vector mat(r, vector<mint>(r * 2));
auto add_edge = [&](int u, int v, int w) {
vector<pair<int, mint>> bret, cret;
if (Label(u)) {
bret.emplace_back(truncate(u * 2 + 1), -Label(u));
} else {
bret.emplace_back(truncate(u * 2), 1);
}
if (Label(v)) {
bret.emplace_back(truncate(v * 2 + 1), Label(v));
} else {
bret.emplace_back(truncate(v * 2), -1);
}
cret.emplace_back(truncate(u * 2 + 1), 1);
cret.emplace_back(truncate(v * 2 + 1), -1);
mint x = rndgen(mt);
for (auto [i, bi] : bret) {
for (auto [j, cj] : cret) {
auto v = bi * cj * x;
if (w == 0) {
mat[i][r + j] += v;
mat[j][r + i] -= v;
mat[i][j] += v;
mat[j][i] -= v;
}
if (w == 1) {
mat[i][j] += v;
mat[j][i] -= v;
}
}
}
};
for (int i = 0; i < N + 3; ++i) {
if (!is_alive[i]) continue;
for (int j = 0; j < N + 3; ++j) {
if (!conn[i][j] or i > j or !is_alive[j]) continue;
add_edge(i, j, 1);
}
if (!Label(i)) add_edge(i, N + 3, 0);
}
add_edge(N + 3, N + 4, 0);
// det(x mat1 + mat0) を x の多項式として求めたい
// mat1 を掃き出して det(x \lambda I - A) の形に帰着させる
for (int i = 0; i < r; ++i) {
int piv = -1;
for (int h = i; h < r; ++h) {
if (mat[h][i] != 0) piv = h;
}
if (piv < 0) return -1;
assert(piv >= i);
swap(mat[i], mat[piv]);
mint inv = mat[i][i].inv();
for (auto &x : mat[i]) x *= inv;
for (int h = 0; h < r; ++h) {
if (h == i) continue;
if (mat[h][i] == 0) continue;
const mint coeff = mat[h][i];
for (int w = 0; w < r * 2; ++w) {
mat[h][w] -= coeff * mat[i][w];
}
}
}
// det(xI - A) : https://judge.yosupo.jp/problem/characteristic_polynomial
for (auto &v : mat) {
v.erase(v.begin(), v.begin() + r);
for (auto &x : v) x = -x;
}
auto det_poly = characteristic_poly<mint>(mat);
int ret = 0;
while (ret < int(det_poly.size()) and det_poly[ret] == 0) ++ret;
if (ret < int(det_poly.size())) {
return ret / 2;
} else {
return -1;
}
}
// used_vs に含まれる頂点は使わずに，from -> to1 と from->to2 の点素なパスを構成する．
// 両方のパスが構築できなければ empty vector の組を返す．
pair<vector<int>, vector<int>> twopaths(int N, const vector<vector<int>> &to, const vector<int> &used_vs, int from, int to1, int to2) {
const int gt = N * 2;
atcoder::mcf_graph<int, int> graph(gt + 1);
vector<int> valid_v(N, 1);
for (auto i : used_vs) valid_v[i] = 0;
valid_v[to1] = valid_v[to2] = 1;
for (int i = 0; i < N; ++i) { graph.add_edge(i, i + N, valid_v[i], 0); }
for (int i = 0; i < N; ++i) {
for (auto j : to[i]) {
graph.add_edge(i + N, j, 1, 1);
}
}
graph.add_edge(to1 + N, gt, 1, 0);
graph.add_edge(to2 + N, gt, 1, 0);
auto f = graph.flow(from + N, gt, 2);
if (f.first < 2) return {{}, {}};
vector<int> conn(N);
for (auto e : graph.edges()) {
if (e.flow) {
if (e.to == gt) continue;
int s = e.from % N, t = e.to % N;
conn[s] ^= t;
conn[t] ^= s; // ループがないので生えている辺の xor だけ持っておけば後で解が復元できる
}
}
vector<int> ret1, ret2;
while (to1 != from) {
ret1.push_back(to1);
to1 = conn[to1];
conn[to1] ^= ret1.back();
}
while (to2 != from) {
ret2.push_back(to2);
to2 = conn[to2];
conn[to2] ^= ret2.back();
}
ret1.push_back(from);
ret2.push_back(from);
reverse(ret1.begin(), ret1.end());
reverse(ret2.begin(), ret2.end());
return {ret1, ret2};
}
template <typename T, T INF = std::numeric_limits<T>::max() / 2, int INVALID = -1> struct ShortestPath {
int V, E;
bool single_positive_weight;
T wmin, wmax;
std::vector<std::vector<std::pair<int, T>>> to;
ShortestPath(int V = 0) : V(V), E(0), single_positive_weight(true), wmin(0), wmax(0), to(V) {}
void add_edge(int s, int t, T w) {
assert(0 <= s and s < V);
assert(0 <= t and t < V);
to[s].emplace_back(t, w);
E++;
if (w > 0 and wmax > 0 and wmax != w) single_positive_weight = false;
wmin = std::min(wmin, w);
wmax = std::max(wmax, w);
}
std::vector<T> dist;
std::vector<int> prev;
void ZeroOneBFS(int s) {
assert(0 <= s and s < V);
dist.assign(V, INF), prev.assign(V, INVALID);
dist[s] = 0;
std::deque<int> que;
que.push_back(s);
while (!que.empty()) {
int v = que.front();
que.pop_front();
for (auto nx : to[v]) {
T dnx = dist[v] + nx.second;
if (dist[nx.first] > dnx) {
dist[nx.first] = dnx, prev[nx.first] = v;
if (nx.second) {
que.push_back(nx.first);
} else {
que.push_front(nx.first);
}
}
}
}
}
};
// a から b を経由し c に辿り着く点素なパスで，banned にあるものを使わないもののうち
// 最短で辞書順最小のものを求め，{[a, ..., b], [b, ..., c]} という pair of vector を返す．
pair<vector<int>, vector<int>> refine_path(int N, const vector<vector<int>> &to, int a, int b, int c, vector<int> banned) {
vector<int> is_banned(N);
for (auto i : banned) is_banned[i] = 1;
auto [v1, v2] = twopaths(N, to, banned, b, a, c);
const int len = v1.size() + v2.size();
vector<int> ab{a};
banned.push_back(a);
is_banned[a] = 1;
while (ab.back() != b) {
int cur = ab.back();
for (auto j : to[cur]) {
if (j == b) {
ab.push_back(b);
break;
}
if (is_banned[j]) continue;
if (j == c) continue;
auto [p1, p2] = twopaths(N, to, banned, b, j, c);
if (p1.empty()) continue;
if (int(p1.size() + p2.size() + ab.size()) != len) continue;
ab.push_back(j);
banned.push_back(j);
is_banned[j] = 1;
break;
}
}
ShortestPath<int, 1 << 20> sssp(N);
for (int i = 0; i < N; ++i) {
if (is_banned[i]) continue;
for (auto j : to[i]) {
if (is_banned[j]) continue;
sssp.add_edge(i, j, 1);
}
}
sssp.ZeroOneBFS(b);
const auto Db = sssp.dist;
sssp.ZeroOneBFS(c);
const auto Dc = sssp.dist;
vector<int> bc{b};
int cur = b;
while (cur != c) {
for (auto nxt : to[cur]) {
if (Db[nxt] == Db[cur] + 1 and Dc[nxt] == Dc[cur] - 1) {
cur = nxt;
bc.push_back(cur);
break;
}
}
}
return {ab, bc};
}
int main() {
cin.tie(nullptr), ios::sync_with_stdio(false);
int N, M, X, Y, Z;
cin >> N >> M >> X >> Y >> Z;
--X, --Y, --Z;
vector conn(N + 3, vector<int>(N + 3, 1));
for (int i = 0; i < N + 3; ++i) conn[i][i] = 0;
while (M--) {
int u, v;
cin >> u >> v;
--u, --v;
vector<int> us{u}, vs{v};
if (u == X) us.push_back(N);
if (u == Y) us.push_back(N + 1);
if (u == Z) us.push_back(N + 2);
if (v == X) vs.push_back(N);
if (v == Y) vs.push_back(N + 1);
if (v == Z) vs.push_back(N + 2);
for (auto i : us) {
for (auto j : vs) conn[i][j] = conn[j][i] = 0;
}
}
int ans_len = shortest_len(N, X, Y, Z, conn, {});
cout << ans_len << '\n';
if (ans_len < 0) return 0;
vector<int> ret;
vector<int> is_used(N + 3);
int cur = X;
while (true) {
vector<int> next_steps;
for (int i = 0; i < N + 3; ++i) {
if (conn[cur][i] and !is_used[i]) next_steps.push_back(i);
}
int ng = 0, ok = next_steps.size();
while (ok - ng > 1) {
int c = (ok + ng) / 2;
for (int i = c; i < int(next_steps.size()); ++i) {
conn[cur][next_steps[i]] = conn[next_steps[i]][cur] = 0;
}
auto d = shortest_len(N, cur, Y, Z, conn, ret);
if (d < 0 or d + int(ret.size()) > ans_len) {
ng = c;
} else {
ok = c;
}
for (int i = c; i < int(next_steps.size()); ++i) {
conn[cur][next_steps[i]] = conn[next_steps[i]][cur] = 1;
}
}
int nxt = next_steps.at(ng);
ret.push_back(cur);
is_used[cur] = 1;
cur = nxt;
if (cur == Y or cur == Z) break;
}
if (cur != Z) swap(Y, Z);
// Z -> Y -> X 辞書順最小
vector<vector<int>> to(N);
for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) {
if (conn[i][j]) to[i].push_back(j);
}
}
auto [ZY, YX] = refine_path(N, to, Z, Y, X, vector<int>(ret.begin() + 1, ret.end()));
for (auto v : ZY) ret.push_back(v);
ret.pop_back();
for (auto v : YX) ret.push_back(v);
for (auto x : ret) cout << x + 1 << ' ';
cout << '\n';
}