#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 #pragma GCC target( \ "sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,avx2,tune=native") #pragma GCC optimize("Ofast") using namespace std; using ll = long long; constexpr double PI = 3.1415926535897932; template inline bool chmin(T& m, S q) { if (m > q) { m = q; return true; } else return false; } template inline bool chmax(T& m, const S q) { if (m < q) { m = q; return true; } else return false; } // 2 次元ベクトル template struct Vec2 { /* y 軸正は下方向 x 軸正は右方向 回転は時計回りが正(y 軸正を上と考えると反時計回りになる) */ using value_type = T; T y, x; constexpr inline Vec2() = default; constexpr Vec2(const T& arg_y, const T& arg_x) : y(arg_y), x(arg_x) {} inline Vec2(const Vec2&) = default; // コピー inline Vec2(Vec2&&) = default; // ムーブ inline Vec2& operator=(const Vec2&) = default; // 代入 inline Vec2& operator=(Vec2&&) = default; // ムーブ代入 template constexpr inline Vec2(const Vec2& v) : y((T)v.y), x((T)v.x) {} inline Vec2 operator+(const Vec2& rhs) const { return Vec2(y + rhs.y, x + rhs.x); } inline Vec2 operator+(const T& rhs) const { return Vec2(y + rhs, x + rhs); } inline Vec2 operator-(const Vec2& rhs) const { return Vec2(y - rhs.y, x - rhs.x); } template inline Vec2 operator*(const S& rhs) const { return Vec2(y * rhs, x * rhs); } inline Vec2 operator*(const Vec2& rhs) const { // x + yj とみなす return Vec2(x * rhs.y + y * rhs.x, x * rhs.x - y * rhs.y); } template inline Vec2 operator/(const S& rhs) const { assert(rhs != 0.0); return Vec2(y / rhs, x / rhs); } inline Vec2 operator/(const Vec2& rhs) const { // x + yj とみなす return (*this) * rhs.inv(); } inline Vec2& operator+=(const Vec2& rhs) { y += rhs.y; x += rhs.x; return *this; } inline Vec2& operator-=(const Vec2& rhs) { y -= rhs.y; x -= rhs.x; return *this; } template inline Vec2& operator*=(const S& rhs) const { y *= rhs; x *= rhs; return *this; } inline Vec2& operator*=(const Vec2& rhs) { return *this = (*this) * rhs; } inline Vec2& operator/=(const Vec2& rhs) { return *this = (*this) / rhs; } inline bool operator!=(const Vec2& rhs) const { return x != rhs.x || y != rhs.y; } inline bool operator==(const Vec2& rhs) const { return x == rhs.x && y == rhs.y; } inline void rotate(const double& rad) { *this = rotated(rad); } inline Vec2 rotated(const double& rad) const { return (*this) * rotation(rad); } static inline Vec2 rotation(const double& rad) { return Vec2(sin(rad), cos(rad)); } static inline Vec2 rotation_deg(const double& deg) { return rotation(PI * deg / 180.0); } inline Vec2 rounded() const { return Vec2(round(y), round(x)); } inline Vec2 inv() const { // x + yj とみなす const double norm_sq = l2_norm_square(); assert(norm_sq != 0.0); return Vec2(-y / norm_sq, x / norm_sq); } inline double l2_norm() const { return sqrt(x * x + y * y); } inline double l2_norm_square() const { return x * x + y * y; } inline T l1_norm() const { return std::abs(x) + std::abs(y); } inline double abs() const { return l2_norm(); } inline double phase() const { // [-PI, PI) のはず return atan2(y, x); } inline double phase_deg() const { // [-180, 180) のはず return phase() * (180.0 / PI); } }; template inline Vec2 operator*(const S& lhs, const Vec2& rhs) { return rhs * lhs; } template ostream& operator<<(ostream& os, const Vec2& vec) { os << vec.y << ' ' << vec.x; return os; } // 2 次元配列 template struct Board { array data; template constexpr inline auto& operator[](const Vec2& p) { return data[width * p.y + p.x]; } template constexpr inline const auto& operator[](const Vec2& p) const { return data[width * p.y + p.x]; } template constexpr inline auto& operator[](const initializer_list& p) { return data[width * *p.begin() + *(p.begin() + 1)]; } template constexpr inline const auto& operator[](const initializer_list& p) const { return data[width * *p.begin() + *(p.begin() + 1)]; } constexpr inline void Fill(const T& fill_value) { fill(data.begin(), data.end(), fill_value); } void Print() const { for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { cout << data[width * y + x] << " \n"[x == width - 1]; } } } }; inline double Time() { return static_cast( chrono::duration_cast( chrono::steady_clock::now().time_since_epoch()) .count()) * 1e-9; } // まず、全点間距離を求める // 次に、 2-opt する static constexpr auto N = 100; static constexpr auto M = 8; static constexpr auto NM = N + M; using Point = Vec2; // struct Point { // enum PointType { // kPlanet = 1, // kStation = 2, // }; // Vec2 yx; // PointType type; // short id; // array // }; auto rng = mt19937(); struct State { array points; double score; array order; double entire_best_score; array, M> entire_best_points; vector solution; void Read() { auto dummy = 0; cin >> dummy >> dummy; for (auto i = 0; i < N; i++) { cin >> points[i].y >> points[i].x; } for (auto i = 0; i < NM; i++) { order[i] = i; } order.back() = 0; entire_best_score = 1e30; } void Randomize() { for (auto i = 0; i < M; i++) { points[N + i].y = uniform_int_distribution<>(100, 900)(rng); points[N + i].x = uniform_int_distribution<>(100, 900)(rng); } shuffle(order.begin() + 1, order.end() - 1, rng); }; void PrintBest() const { for (auto i = 0; i < M; i++) { cout << entire_best_points[i] << endl; } cout << solution.size() << endl; for (const auto p : solution) { if (p < 100) { cout << "1 " << p + 1 << endl; } else { cout << "2 " << p - 99 << endl; } } } }; void TwoOpt(State& state) { auto& points = state.points; auto& order = state.order; auto distances = Board(); auto warshall_floyd_next = Board(); for (auto i = 0; i < NM; i++) { for (auto j = 0; j < NM; j++) { distances[{i, j}] = (points[i] - points[j]).l2_norm_square(); if (i < N) distances[{i, j}] *= 5.0f; if (j < N) distances[{i, j}] *= 5.0f; warshall_floyd_next[{i, j}] = j; } } for (int k = 0; k < NM; k++) for (int i = 0; i < NM; i++) for (int j = 0; j < NM; j++) if (distances[{i, k}] + distances[{k, j}] < distances[{i, j}]) { distances[{i, j}] = distances[{i, k}] + distances[{k, j}]; warshall_floyd_next[{i, j}] = warshall_floyd_next[{i, k}]; } // スコア初期化 auto& current_score = state.score; current_score = 0.0; for (auto i = 0; i < NM; i++) { current_score += distances[{order[i], order[i + 1]}]; } for (auto iteration = 0; iteration < 100000; iteration++) { // 0 -> a -> b ... // 0 <- d <- c ... // vvv // 0 -> a -> c ... // 0 <- d <- b ... auto idx_a = uniform_int_distribution<>(0, NM - 1)(rng); auto idx_c = uniform_int_distribution<>(0, NM - 1)(rng); if (abs(idx_a - idx_c) <= 1) continue; if (idx_a > idx_c) { swap(idx_a, idx_c); } const auto a = order[idx_a]; const auto b = order[idx_a + 1]; const auto c = order[idx_c]; const auto d = order[idx_c + 1]; const auto distance_before = distances[{a, b}] + distances[{c, d}]; const auto distance_after = distances[{a, c}] + distances[{b, d}]; const auto delta = distance_after - distance_before; if (delta < 0.0f) { reverse(order.begin() + idx_a + 1, order.begin() + idx_c + 1); current_score += delta; } } if (current_score >= state.entire_best_score) return; cerr << current_score << endl; state.entire_best_score = current_score; for (auto i = 0; i < M; i++) { state.entire_best_points[i] = {(short)points[N + i].y, (short)points[N + i].x}; } // 経路復元 // from https://zeosutt.hatenablog.com/entry/2015/05/05/045943 state.solution.clear(); for (auto idx_order = 0; idx_order < N + M; idx_order++) { const auto start = order[idx_order]; const auto goal = order[idx_order + 1]; for (short cur = start; cur != goal; cur = warshall_floyd_next[{cur, goal}]) state.solution.push_back(cur); } state.solution.push_back(0); } void Solve() { auto state = State(); state.Read(); auto t0 = Time(); auto outer_iteraiton = 0; while (Time() - t0 < 0.9) { state.Randomize(); TwoOpt(state); outer_iteraiton++; } cerr << "outer_iteraiton=" << outer_iteraiton << endl; state.PrintBest(); } int main() { Solve(); }