#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class JsonDumper { struct KeyValue { std::string key; std::string value; }; std::vector _items; bool dump_at_end = false; public: JsonDumper(bool dump_at_end_ = false) : dump_at_end(dump_at_end_) {} ~JsonDumper() { if (dump_at_end) std::cout << dump() << std::endl; } void set_dump_at_end() { dump_at_end = true; } void operator()(const std::string &key, const std::string &value) { _items.push_back(KeyValue{key, "\"" + value + "\""}); } template void operator()(const std::string &key, T value) { _items.push_back(KeyValue{key, std::to_string(value)}); } std::string dump() const { std::string ret = "{\n"; if (!_items.empty()) { for (const auto &[k, v] : _items) ret += " \"" + k + "\": " + v + ",\n"; ret.erase(ret.end() - 2); } ret += "}"; return ret; } } jdump; #define ALL(x) (x).begin(), (x).end() #define FOR(i, begin, end) for (int i = (begin), i##_end_ = (end); i < i##_end_; i++) #define IFOR(i, begin, end) for (int i = (end)-1, i##_begin_ = (begin); i >= i##_begin_; i--) #define REP(i, n) FOR(i, 0, n) #define IREP(i, n) IFOR(i, 0, n) template bool chmax(T &m, const T q) { return m < q ? (m = q, true) : false; } template bool chmin(T &m, const T q) { return m > q ? (m = q, true) : false; } int floor_lg(long long x) { return x <= 0 ? -1 : 63 - __builtin_clzll(x); } template T1 floor_div(T1 num, T2 den) { return (num > 0 ? num / den : -((-num + den - 1) / den)); } template std::pair operator+(const std::pair &l, const std::pair &r) { return std::make_pair(l.first + r.first, l.second + r.second); } template std::pair operator-(const std::pair &l, const std::pair &r) { return std::make_pair(l.first - r.first, l.second - r.second); } template std::vector sort_unique(std::vector vec) { sort(vec.begin(), vec.end()), vec.erase(unique(vec.begin(), vec.end()), vec.end()); return vec; } template int arglb(const std::vector &v, const T &x) { return std::distance(v.begin(), std::lower_bound(v.begin(), v.end(), x)); } template int argub(const std::vector &v, const T &x) { return std::distance(v.begin(), std::upper_bound(v.begin(), v.end(), x)); } template IStream &operator>>(IStream &is, std::vector &vec) { for (auto &v : vec) is >> v; return is; } template OStream &operator<<(OStream &os, const std::vector &vec); template OStream &operator<<(OStream &os, const std::array &arr); template OStream &operator<<(OStream &os, const std::unordered_set &vec); template OStream &operator<<(OStream &os, const std::pair &pa); template OStream &operator<<(OStream &os, const std::deque &vec); template OStream &operator<<(OStream &os, const std::set &vec); template OStream &operator<<(OStream &os, const std::multiset &vec); template OStream &operator<<(OStream &os, const std::unordered_multiset &vec); template OStream &operator<<(OStream &os, const std::pair &pa); template OStream &operator<<(OStream &os, const std::map &mp); template OStream &operator<<(OStream &os, const std::unordered_map &mp); template OStream &operator<<(OStream &os, const std::tuple &tpl); template OStream &operator<<(OStream &os, const std::vector &vec) { os << '['; for (auto v : vec) os << v << ','; os << ']'; return os; } template OStream &operator<<(OStream &os, const std::array &arr) { os << '['; for (auto v : arr) os << v << ','; os << ']'; return os; } template std::istream &operator>>(std::istream &is, std::tuple &tpl) { std::apply([&is](auto &&...args) { ((is >> args), ...); }, tpl); return is; } template OStream &operator<<(OStream &os, const std::tuple &tpl) { os << '('; std::apply([&os](auto &&...args) { ((os << args << ','), ...); }, tpl); return os << ')'; } template OStream &operator<<(OStream &os, const std::unordered_set &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::deque &vec) { os << "deq["; for (auto v : vec) os << v << ','; os << ']'; return os; } template OStream &operator<<(OStream &os, const std::set &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::multiset &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::unordered_multiset &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::pair &pa) { return os << '(' << pa.first << ',' << pa.second << ')'; } template OStream &operator<<(OStream &os, const std::map &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; } template OStream &operator<<(OStream &os, const std::unordered_map &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; } #ifdef HITONANODE_LOCAL const std::string COLOR_RESET = "\033[0m", BRIGHT_GREEN = "\033[1;32m", BRIGHT_RED = "\033[1;31m", BRIGHT_CYAN = "\033[1;36m", NORMAL_CROSSED = "\033[0;9;37m", RED_BACKGROUND = "\033[1;41m", NORMAL_FAINT = "\033[0;2m"; #define dbg(x) \ std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " \ << __FILE__ << COLOR_RESET << std::endl #define dbgif(cond, x) \ ((cond) ? std::cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ \ << ") " << __FILE__ << COLOR_RESET << std::endl \ : std::cerr) #else #define dbg(x) 0 #define dbgif(cond, x) 0 #endif #ifdef BENCHMARK #define dump_onlinejudge(x) 0 struct setenv { setenv() { jdump.set_dump_at_end(); } } setenv_; #else #define dump_onlinejudge(x) (std::cout << (x) << std::endl) #endif #include #include uint32_t rand_int() // XorShift random integer generator { static uint32_t x = 123456789, y = 362436069, z = 521288629, w = 88675123; uint32_t t = x ^ (x << 11); x = y; y = z; z = w; return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); } double rand_double() { return (double)rand_int() / UINT32_MAX; } #include #include #include #include #include #include #include #include #include #include #include template ::max() / 2, int INVALID = -1> struct shortest_path { int V, E; bool single_positive_weight; T wmin, wmax; std::vector> tos; std::vector head; std::vector> edges; void build_() { if (int(tos.size()) == E and int(head.size()) == V + 1) return; tos.resize(E); head.assign(V + 1, 0); for (const auto &e : edges) ++head[std::get<0>(e) + 1]; for (int i = 0; i < V; ++i) head[i + 1] += head[i]; auto cur = head; for (const auto &e : edges) { tos[cur[std::get<0>(e)]++] = std::make_pair(std::get<1>(e), std::get<2>(e)); } } shortest_path(int V = 0) : V(V), E(0), single_positive_weight(true), wmin(0), wmax(0) {} void add_edge(int s, int t, T w) { assert(0 <= s and s < V); assert(0 <= t and t < V); edges.emplace_back(s, 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); } void add_bi_edge(int u, int v, T w) { add_edge(u, v, w); add_edge(v, u, w); } std::vector dist; std::vector prev; using Pque = std::priority_queue, std::vector>, std::greater>>; template void dijkstra(int s, int t = INVALID) { assert(0 <= s and s < V); build_(); dist.assign(V, INF); prev.assign(V, INVALID); dist[s] = 0; Heap pq; pq.emplace(0, s); while (!pq.empty()) { T d; int v; std::tie(d, v) = pq.top(); pq.pop(); if (t == v) return; if (dist[v] < d) continue; for (int e = head[v]; e < head[v + 1]; ++e) { const auto &nx = tos[e]; T dnx = d + nx.second; if (dist[nx.first] > dnx) { dist[nx.first] = dnx, prev[nx.first] = v; pq.emplace(dnx, nx.first); } } } } void dijkstra_vquad(int s, int t = INVALID) { assert(0 <= s and s < V); build_(); dist.assign(V, INF); prev.assign(V, INVALID); dist[s] = 0; std::vector fixed(V, false); while (true) { int r = INVALID; T dr = INF; for (int i = 0; i < V; i++) { if (!fixed[i] and dist[i] < dr) r = i, dr = dist[i]; } if (r == INVALID or r == t) break; fixed[r] = true; int nxt; T dx; for (int e = head[r]; e < head[r + 1]; ++e) { std::tie(nxt, dx) = tos[e]; if (dist[nxt] > dist[r] + dx) dist[nxt] = dist[r] + dx, prev[nxt] = r; } } } bool bellman_ford(int s, int nb_loop) { assert(0 <= s and s < V); build_(); dist.assign(V, INF), prev.assign(V, INVALID); dist[s] = 0; for (int l = 0; l < nb_loop; l++) { bool upd = false; for (int v = 0; v < V; v++) { if (dist[v] == INF) continue; for (int e = head[v]; e < head[v + 1]; ++e) { const auto &nx = tos[e]; T dnx = dist[v] + nx.second; if (dist[nx.first] > dnx) dist[nx.first] = dnx, prev[nx.first] = v, upd = true; } } if (!upd) return true; } return false; } void spfa(int s) { assert(0 <= s and s < V); build_(); dist.assign(V, INF); prev.assign(V, INVALID); dist[s] = 0; std::deque q; std::vector in_queue(V); q.push_back(s), in_queue[s] = 1; while (!q.empty()) { int now = q.front(); q.pop_front(), in_queue[now] = 0; for (int e = head[now]; e < head[now + 1]; ++e) { const auto &nx = tos[e]; T dnx = dist[now] + nx.second; int nxt = nx.first; if (dist[nxt] > dnx) { dist[nxt] = dnx; if (!in_queue[nxt]) { if (q.size() and dnx < dist[q.front()]) { // Small label first optimization q.push_front(nxt); } else { q.push_back(nxt); } prev[nxt] = now, in_queue[nxt] = 1; } } } } } void zero_one_bfs(int s, int t = INVALID) { assert(0 <= s and s < V); build_(); dist.assign(V, INF), prev.assign(V, INVALID); dist[s] = 0; std::vector q(V * 4); int ql = V * 2, qr = V * 2; q[qr++] = s; while (ql < qr) { int v = q[ql++]; if (v == t) return; for (int e = head[v]; e < head[v + 1]; ++e) { const auto &nx = tos[e]; T dnx = dist[v] + nx.second; if (dist[nx.first] > dnx) { dist[nx.first] = dnx, prev[nx.first] = v; if (nx.second) { q[qr++] = nx.first; } else { q[--ql] = nx.first; } } } } } void dial(int s, int t = INVALID) { assert(0 <= s and s < V); build_(); dist.assign(V, INF), prev.assign(V, INVALID); dist[s] = 0; std::vector>> q(wmax + 1); q[0].emplace_back(s, dist[s]); int ninq = 1; int cur = 0; T dcur = 0; for (; ninq; ++cur, ++dcur) { if (cur == wmax + 1) cur = 0; while (!q[cur].empty()) { int v = q[cur].back().first; T dnow = q[cur].back().second; q[cur].pop_back(), --ninq; if (v == t) return; if (dist[v] < dnow) continue; for (int e = head[v]; e < head[v + 1]; ++e) { const auto &nx = tos[e]; T dnx = dist[v] + nx.second; if (dist[nx.first] > dnx) { dist[nx.first] = dnx, prev[nx.first] = v; int nxtcur = cur + int(nx.second); if (nxtcur >= int(q.size())) nxtcur -= q.size(); q[nxtcur].emplace_back(nx.first, dnx), ++ninq; } } } } } bool dag_solver(int s) { assert(0 <= s and s < V); build_(); dist.assign(V, INF), prev.assign(V, INVALID); dist[s] = 0; std::vector indeg(V, 0); std::vector q(V * 2); int ql = 0, qr = 0; q[qr++] = s; while (ql < qr) { int now = q[ql++]; for (int e = head[now]; e < head[now + 1]; ++e) { const auto &nx = tos[e]; ++indeg[nx.first]; if (indeg[nx.first] == 1) q[qr++] = nx.first; } } ql = qr = 0; q[qr++] = s; while (ql < qr) { int now = q[ql++]; for (int e = head[now]; e < head[now + 1]; ++e) { const auto &nx = tos[e]; --indeg[nx.first]; if (dist[nx.first] > dist[now] + nx.second) dist[nx.first] = dist[now] + nx.second, prev[nx.first] = now; if (indeg[nx.first] == 0) q[qr++] = nx.first; } } return *max_element(indeg.begin(), indeg.end()) == 0; } std::vector retrieve_path(int goal) const { assert(int(prev.size()) == V); assert(0 <= goal and goal < V); if (dist[goal] == INF) return {}; std::vector ret{goal}; while (prev[goal] != INVALID) { goal = prev[goal]; ret.push_back(goal); } std::reverse(ret.begin(), ret.end()); return ret; } void solve(int s, int t = INVALID) { if (wmin >= 0) { if (single_positive_weight) { zero_one_bfs(s, t); } else if (wmax <= 10) { dial(s, t); } else { if ((long long)V * V < (E << 4)) { dijkstra_vquad(s, t); } else { dijkstra(s, t); } } } else { bellman_ford(s, V); } } std::vector> floyd_warshall() { build_(); std::vector> dist2d(V, std::vector(V, INF)); for (int i = 0; i < V; i++) { dist2d[i][i] = 0; for (const auto &e : edges) { int s = std::get<0>(e), t = std::get<1>(e); dist2d[s][t] = std::min(dist2d[s][t], std::get<2>(e)); } } for (int k = 0; k < V; k++) { for (int i = 0; i < V; i++) { if (dist2d[i][k] == INF) continue; for (int j = 0; j < V; j++) { if (dist2d[k][j] == INF) continue; dist2d[i][j] = std::min(dist2d[i][j], dist2d[i][k] + dist2d[k][j]); } } } return dist2d; } void to_dot(std::string filename = "shortest_path") const { std::ofstream ss(filename + ".DOT"); ss << "digraph{\n"; build_(); for (int i = 0; i < V; i++) { for (int e = head[i]; e < head[i + 1]; ++e) { ss << i << "->" << tos[e].first << "[label=" << tos[e].second << "];\n"; } } ss << "}\n"; ss.close(); return; } }; #include #include #include #include #include #include #include #include template std::vector> mst_edges(const DistanceMatrix &dist) { using T = decltype((*dist.adjacents(0).begin()).second); if (dist.n() <= 1) return {}; std::vector dp(dist.n(), std::numeric_limits::max()); std::vector prv(dist.n(), -1); std::vector used(dist.n()); std::vector> ret(dist.n() - 1); for (int t = 0; t < dist.n(); ++t) { int x = std::min_element(dp.cbegin(), dp.cend()) - dp.cbegin(); dp.at(x) = std::numeric_limits::max(); used.at(x) = 1; if (t > 0) ret.at(t - 1) = {prv.at(x), x}; for (auto [y, len] : dist.adjacents(x)) { if (!used.at(y) and len < dp.at(y)) dp.at(y) = len, prv.at(y) = x; } } return ret; } template auto calc_lkh_alpha(const DistanceMatrix &dist) { using T = decltype((*dist.adjacents(0).begin()).second); std::vector> to(dist.n()); for (auto [s, t] : mst_edges(dist)) { to.at(s).push_back(t); to.at(t).push_back(s); } std::vector ret(dist.n(), std::vector(dist.n())); for (int s = 0; s < dist.n(); ++s) { auto rec = [&](auto &&self, int now, int prv, T hi) -> void { ret.at(s).at(now) = dist(s, now) - hi; for (int nxt : to.at(now)) { if (nxt == prv) continue; self(self, nxt, now, std::max(hi, dist(now, nxt))); } }; rec(rec, s, -1, T()); } int best_one = -1; T longest_2nd_nearest = T(); for (int one = 0; one < dist.n(); ++one) { if (to.at(one).size() != 1) continue; const int ng = to.at(one).front(); bool found = false; T second_nearest = T(); for (auto [v, len] : dist.adjacents(one)) { if (v == ng) continue; if (!found) { found = true, second_nearest = len; } else if (len < second_nearest) { second_nearest = len; } } if (found and (best_one < 0 or second_nearest > longest_2nd_nearest)) best_one = one, longest_2nd_nearest = second_nearest; } if (best_one != -1) { for (auto [v, len] : dist.adjacents(best_one)) { if (v == to.at(best_one).front()) continue; ret.at(best_one).at(v) = ret.at(v).at(best_one) = len - longest_2nd_nearest; } } return ret; } template class csr_distance_matrix { int _rows = 0; std::vector begins; std::vector> vals; public: csr_distance_matrix() : csr_distance_matrix({}, 0) {} csr_distance_matrix(const std::vector> &init, int rows) : _rows(rows), begins(rows + 1) { std::vector degs(rows); for (const auto &p : init) ++degs.at(std::get<0>(p)); for (int i = 0; i < rows; ++i) begins.at(i + 1) = begins.at(i) + degs.at(i); vals.resize(init.size(), std::make_pair(-1, T())); for (auto [i, j, w] : init) vals.at(begins.at(i + 1) - (degs.at(i)--)) = {j, w}; } void apply_pi(const std::vector &pi) { for (int i = 0; i < n(); ++i) { for (auto &[j, d] : adjacents(i)) d += pi.at(i) + pi.at(j); } } int n() const noexcept { return _rows; } struct adjacents_sequence { csr_distance_matrix *ptr; int from; using iterator = typename std::vector>::iterator; iterator begin() { return std::next(ptr->vals.begin(), ptr->begins.at(from)); } iterator end() { return std::next(ptr->vals.begin(), ptr->begins.at(from + 1)); } }; struct const_adjacents_sequence { const csr_distance_matrix *ptr; const int from; using const_iterator = typename std::vector>::const_iterator; const_iterator begin() const { return std::next(ptr->vals.cbegin(), ptr->begins.at(from)); } const_iterator end() const { return std::next(ptr->vals.cbegin(), ptr->begins.at(from + 1)); } }; adjacents_sequence adjacents(int from) { return {this, from}; } const_adjacents_sequence adjacents(int from) const { return {this, from}; } const_adjacents_sequence operator()(int from) const { return {this, from}; } }; template auto build_adjacent_info(const DistanceMatrix &dist, int sz) { using T = decltype((*dist.adjacents(0).begin()).second); const std::vector> alpha = calc_lkh_alpha(dist); std::vector> adjacent_edges; std::vector> candidates; for (int i = 0; i < dist.n(); ++i) { candidates.clear(); for (auto [j, d] : dist.adjacents(i)) { if (i != j) candidates.emplace_back(alpha.at(i).at(j), d, j); } const int final_sz = std::min(sz, candidates.size()); std::nth_element(candidates.begin(), candidates.begin() + final_sz, candidates.end()); candidates.resize(final_sz); std::sort(candidates.begin(), candidates.end(), [&](const auto &l, const auto &r) { return std::get<1>(l) < std::get<1>(r); }); for (auto [alpha, dij, j] : candidates) adjacent_edges.emplace_back(i, j, dij); } return csr_distance_matrix(adjacent_edges, dist.n()); } #include #include #include #include template struct SymmetricTSP { DistanceMatrix dist; Adjacents adjacents; using T = decltype((*dist.adjacents(0).begin()).second); struct Solution { T cost; std::vector path; template friend OStream &operator<<(OStream &os, const Solution &x) { os << "[cost=" << x.cost << ", path=("; for (int i : x.path) os << i << ","; return os << x.path.front() << ")]"; } }; T eval(const Solution &sol) const { T ret = T(); int now = sol.path.back(); for (int nxt : sol.path) ret += dist(now, nxt), now = nxt; return ret; } SymmetricTSP(const DistanceMatrix &distance_matrix, const Adjacents &adjacents) : dist(distance_matrix), adjacents(adjacents) {} Solution nearest_neighbor(int init) const { if (n() == 0) return {T(), {}}; int now = init; std::vector ret{now}, alive(n(), 1); T len = T(); ret.reserve(n()); alive.at(now) = 0; while (int(ret.size()) < n()) { int nxt = -1; for (int i = 0; i < n(); ++i) { if (alive.at(i) and (nxt < 0 or dist(now, i) < dist(now, nxt))) nxt = i; } ret.push_back(nxt); alive.at(nxt) = 0; len += dist(now, nxt); now = nxt; } len += dist(ret.back(), ret.front()); return Solution{len, ret}; } void two_opt(Solution &sol) const { static std::vector v_to_i; v_to_i.resize(n()); for (int i = 0; i < n(); ++i) v_to_i.at(sol.path.at(i)) = i; while (true) { bool updated = false; for (int i = 0; i < n() and !updated; ++i) { const int u = sol.path.at(i), nxtu = sol.path.at(modn(i + 1)); const T dunxtu = dist(u, nxtu); for (auto [v, duv] : adjacents(u)) { if (duv >= dunxtu) break; int j = v_to_i.at(v), nxtv = sol.path.at(modn(j + 1)); T diff = duv + dist(nxtu, nxtv) - dunxtu - dist(v, nxtv); if (diff < 0) { sol.cost += diff; int l, r; if (modn(j - i) < modn(i - j)) { l = modn(i + 1), r = j; } else { l = modn(j + 1), r = i; } while (l != r) { std::swap(sol.path.at(l), sol.path.at(r)); v_to_i.at(sol.path.at(l)) = l; v_to_i.at(sol.path.at(r)) = r; l = modn(l + 1); if (l == r) break; r = modn(r - 1); } updated = true; break; } } if (updated) break; for (auto [nxtv, dnxtunxtv] : adjacents(nxtu)) { if (dnxtunxtv >= dunxtu) break; int j = modn(v_to_i.at(nxtv) - 1), v = sol.path.at(j); T diff = dist(u, v) + dnxtunxtv - dunxtu - dist(v, nxtv); if (diff < 0) { sol.cost += diff; int l, r; if (modn(j - i) < modn(i - j)) { l = modn(i + 1), r = j; } else { l = modn(j + 1), r = i; } while (l != r) { std::swap(sol.path.at(l), sol.path.at(r)); v_to_i.at(sol.path.at(l)) = l; v_to_i.at(sol.path.at(r)) = r; l = modn(l + 1); if (l == r) break; r = modn(r - 1); } updated = true; break; } } } if (!updated) break; } } bool three_opt(Solution &sol) const { static std::vector v_to_i; v_to_i.resize(n()); for (int i = 0; i < n(); ++i) v_to_i.at(sol.path.at(i)) = i; auto check_uvw_order = [](int u, int v, int w) { int i = v_to_i.at(u); int j = v_to_i.at(v); int k = v_to_i.at(w); if (i < j and j < k) return true; if (j < k and k < i) return true; if (k < i and i < j) return true; return false; }; auto rev = [&](const int u, const int v) -> void { int l = v_to_i.at(u), r = v_to_i.at(v); while (l != r) { std::swap(sol.path.at(l), sol.path.at(r)); l = modn(l + 1); if (l == r) break; r = modn(r - 1); } }; static int i = 0; for (int nseen = 0; nseen < n(); ++nseen, i = modn(i + 1)) { const int u = sol.path.at(modn(i - 1)), nxtu = sol.path.at(i); const T dunxtu = dist(u, nxtu); for (const auto &[nxtv, dunxtv] : adjacents(u)) { if (dunxtv >= dunxtu) break; const int v = sol.path.at(modn(v_to_i.at(nxtv) - 1)); const T dvnxtv = dist(v, nxtv); for (const auto &[nxtw, dvnxtw] : adjacents(v)) { if (nxtw == nxtv or nxtw == nxtu) continue; if (dunxtv + dvnxtw >= dunxtu + dvnxtv) break; const int w = sol.path.at(modn(v_to_i.at(nxtw) - 1)); if (!check_uvw_order(u, v, w)) continue; const T current = dunxtu + dvnxtv + dist(w, nxtw); if (T diff = dunxtv + dist(w, nxtu) + dvnxtw - current; diff < T()) { sol.cost += diff; rev(nxtu, v); rev(nxtv, w); rev(nxtw, u); return true; } } for (const auto &[w, dvw] : adjacents(v)) { if (dunxtv + dvw >= dunxtu + dvnxtv) break; if (!check_uvw_order(u, v, w)) continue; const int nxtw = sol.path.at(modn(v_to_i.at(w) + 1)); const T current = dunxtu + dvnxtv + dist(w, nxtw); if (T diff = dunxtv + dvw + dist(nxtu, nxtw) - current; diff < T()) { sol.cost += diff; rev(nxtw, u); rev(nxtv, w); return true; } } } for (const auto &[nxtv, dnxtunxtv] : adjacents(nxtu)) { if (dnxtunxtv >= dunxtu) break; const int v = sol.path.at(modn(v_to_i.at(nxtv) - 1)); const T dvnxtv = dist(v, nxtv); for (const auto &[nxtw, dvnxtw] : adjacents(v)) { const int w = sol.path.at(modn(v_to_i.at(nxtw) - 1)); if (dnxtunxtv + dvnxtw >= dunxtu + dvnxtv) break; if (!check_uvw_order(u, v, w)) continue; const T current = dunxtu + dvnxtv + dist(w, nxtw); if (T diff = dist(u, w) + dnxtunxtv + dvnxtw - current; diff < T()) { sol.cost += diff; rev(nxtu, v); rev(nxtw, u); return true; } } } for (const auto &[v, duv] : adjacents(u)) { if (duv >= dunxtu) break; const int nxtv = sol.path.at(modn(v_to_i.at(v) + 1)); const T dvnxtv = dist(v, nxtv); for (const auto &[nxtw, dnxtvnxtw] : adjacents(nxtv)) { const int w = sol.path.at(modn(v_to_i.at(nxtw) - 1)); if (duv + dnxtvnxtw >= dunxtu + dvnxtv) break; if (!check_uvw_order(u, v, w)) continue; const T current = dunxtu + dvnxtv + dist(w, nxtw); if (T diff = duv + dist(nxtu, w) + dnxtvnxtw - current; diff < T()) { sol.cost += diff; rev(nxtu, v); rev(nxtv, w); return true; } } } } return false; } template bool double_bridge(Solution &sol, Rng &rng) const { if (n() < 8) return false; std::vector &p = sol.path; int rand_rot = std::uniform_int_distribution(0, n() - 1)(rng); std::rotate(p.begin(), p.begin() + rand_rot, p.end()); static std::array arr; for (int &y : arr) y = std::uniform_int_distribution(2, n() - 6)(rng); std::sort(arr.begin(), arr.end()); const int i = arr.at(0), j = arr.at(1) + 2, k = arr.at(2) + 4; static std::array diffs; for (int d = 0; d < 2; ++d) { int u = p.at(n() - 1), nxtu = p.at(0); int v = p.at(i - 1), nxtv = p.at(i); int w = p.at(j - 1), nxtw = p.at(j); int x = p.at(k - 1), nxtx = p.at(k); diffs.at(d) = dist(u, nxtu) + dist(v, nxtv) + dist(w, nxtw) + dist(x, nxtx); if (d == 1) break; std::reverse(p.begin(), p.begin() + i); std::reverse(p.begin() + i, p.begin() + j); std::reverse(p.begin() + j, p.begin() + k); std::reverse(p.begin() + k, p.end()); } sol.cost += diffs.at(1) - diffs.at(0); return true; } int n() const noexcept { return dist.n(); } int modn(int x) const noexcept { if (x < 0) return x + n(); if (x >= n()) return x - n(); return x; } }; template class dense_distance_matrix { int _n; std::vector _d; public: dense_distance_matrix(const std::vector> &distance_vecvec) : _n(distance_vecvec.size()) { _d.reserve(n() * n()); for (const auto &vec : distance_vecvec) _d.insert(end(_d), begin(vec), end(vec)); } template void apply_pi(const std::vector &pi) { for (int i = 0; i < n(); ++i) { for (int j = 0; j < n(); ++j) _d.at(i * n() + j) += pi.at(i) + pi.at(j); } } int n() const noexcept { return _n; } T dist(int s, int t) const { return _d.at(s * n() + t); } T operator()(int s, int t) const { return dist(s, t); } struct adjacents_sequence { const dense_distance_matrix *ptr; int from; struct iterator { const dense_distance_matrix *ptr; int from; int to; iterator operator++() { return {ptr, from, to++}; } std::pair operator*() const { return {to, ptr->dist(from, to)}; } bool operator!=(const iterator &rhs) const { return to != rhs.to or ptr != rhs.ptr or from != rhs.from; } }; iterator begin() const { return iterator{ptr, from, 0}; } iterator end() const { return iterator{ptr, from, ptr->n()}; } }; adjacents_sequence adjacents(int from) const { return {this, from}; } }; struct Solver { static constexpr int N = 100, M = 8; static constexpr int alpha = 5; static constexpr int UB = 1000; std::vector X, Y; std::vector C, D; std::vector> direct_distances; int calc_direct_dist(int i, int j) const { int dx = X.at(i) - X.at(j); int dy = Y.at(i) - Y.at(j); return (dx * dx + dy * dy) * (alpha * alpha); } int calc_s2s_dist(int p, int q) const { int dx = C.at(p) - C.at(q); int dy = D.at(p) - D.at(q); return (dx * dx + dy * dy); } int calc_p2s_dist(int i, int p) const { int dx = X.at(i) - C.at(p); int dy = Y.at(i) - D.at(p); return (dx * dx + dy * dy) * alpha; } using Matrix = dense_distance_matrix; using TSP = SymmetricTSP>; void three_opt(const TSP &tsp, TSP::Solution &sol) { do { tsp.two_opt(sol); } while (tsp.three_opt(sol)); // three_opt() は一度 3-opt による改善に成功した時点で return } TSP::Solution sol; Solver(const std::vector &a, const std::vector &b) : X(a), Y(b), C(M), D(M) { direct_distances.assign(N, std::vector(N)); for (int i = 0; i < N; ++i) { direct_distances.at(i).at(i) = 0; for (int j = i + 1; j < N; ++j) { direct_distances.at(i).at(j) = direct_distances.at(j).at(i) = calc_direct_dist(i, j); } } for (int k = 0; k < N; ++k) { for (int i = 0; i < N; ++i) { for (int j = 0; j < N; ++j) { direct_distances.at(i).at(j) = std::min( direct_distances.at(i).at(j), direct_distances.at(i).at(k) + direct_distances.at(k).at(j)); } } } for (int i = 0; i < M; ++i) { C.at(i) = rand_int() % (UB + 1); D.at(i) = rand_int() % (UB + 1); } auto best_mat = gen_best_mat(); csr_distance_matrix adj = build_adjacent_info(best_mat, 30); auto tsp = SymmetricTSP(best_mat, adj); sol = tsp.nearest_neighbor(0); dbg(sol.cost); three_opt(tsp, sol); dbg(sol.cost); } Matrix gen_best_mat() const { auto ret = direct_distances; std::vector pqwf(M, std::vector(M)); for (int p = 0; p < M; ++p) { for (int q = 0; q < M; ++q) { pqwf.at(p).at(q) = calc_s2s_dist(p, q); } } for (int k = 0; k < M; ++k) { for (int p = 0; p < M; ++p) { for (int q = 0; q < M; ++q) { pqwf.at(p).at(q) = std::min(pqwf.at(p).at(q), pqwf.at(p).at(k) + pqwf.at(k).at(q)); } } } for (int i = 0; i < N; ++i) { for (int j = 0; j < i; ++j) { ret.at(i).at(j) = calc_direct_dist(i, j); for (int p = 0; p < M; ++p) { for (int q = 0; q < M; ++q) { ret.at(i).at(j) = std::min( ret.at(i).at(j), calc_p2s_dist(i, p) + pqwf.at(p).at(q) + calc_p2s_dist(j, q)); } } ret.at(j).at(i) = ret.at(i).at(j); assert(ret.at(i).at(j) >= 0); } } return Matrix(ret); } void step(int m, int dmax) { const int dx = rand_int() % (dmax * 2 + 1) - dmax; const int dy = rand_int() % (dmax * 2 + 1) - dmax; C.at(m) += dx; D.at(m) += dy; auto tmpsol = sol; auto best_mat = gen_best_mat(); csr_distance_matrix adj = build_adjacent_info(best_mat, 30); const auto tsp = SymmetricTSP(best_mat, adj); tmpsol.cost = tsp.eval(tmpsol); tsp.two_opt(tmpsol); if (tmpsol.cost < sol.cost) { dbg(tmpsol.cost); sol = tmpsol; } else { C.at(m) -= dx; D.at(m) -= dy; } } void finalize() { auto best_mat = gen_best_mat(); csr_distance_matrix adj = build_adjacent_info(best_mat, 30); const auto tsp = SymmetricTSP(best_mat, adj); three_opt(tsp, sol); } void exp_step() { const int m = rand_int() % M; const int cm = C.at(m), dm = D.at(m); int dx = 0, dy = 0; while (!dx and !dy) { dx = rand_int() % 31 - 15; dy = rand_int() % 31 - 15; } bool init = true; while (true) { C.at(m) += dx; D.at(m) += dy; auto tmpsol = sol; auto best_mat = gen_best_mat(); csr_distance_matrix adj = build_adjacent_info(best_mat, 30); const auto tsp = SymmetricTSP(best_mat, adj); tmpsol.cost = tsp.eval(tmpsol); tsp.two_opt(tmpsol); if (tmpsol.cost < sol.cost) { dbg(tmpsol.cost); sol = tmpsol; if (init) { init = false; } else { dx *= 2; dy *= 2; } } else { C.at(m) -= dx; D.at(m) -= dy; dbg(std::make_tuple(m, dx, dy)); break; } } } int score() const { return round(1e9 / (1000 + sqrt(sol.cost))); } std::vector> dump_solution() const { shortest_path graph(N + M); for (int i = 0; i < N; ++i) { for (int j = i + 1; j < N; ++j) graph.add_bi_edge(i, j, calc_direct_dist(i, j)); for (int p = 0; p < M; ++p) graph.add_bi_edge(i, N + p, calc_p2s_dist(i, p)); } for (int p = 0; p < M; ++p) { for (int q = p + 1; q < M; ++q) graph.add_bi_edge(N + p, N + q, calc_s2s_dist(p, q)); } auto path = sol.path; std::rotate(path.begin(), std::find(path.begin(), path.end(), 0), path.end()); assert(path.front() == 0); path.push_back(0); dbg(path); std::vector> ret; ret.emplace_back(1, 0); for (int c = 1; c < (int)path.size(); ++c) { int prv = path.at(c - 1), cur = path.at(c); graph.solve(prv, cur); auto p = graph.retrieve_path(cur); for (int d = 1; d < (int)p.size(); ++d) { int v = p.at(d); if (v < N) ret.emplace_back(1, v); else ret.emplace_back(2, v - N); } } return ret; } std::string build_dump_str() const { std::string ret; for (int p = 0; p < M; ++p) ret += std::to_string(C.at(p)) + " " + std::to_string(D.at(p)) + "\n"; auto vec = dump_solution(); ret += std::to_string(vec.size()) + "\n"; for (auto [k, v] : vec) { ret += std::to_string(k) + " " + std::to_string(v + 1) + "\n"; } return ret; } }; using namespace std; using lint = long long; using pint = std::pair; using plint = std::pair; struct fast_ios { fast_ios() { std::cin.tie(nullptr), std::ios::sync_with_stdio(false), std::cout << std::fixed << std::setprecision(20); }; } fast_ios_; int main(int argc, char *argv[]) { int X = 0; if (argc >= 2) { X = std::stoi(argv[1]); } int N, M; cin >> N >> M; vector A(N), B(N); REP(i, N) cin >> A.at(i) >> B.at(i); dbg(make_tuple(N, M, A, B)); dbg(A); dbg(B); Solver solver(A, B); for (int d = 350; d > 0; d -= 8) { for (int m = 0; m < M; ++m) solver.step(m, d); } solver.finalize(); dbg(solver.sol.cost); dbg(solver.score()); // dump_onlinejudge("solution"); auto sol = solver.dump_solution(); dbg(sol); dump_onlinejudge(solver.build_dump_str()); jdump("score", solver.score()); }