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