namespace atcoder {} #ifdef LOCAL #define dbg(x) cerr << __LINE__ << " : " << #x << " = " << (x) << endl; #else #define NDEBUG #define dbg(x) true; #pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") #endif #ifdef GTEST #include #endif #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef PERF #include #endif using namespace std; using namespace atcoder; #define fast_io \ ios_base::sync_with_stdio(false); \ cin.tie(0); \ cout.tie(0); #define ll long long int #define rep(i, n) for (int i = 0; i < (int)(n); i++) #define reps(i, n) for (int i = 1; i <= (int)(n); i++) #define REP(i, n) for (int i = n - 1; i >= 0; i--) #define REPS(i, n) for (int i = n; i > 0; i--) #define MOD (long long int)(1e9 + 7) #define INF (int)(1e9) #define LINF (long long int)(1e18) #define chmax(a, b) a = (((a) < (b)) ? (b) : (a)) #define chmin(a, b) a = (((a) > (b)) ? (b) : (a)) #define all(v) v.begin(), v.end() typedef pair Pii; typedef pair Pll; constexpr double PI = acos(-1); #ifdef NDEBUG #define CHECK(v1, op, v2) #else #define CHECK(v1, op, v2) \ if (!((v1)op(v2))) { \ cerr << "ERROR:" << (v1) << " " << (v2) << endl; \ assert((v1)op(v2)); \ } #endif long double nCr(const int n, const int r) { long double ret = 1; rep(t, r) { ret *= (n - t); ret /= (r - t); } return ret; } template string to_string(const vector& vec) { string ret = ""; rep(i, vec.size()) { ret += vec[i].to_string(); if (i + 1 != vec.size()) { ret += ","; } } return ret; } template ostream& operator<<(ostream& os, const vector& vec) { os << to_string(vec); return os; } uint32_t xorshift() { static uint32_t x = 12345789; static uint32_t y = 362436069; static uint32_t z = 521288629; static uint32_t w = 88675123; uint32_t t; t = x ^ (x << 11); x = y; y = z; z = w; w ^= t ^ (t >> 8) ^ (w >> 19); return w; } int rand(const uint32_t l, const uint32_t r) { return xorshift() % (r - l) + l; } uint32_t rand_other_than(const uint32_t l, const uint32_t r, const uint32_t other) { const uint32_t num = rand(l, r - 1); return num + (num >= other); } template const T& rand_vec(const vector& vec) { assert(vec.size() > 0); return vec[rand(0, vec.size())]; } template void shuffle(vector& vec) { rep(l, (int)vec.size() - 1) { const int idx = rand(l, vec.size()); swap(vec[idx], vec[l]); } } /* start */ vector PARAMS = {}; class Timer { chrono::system_clock::time_point _start, _end; ll _sum = 0, _count = 0; public: void start() { _start = chrono::system_clock::now(); } void stop() { _end = chrono::system_clock::now(); } void add() { const chrono::system_clock::time_point now = chrono::system_clock::now(); _sum += static_cast( chrono::duration_cast(now - _start).count()); _count++; } ll sum() const { return _sum / 1000; } int count() const { return _count; } string average() const { if (_count == 0) { return "NaN"; } return to_string(_sum / 1000 / _count); } void reset() { _start = chrono::system_clock::now(); _sum = 0; _count = 0; } inline int ms() const { const chrono::system_clock::time_point now = chrono::system_clock::now(); return static_cast( chrono::duration_cast(now - _start).count() / 1000); } inline int ns() const { const chrono::system_clock::time_point now = chrono::system_clock::now(); return static_cast( chrono::duration_cast(now - _start).count()); } }; #ifdef LOCAL struct Timers : unordered_map { friend ostream& operator<<(ostream& os, const Timers& timers) { for (const auto& pa : timers) { os << pa.first << " time: " << pa.second.sum() / 1000 << " count: " << pa.second.count() << endl; } return os; } }; #else struct Timers { struct Dummy { void start() const {} void add() const {} }; Dummy dummy; const Dummy& operator[](const std::string& str) { return dummy; } friend ostream& operator<<(ostream& os, const Timers& timers) { return os; } }; #endif Timers global_timers; // NBD = neighborhood template struct SimulatedAnnealing { public: SimulatedAnnealing(double end_time, double start_temp, double end_temp) : end_time_(end_time), inv_time_(1.0 / end_time_), cur_time_(0.0), start_temp_(start_temp), end_temp_(end_temp), diff_temp_(end_temp - start_temp), cur_temp_(start_temp) { assert(start_temp >= end_temp); timer_.start(); } Solution run(const Solution& initial_sol, NBDGenerator& nbd_generator) { Solution sol(initial_sol); int acc = 0; int g = 0; for (;; ++g) { if ((g & 0xf) == 0) { #ifdef SA_MAX_G if (g >= SA_MAX_G) { break; } #endif UpdateTime(); UpdateTemp(); if (cur_time_ > end_time_) { break; } } const auto nbd = nbd_generator.Generate(&sol, g); constexpr int mod = 1'000'000'000; const int r = rand(0, mod); if (static_cast(r) / mod < exp(-nbd->Diff() / cur_temp_)) { acc++; nbd->Update(&sol); } else { nbd->Restore(&sol); } } cerr << "g,acc: " << g << " " << acc << endl; return sol; } private: double end_time_, inv_time_, cur_time_; double start_temp_, end_temp_, diff_temp_, cur_temp_; Timer timer_; void UpdateTime() { cur_time_ = timer_.ms(); } void UpdateTemp() { cur_temp_ = max(0.0, start_temp_ - cur_time_ * diff_temp_ * inv_time_); } }; /* start */ template class Graph_T { public: struct Edge { Point to; Dist_T dist; Edge(const Point _to, const Dist_T _dist) : to(_to), dist(_dist) {} }; Graph_T() {} explicit Graph_T(const int n) : n_(n), adjs_(n), dists_(n, vector(n)) {} int N() const { return n_; } Dist_T Dist(const Point& from, const Point& to) const { return dists_[from][to]; } void AddEdge(const Point a, const Point b, const Dist_T dist) { adjs_[a].emplace_back(b, dist); dists_[a][b] = dist; } const vector& Adj(const Point a) const { return adjs_[a]; } void RemoveEdge(const Point a, const Point b) { { const auto itr = std::find_if( all(adjs_[a]), [&](const Edge& edge) { return edge.to == b; }); assert(itr != adjs_[a].end()); adjs_[a].erase(itr); dists_[a][b] = INF; } } private: int n_; vector> adjs_; vector> dists_; }; enum Dir { // y-1, y+1, x-1, x+1 kU, kD, kL, kR }; struct Pos { int idx_; Pos() {} explicit Pos(const int _idx) : idx_(_idx) {} Pos(int _x, int _y) : idx_(_y * N + _x) { assert(N > 0); } int X() const { return pos_2_x[*this]; } int Y() const { return pos_2_y[*this]; } int Idx() const { return idx_; } operator int() const { return Idx(); } operator size_t() const { return Idx(); } const vector& Adj() const { return adj_poses[*this]; } const vector& AdjDirs() const { return adj_dirs[*this]; } int Manhattan(const Pos& other) const { return abs(X() - other.X()) + abs(Y() - other.Y()); } int Euclid2(const Pos& other) const { const int dx = X() - other.X(); const int dy = Y() - other.Y(); return dx * dx + dy * dy; } bool Move(const Dir dir) { if (move_to[dir][*this].IsDummy()) { return false; } else { *this = move_to[dir][*this]; return true; } } bool operator<(const Pos& other) const { return this->Idx() < other.Idx(); } bool operator==(const Pos& other) const { return this->Idx() == other.Idx(); } bool operator!=(const Pos& other) const { return this->Idx() != other.Idx(); } friend ostream& operator<<(ostream& os, const Pos& pos) { os << pos.X() << " " << pos.Y(); return os; } bool IsDummy() const { return this->idx_ < 0; } static Pos Dummy() { Pos p; p.idx_ = -1; return p; } static Pos TryCreate(const int x, const int y, const int z) { if (y < 0 || y >= N || x < 0 || x >= N) { return Pos::Dummy(); } else { return Pos(x, y); } } static void StaticInit(const int n) { N = n; N2 = N * N; N4 = N2 * N2; adj_poses = vector>(N2); adj_dirs = vector>(N2); move_to = vector>(4, vector(N2, Pos::Dummy())); pos_2_y.resize(N2); pos_2_x.resize(N2); rep(y, N) { rep(x, N) { const Pos p{x, y}; pos_2_y[p] = y; pos_2_x[p] = x; for (int dy = -1; dy <= 1; ++dy) { for (int dx = -1; dx <= 1; ++dx) { if (abs(dy) + abs(dx) != 1) continue; const int adj_y = y + dy; const int adj_x = x + dx; if (adj_y < 0 || adj_y >= N) continue; if (adj_x < 0 || adj_x >= N) continue; adj_poses[p].emplace_back(adj_x, adj_y); Dir dir; if (dy == -1) { dir = kU; } else if (dy == 1) { dir = kD; } else if (dx == -1) { dir = kL; } else if (dx == 1) { dir = kR; } else { assert(false); } move_to[dir][p] = {adj_x, adj_y}; adj_dirs[p].emplace_back(dir); } } } } } static int N, N2, N4; static vector> adj_poses; static vector> adj_dirs; static vector> move_to; static vector pos_2_x, pos_2_y; }; vector> Pos::adj_poses; vector> Pos::adj_dirs; vector> Pos::move_to; vector Pos::pos_2_x, Pos::pos_2_y; int Pos::N, Pos::N2, Pos::N4; /* start */ constexpr int kPlanetType = 1; constexpr int kStationType = 2; struct Solution { using Graph = Graph_T; Solution() { orders.resize(N); std::iota(all(orders), 0); stations.resize(M); rep(i, M) { stations[i] = Pos{rand(0, Pos::N + 1), rand(0, Pos::N + 1)}; } } static void StaticInit(const int _N, const int _M, const vector& _planets, const Graph& _graph) { N = _N; M = _M; planets = _planets; graph = _graph; } // [l, r)のreverse int DiffSwap(const int l, const int r) const { // l-1, l, r-1, r // l-1, r-1, l, r const int i0 = (l - 1 + N) % N; const int i1 = l % N; const int i2 = (r - 1 + N) % N; const int i3 = r % N; const int prev_dist = graph.Dist(orders[i0], orders[i1]) + graph.Dist(orders[i2], orders[i3]); const int next_dist = graph.Dist(orders[i0], orders[i2]) + graph.Dist(orders[i1], orders[i3]); return next_dist - prev_dist; } void Swap(int l, int r) { if (l > r) { swap(l, r); } reverse(orders.begin() + l, orders.begin() + r); } friend ostream& operator<<(ostream& os, const Solution& sol) { int start_i = 0; for (start_i = 0; start_i < N; start_i++) { if (sol.orders[start_i] == 0) { break; } } #if 0 // その順序固定だとして、10^6通りの駅について、「その駅を通ったほうが嬉しい」辺のbitsetは何通りあるか? using BS = bitset<100>; unordered_set S; rep(y, 1001) { rep(x, 1001) { const Pos station(x, y); BS bs = 0; rep(i, N) { const int next_i = (i + 1) % N; const auto p0 = planets[sol.orders[i]]; const auto p1 = planets[sol.orders[next_i]]; const int default_dist = graph.Dist(sol.orders[i], sol.orders[next_i]); const int station_dist = p0.Euclid2(station) * alpha + station.Euclid2(p1) * alpha; bs[i] = station_dist <= default_dist; } S.insert(bs); } } cerr << endl; cerr << "size: " << S.size() << endl; #endif vector> ans; ans.emplace_back(kPlanetType, sol.orders[start_i]); rep(add_i, N) { const int i = (start_i + add_i) % N; const int next_i = (i + 1) % N; const Pos& from_planet = planets[sol.orders[i]]; const Pos& to_planet = planets[sol.orders[next_i]]; // stationを検討 int min_dist = Solution::graph.Dist(sol.orders[i], sol.orders[next_i]); vector best_stations; // i,jの順に経由 rep(i, M) { const Pos& station_i = sol.stations[i]; rep(j, M) { if (i == j) continue; const Pos& station_j = sol.stations[j]; const int dist = from_planet.Euclid2(station_i) * alpha + station_i.Euclid2(station_j) + station_j.Euclid2(to_planet) * alpha; if (dist < min_dist) { min_dist = dist; best_stations = {i, j}; } } } // iのみ経由 rep(i, M) { const Pos& station = sol.stations[i]; const int dist = from_planet.Euclid2(station) * alpha + station.Euclid2(to_planet) * alpha; if (dist < min_dist) { min_dist = dist; best_stations = {i}; } } if (best_stations.size() == 0) { // 惑星to惑星 int from = sol.orders[i]; int to = sol.orders[next_i]; vector rev_planets; while (from != to) { rev_planets.emplace_back(to); rep(i, N) { if (i == to) continue; if (graph.Dist(from, i) + graph.Dist(i, to) == graph.Dist(from, to)) { to = i; break; } } } REP(i, rev_planets.size()) { ans.emplace_back(kPlanetType, rev_planets[i]); } } else { for (const auto station : best_stations) { ans.emplace_back(kStationType, station); } ans.emplace_back(kPlanetType, sol.orders[next_i]); } } rep(i, M) { os << sol.stations[i] << endl; } os << ans.size() << endl; rep(i, ans.size()) { const auto& [t, r] = ans[i]; os << t << " " << r + 1; if (i + 1 != (int)ans.size()) { os << endl; } } int dist_sum = 0; rep(i, (int)ans.size() - 1) { const auto& [t0, r0] = ans[i]; const auto& [t1, r1] = ans[(i + 1) % ans.size()]; const Pos& from = (t0 == kPlanetType ? planets[r0] : sol.stations[r0]); const Pos& to = (t1 == kPlanetType ? planets[r1] : sol.stations[r1]); dist_sum += from.Euclid2(to) * (t0 == kPlanetType ? alpha : 1) * (t1 == kPlanetType ? alpha : 1); } cerr << "score: " << std::round(1e9 / (1000 + sqrt(dist_sum))) << endl; return os; } static vector planets; static int N, M; static Graph graph; static constexpr int alpha = 5; static constexpr int alpha2 = 25; vector orders; vector stations; }; int Solution::N; int Solution::M; vector Solution::planets; Graph_T Solution::graph; /* start */ struct NBD { int Diff() const { return diff_; } virtual void Update(Solution*) const = 0; void Restore(Solution*) const {} int diff_; }; struct NBD_swap : public NBD { NBD_swap() {} NBD_swap(Solution* const sol, const int l, const int r) { l_ = l; r_ = r; diff_ = sol->DiffSwap(l, r); } void Update(Solution* sol) const override { sol->Swap(l_, r_); } int l_; int r_; }; struct NBDGenerator { NBDGenerator() {} NBD* Generate(Solution* const sol, const int g) { // swap const int l = rand(0, Solution::N); const int r = rand_other_than(0, Solution::N, l); nbd_swap0 = NBD_swap(sol, l, r); return &nbd_swap0; } NBD_swap nbd_swap0; }; class Solver { public: using Graph = Graph_T; Solver(istream& is) { is >> N >> M; Pos::StaticInit(1001); planets.reserve(N); rep(i, N) { int a, b; is >> a >> b; planets.emplace_back(a, b); } graph = InitGraph(); Solution::StaticInit(N, M, planets, graph); } Graph InitGraph() const { vector> dists(N, vector(N, INF)); rep(i, N) { rep(j, N) { dists[i][j] = alpha2 * planets[i].Euclid2(planets[j]); } } rep(k, N) { rep(i, N) { rep(j, N) { chmin(dists[i][j], dists[i][k] + dists[k][j]); } } } Graph graph(N); rep(i, N) { rep(j, N) { const int dist = dists[i][j]; graph.AddEdge(i, j, dist); } } return graph; } using BS = bitset<100>; unordered_map> CalcBestBSs( const vector& orders) const { // その順序固定だとして、10^6通りの駅について、「その駅を通ったほうが嬉しい」辺のbitsetを求めて、 // 各bitsetごとに、距離の改善量が最も大きいものを保持する // TODO: 駅to駅を考慮 // TODO: intにするかも unordered_map> mp; rep(y, 1001) { rep(x, 1001) { const Pos station(x, y); BS bs = 0; ll diff_sum = 0; rep(i, N) { const int next_i = (i + 1) % N; const auto p0 = planets[orders[i]]; const auto p1 = planets[orders[next_i]]; const int default_dist = graph.Dist(orders[i], orders[next_i]); const int station_dist = p0.Euclid2(station) * alpha + station.Euclid2(p1) * alpha; const ll diff = default_dist - station_dist; if (diff > 0) { bs[i] = true; diff_sum += diff; } } if (diff_sum == 0) continue; if (mp[bs].second < diff_sum) { mp[bs] = {Pos{x, y}, diff_sum}; } } } cerr << "size: " << mp.size() << endl; return mp; } Solution DPStations(const Solution& order_sol) const { const auto orders = order_sol.orders; const auto mp = CalcBestBSs(orders); // dp[j][k] := i個目まで見て、j個駅を置いて、bit_countがk個 using P = pair; using IDX = tuple; vector dp(M + 1, vector(N + 1, P{BS{0}, (int)-INF})); vector prev(mp.size() + 1, vector(M + 1, vector(101, IDX{0, 0, 0}))); dp[0][0].second = 0; vector pos_candidates; pos_candidates.reserve(mp.size()); int i = -1; for (const auto& [bs, pos_diff] : mp) { ++i; const int bs_count = bs.count(); const auto& [pos, diff] = pos_diff; pos_candidates.emplace_back(pos); prev[i + 1] = prev[i]; REP(j, M) { const auto next_j = j + 1; REP(k, N + 1 - bs_count) { if (dp[j][k].second < 0) continue; if ((dp[j][k].first & bs) != 0) continue; const auto next_bs = (dp[j][k].first | bs); const auto next_k = next_bs.count(); const auto next_val = dp[j][k].second + diff; if (dp[next_j][next_k].second < next_val) { dp[next_j][next_k].first = next_bs; dp[next_j][next_k].second = next_val; prev[i + 1][next_j][next_k] = make_tuple(i, j, k); } } } } IDX best_idx{0, 0, 0}; ll best_diff = 0; rep(j, M + 1) { rep(k, N + 1) { const auto diff = dp[j][k].second; if (best_diff < diff) { best_diff = diff; best_idx = make_tuple((int)mp.size(), j, k); } } } // 復元 vector stations; IDX idx = best_idx; while (idx != make_tuple(0, 0, 0)) { const auto& [i, j, k] = idx; const auto& prev_idx = prev[i][j][k]; const auto& [prev_i, prev_j, prev_k] = prev_idx; if (prev_j != j) { stations.emplace_back(pos_candidates[prev_i]); } idx = prev_idx; } cerr << "station size: " << stations.size() << endl; while (stations.size() < M) { stations.emplace_back(Pos{500, 500}); } auto sol = order_sol; sol.stations = stations; return sol; } Solution Solve(const int time_limit) const { const Solution initial_sol; SimulatedAnnealing sa(time_limit - 500, 1, 0.01); NBDGenerator nbd_generator; const auto order_sol = sa.run(initial_sol, nbd_generator); const auto final_sol = DPStations(order_sol); return final_sol; } private: int N, M; vector planets; Graph graph; static constexpr int alpha = 5; static constexpr int alpha2 = 25; }; int main(int argc, char* argv[]) { fast_io; if (argc >= 2) { int idx = 0; for (int i = 1; i < argc; ++i) { PARAMS[idx++] = std::stod(argv[i]); } } Timer timer; timer.start(); Solver solver(cin); auto sol = solver.Solve(900 - timer.ms()); cout << sol << endl; return 0; }