fork download 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. using ll = long long;
  4. struct P { int x, y, id; };
  5. static inline ll dist2(const P &a, const P &b) {
  6.     ll dx = (ll)a.x - (ll)b.x;
  7.     ll dy = (ll)a.y - (ll)b.y;
  8.     return dx*dx + dy*dy;
  9. }
  10. unsigned long long hilbertKey(int x, int y, int bits = 21) {
  11.     unsigned long long hx = 0;
  12.     unsigned long long X = (unsigned long long)x, Y = (unsigned long long)y;
  13.     for (int i = bits - 1; i >= 0; --i) {
  14.         unsigned long long rx = (X >> i) & 1ull;
  15.         unsigned long long ry = (Y >> i) & 1ull;
  16.         unsigned long long idx = (rx << 1) | ry;
  17.         hx = (hx << 2) | idx;
  18.     }
  19.     return hx;
  20. }
  21. void two_opt(const vector<P> &pts, vector<int> &tour) {
  22.     int n = (int)tour.size();
  23.     if (n <= 3) return;
  24.     bool improved = true;
  25.     int iter = 0;
  26.     while (improved && iter < 200) {
  27.         improved = false; ++iter;
  28.         for (int i = 0; i < n; ++i) {
  29.             int i1 = (i + 1) % n;
  30.             P A = pts[tour[i]];
  31.             P B = pts[tour[i1]];
  32.             for (int j = i + 2; j < n; ++j) {
  33.                 if (i == 0 && j == n - 1) continue;
  34.                 int j1 = (j + 1) % n;
  35.                 P C = pts[tour[j]];
  36.                 P D = pts[tour[j1]];
  37.                 ll before = dist2(A, B) + dist2(C, D);
  38.                 ll after  = dist2(A, C) + dist2(B, D);
  39.                 if (after < before) {
  40.                     if (i1 <= j) reverse(tour.begin() + i1, tour.begin() + j + 1);
  41.                     improved = true;
  42.                     goto next_outer;
  43.                 }
  44.             }
  45.         }
  46.         next_outer: ;
  47.     }
  48. }
  49. double tour_length(const vector<P> &pts, const vector<int> &tour) {
  50.     int n = (int)tour.size();
  51.     if (n == 0) return 0.0;
  52.     double s = 0.0;
  53.     for (int i = 0; i < n; ++i) {
  54.         int j = (i + 1) % n;
  55.         s += sqrt((double)dist2(pts[tour[i]], pts[tour[j]]));
  56.     }
  57.     return s;
  58. }
  59. int main() {
  60.     ios::sync_with_stdio(false);
  61.     cin.tie(nullptr);
  62.     int N, K;
  63.     if (!(cin >> N >> K)) return 0;
  64.     vector<P> pts(N);
  65.     for (int i = 0; i < N; ++i) { cin >> pts[i].x >> pts[i].y; pts[i].id = i; }
  66.     vector<pair<unsigned long long,int>> ord; ord.reserve(N);
  67.     int bits = 21;
  68.     for (int i = 0; i < N; ++i) ord.emplace_back(hilbertKey(pts[i].x, pts[i].y, bits), i);
  69.     sort(ord.begin(), ord.end());
  70.     vector<int> perm; perm.reserve(N);
  71.     for (auto &pr : ord) perm.push_back(pr.second);
  72.     vector<int> sizes(K, N / K);
  73.     for (int i = 0; i < N % K; ++i) sizes[i]++;
  74.     vector<vector<int>> clusters(K);
  75.     int idx = 0;
  76.     for (int k = 0; k < K; ++k) {
  77.         int s = sizes[k];
  78.         clusters[k].reserve(s);
  79.         for (int t = 0; t < s; ++t) clusters[k].push_back(perm[idx++]);
  80.     }
  81.     for (int k = 0; k < K; ++k) two_opt(pts, clusters[k]);
  82.     vector<pair<double,double>> cent(K);
  83.     for (int k = 0; k < K; ++k) {
  84.         double sx = 0, sy = 0;
  85.         for (int id : clusters[k]) { sx += pts[id].x; sy += pts[id].y; }
  86.         int sz = max(1, (int)clusters[k].size());
  87.         cent[k] = { sx / sz, sy / sz };
  88.     }
  89.     int max_rounds = 1000;
  90.     for (int round = 0; round < max_rounds; ++round) {
  91.         vector<double> lens(K);
  92.         for (int k = 0; k < K; ++k) lens[k] = tour_length(pts, clusters[k]);
  93.         int worst = max_element(lens.begin(), lens.end()) - lens.begin();
  94.         bool moved = false;
  95.         if ((int)clusters[worst].size() <= 1) break;
  96.         int m = (int)clusters[worst].size();
  97.         vector<pair<double,int>> edge_contrib; edge_contrib.reserve(m);
  98.         for (int i = 0; i < m; ++i) {
  99.             int j = (i + 1) % m;
  100.             double d = sqrt((double)dist2(pts[clusters[worst][i]], pts[clusters[worst][j]]));
  101.             edge_contrib.emplace_back(d, i);
  102.         }
  103.         sort(edge_contrib.rbegin(), edge_contrib.rend());
  104.         int tries = min(8, (int)edge_contrib.size());
  105.         for (int tt = 0; tt < tries && !moved; ++tt) {
  106.             int pos = edge_contrib[tt].second;
  107.             int u = clusters[worst][pos];
  108.             double best_delta = 1e100;
  109.             int best_k = -1;
  110.             int best_insert_pos = -1;
  111.             for (int k = 0; k < K; ++k) {
  112.                 if (k == worst) continue;
  113.                 int mk = (int)clusters[k].size();
  114.                 if (mk == 0) continue;
  115.                 double best_ins = 1e100; int bi = -1;
  116.                 for (int i = 0; i < mk; ++i) {
  117.                     int ni = (i + 1) % mk;
  118.                     double before = sqrt((double)dist2(pts[clusters[k][i]], pts[clusters[k][ni]]));
  119.                     double after  = sqrt((double)dist2(pts[clusters[k][i]], pts[u])) + sqrt((double)dist2(pts[u], pts[clusters[k][ni]]));
  120.                     double delta = after - before;
  121.                     if (delta < best_ins) { best_ins = delta; bi = ni; }
  122.                 }
  123.                 int prev = (pos - 1 + m) % m;
  124.                 int next = (pos + 1) % m;
  125.                 double rem_before = sqrt((double)dist2(pts[clusters[worst][prev]], pts[clusters[worst][pos]]))
  126.                                   + sqrt((double)dist2(pts[clusters[worst][pos]], pts[clusters[worst][next]]));
  127.                 double rem_after = sqrt((double)dist2(pts[clusters[worst][prev]], pts[clusters[worst][next]]));
  128.                 double delta_total = best_ins + (rem_after - rem_before);
  129.                 if (delta_total < best_delta) { best_delta = delta_total; best_k = k; best_insert_pos = bi; }
  130.             }
  131.             if (best_k != -1 && best_delta < -1e-7) {
  132.                 clusters[worst].erase(clusters[worst].begin() + pos);
  133.                 clusters[best_k].insert(clusters[best_k].begin() + best_insert_pos, u);
  134.                 two_opt(pts, clusters[worst]);
  135.                 two_opt(pts, clusters[best_k]);
  136.                 for (int t = 0; t < 2; ++t) {
  137.                     int kk = (t == 0 ? worst : best_k);
  138.                     double sx = 0, sy = 0;
  139.                     for (int id2 : clusters[kk]) { sx += pts[id2].x; sy += pts[id2].y; }
  140.                     int sz = max(1, (int)clusters[kk].size());
  141.                     cent[kk] = { sx / sz, sy / sz };
  142.                 }
  143.                 moved = true;
  144.             }
  145.         }
  146.         if (!moved) break;
  147.     }
  148.     for (int k = 0; k < K; ++k) {
  149.         cout << clusters[k].size();
  150.         for (int id : clusters[k]) cout << ' ' << pts[id].id + 1;
  151.         cout << "\n";
  152.     }
  153.     return 0;
  154. }
Success #stdin #stdout 0s 5320KB
comments (?)
stdin 
3 2
0 0
3 0
3 1
stdout 
1 1
2 2 3
https://ideone.com/utFIq0
language:
C++ (gcc 8.3)
created:
8 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!