#721943 (C++17(clang)) No.1776 Love Triangle 2 (Hard)

• 問題: No.1776 Love Triangle 2 (Hard)
• ユーザー: hitonanode
• 提出日時: 2021-12-05 10:53:04
• 言語: C++17(clang) (17.0.6 + boost 1.83.0)
• 結果: AC
• 実行時間: 3,669 ms / 10,000 ms
• コード長: 13,395 bytes
• コンパイル時間: 7,573 ms
• コンパイル使用メモリ: 120,848 KB
• 実行使用メモリ: 6,528 KB
• 最終ジャッジ日時: 2023-09-21 13:28:03
• 合計ジャッジ時間: 187,038 ms
• ジャッジサーバーID （参考情報）: judge11 / judge15

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	2 ms; 4,376 KB
testcase_01	AC	2 ms; 4,380 KB
testcase_02	AC	1 ms; 4,380 KB
testcase_03	AC	2 ms; 4,376 KB
testcase_04	AC	30 ms; 4,380 KB
testcase_05	AC	30 ms; 4,380 KB
testcase_06	AC	31 ms; 4,376 KB
testcase_07	AC	27 ms; 4,380 KB
testcase_08	AC	29 ms; 4,380 KB
testcase_09	AC	214 ms; 4,380 KB
testcase_10	AC	358 ms; 4,380 KB
testcase_11	AC	440 ms; 4,380 KB
testcase_12	AC	305 ms; 4,376 KB
testcase_13	AC	461 ms; 4,380 KB
testcase_14	AC	490 ms; 4,380 KB
testcase_15	AC	556 ms; 4,380 KB
testcase_16	AC	663 ms; 4,376 KB
testcase_17	AC	537 ms; 4,380 KB
testcase_18	AC	654 ms; 4,472 KB
testcase_19	AC	157 ms; 4,380 KB
testcase_20	AC	323 ms; 4,376 KB
testcase_21	AC	517 ms; 4,380 KB
testcase_22	AC	384 ms; 4,376 KB
testcase_23	AC	473 ms; 4,380 KB
testcase_24	AC	411 ms; 4,380 KB
testcase_25	AC	342 ms; 4,380 KB
testcase_26	AC	448 ms; 4,380 KB
testcase_27	AC	390 ms; 4,380 KB
testcase_28	AC	349 ms; 4,380 KB
testcase_29	AC	750 ms; 4,808 KB
testcase_30	AC	572 ms; 4,580 KB
testcase_31	AC	534 ms; 4,380 KB
testcase_32	AC	497 ms; 4,380 KB
testcase_33	AC	383 ms; 4,376 KB
testcase_34	AC	378 ms; 4,380 KB
testcase_35	AC	282 ms; 4,380 KB
testcase_36	AC	298 ms; 4,380 KB
testcase_37	AC	204 ms; 4,376 KB
testcase_38	AC	78 ms; 4,380 KB
testcase_39	AC	590 ms; 4,560 KB
testcase_40	AC	770 ms; 4,380 KB
testcase_41	AC	605 ms; 4,380 KB
testcase_42	AC	552 ms; 4,376 KB
testcase_43	AC	552 ms; 4,380 KB
testcase_44	AC	486 ms; 4,380 KB
testcase_45	AC	395 ms; 4,380 KB
testcase_46	AC	365 ms; 4,376 KB
testcase_47	AC	411 ms; 4,376 KB
testcase_48	AC	328 ms; 4,376 KB
testcase_49	AC	658 ms; 4,580 KB
testcase_50	AC	631 ms; 4,488 KB
testcase_51	AC	648 ms; 4,380 KB
testcase_52	AC	628 ms; 4,380 KB
testcase_53	AC	562 ms; 4,380 KB
testcase_54	AC	487 ms; 4,376 KB
testcase_55	AC	488 ms; 4,376 KB
testcase_56	AC	634 ms; 4,380 KB
testcase_57	AC	502 ms; 4,380 KB
testcase_58	AC	239 ms; 4,376 KB
testcase_59	AC	639 ms; 4,520 KB
testcase_60	AC	903 ms; 4,376 KB
testcase_61	AC	652 ms; 4,376 KB
testcase_62	AC	643 ms; 4,380 KB
testcase_63	AC	590 ms; 4,376 KB
testcase_64	AC	549 ms; 4,376 KB
testcase_65	AC	382 ms; 4,376 KB
testcase_66	AC	429 ms; 4,380 KB
testcase_67	AC	469 ms; 4,376 KB
testcase_68	AC	220 ms; 4,376 KB
testcase_69	AC	94 ms; 4,376 KB
testcase_70	AC	89 ms; 4,376 KB
testcase_71	AC	89 ms; 4,380 KB
testcase_72	AC	59 ms; 4,376 KB
testcase_73	AC	90 ms; 4,376 KB
testcase_74	AC	1,324 ms; 4,500 KB
testcase_75	AC	1,105 ms; 4,508 KB
testcase_76	AC	1,497 ms; 4,532 KB
testcase_77	AC	1,472 ms; 4,476 KB
testcase_78	AC	1,866 ms; 4,536 KB
testcase_79	AC	2,116 ms; 4,604 KB
testcase_80	AC	1,931 ms; 4,504 KB
testcase_81	AC	2,496 ms; 4,540 KB
testcase_82	AC	2,038 ms; 4,804 KB
testcase_83	AC	2,486 ms; 5,748 KB
testcase_84	AC	1,343 ms; 4,540 KB
testcase_85	AC	1,720 ms; 4,608 KB
testcase_86	AC	997 ms; 4,552 KB
testcase_87	AC	1,520 ms; 4,536 KB
testcase_88	AC	1,944 ms; 4,536 KB
testcase_89	AC	1,675 ms; 4,508 KB
testcase_90	AC	1,382 ms; 4,552 KB
testcase_91	AC	1,714 ms; 4,540 KB
testcase_92	AC	1,456 ms; 4,536 KB
testcase_93	AC	2,257 ms; 4,536 KB
testcase_94	AC	2,895 ms; 6,528 KB
testcase_95	AC	2,560 ms; 6,200 KB
testcase_96	AC	2,264 ms; 5,216 KB
testcase_97	AC	1,971 ms; 4,992 KB
testcase_98	AC	1,588 ms; 4,712 KB
testcase_99	AC	1,382 ms; 4,648 KB
testcase_100	AC	1,262 ms; 4,536 KB
testcase_101	AC	898 ms; 4,536 KB
testcase_102	AC	749 ms; 4,532 KB
testcase_103	AC	264 ms; 4,544 KB
testcase_104	AC	2,638 ms; 6,476 KB
testcase_105	AC	2,618 ms; 5,244 KB
testcase_106	AC	2,581 ms; 5,168 KB
testcase_107	AC	2,303 ms; 4,616 KB
testcase_108	AC	1,976 ms; 4,680 KB
testcase_109	AC	2,050 ms; 4,664 KB
testcase_110	AC	1,806 ms; 4,504 KB
testcase_111	AC	1,654 ms; 4,572 KB
testcase_112	AC	1,305 ms; 4,504 KB
testcase_113	AC	1,106 ms; 4,540 KB
testcase_114	AC	2,891 ms; 6,504 KB
testcase_115	AC	2,671 ms; 6,236 KB
testcase_116	AC	2,956 ms; 5,216 KB
testcase_117	AC	2,863 ms; 4,952 KB
testcase_118	AC	2,091 ms; 4,680 KB
testcase_119	AC	1,723 ms; 4,672 KB
testcase_120	AC	2,473 ms; 4,540 KB
testcase_121	AC	1,939 ms; 4,508 KB
testcase_122	AC	1,699 ms; 4,508 KB
testcase_123	AC	1,058 ms; 4,576 KB
testcase_124	AC	2,665 ms; 6,352 KB
testcase_125	AC	2,386 ms; 5,212 KB
testcase_126	AC	2,565 ms; 5,116 KB
testcase_127	AC	2,451 ms; 4,828 KB
testcase_128	AC	2,267 ms; 4,620 KB
testcase_129	AC	1,768 ms; 4,716 KB
testcase_130	AC	1,908 ms; 4,552 KB
testcase_131	AC	1,683 ms; 4,504 KB
testcase_132	AC	2,138 ms; 4,504 KB
testcase_133	AC	2,365 ms; 4,532 KB
testcase_134	AC	192 ms; 4,376 KB
testcase_135	AC	180 ms; 4,380 KB
testcase_136	AC	132 ms; 4,376 KB
testcase_137	AC	123 ms; 4,380 KB
testcase_138	AC	96 ms; 4,380 KB
testcase_139	AC	95 ms; 4,380 KB
testcase_140	AC	102 ms; 4,376 KB
testcase_141	AC	112 ms; 4,376 KB
testcase_142	AC	96 ms; 4,380 KB
testcase_143	AC	107 ms; 4,376 KB
testcase_144	AC	283 ms; 4,376 KB
testcase_145	AC	313 ms; 4,376 KB
testcase_146	AC	418 ms; 4,376 KB
testcase_147	AC	488 ms; 4,376 KB
testcase_148	AC	637 ms; 4,376 KB
testcase_149	AC	440 ms; 4,380 KB
testcase_150	AC	463 ms; 4,376 KB
testcase_151	AC	610 ms; 4,376 KB
testcase_152	AC	539 ms; 4,376 KB
testcase_153	AC	565 ms; 4,380 KB
testcase_154	AC	131 ms; 4,380 KB
testcase_155	AC	467 ms; 4,568 KB
testcase_156	AC	368 ms; 4,552 KB
testcase_157	AC	475 ms; 4,544 KB
testcase_158	AC	478 ms; 4,520 KB
testcase_159	AC	489 ms; 4,488 KB
testcase_160	AC	603 ms; 4,504 KB
testcase_161	AC	753 ms; 4,536 KB
testcase_162	AC	396 ms; 4,572 KB
testcase_163	AC	295 ms; 4,508 KB
testcase_164	AC	324 ms; 4,488 KB
testcase_165	AC	1,750 ms; 4,540 KB
testcase_166	AC	1,328 ms; 4,528 KB
testcase_167	AC	1,792 ms; 4,540 KB
testcase_168	AC	1,361 ms; 4,616 KB
testcase_169	AC	1,585 ms; 4,492 KB
testcase_170	AC	2,552 ms; 4,540 KB
testcase_171	AC	3,669 ms; 4,504 KB
testcase_172	AC	1,744 ms; 4,516 KB
testcase_173	AC	2,053 ms; 4,544 KB
testcase_174	AC	2,548 ms; 4,504 KB
testcase_175	AC	409 ms; 4,508 KB

ソースコード

// 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';
}