#include #include #include #include #include #include using namespace std; namespace myrandom { uint64_t seed = 987654321ULL; uint64_t xorshift64() { seed ^= seed << 13; seed ^= seed >> 7; seed ^= seed << 17; return seed; } int next_int(int l, int r) { return l + int(xorshift64() % (r - l)); } double next_double(double l, double r) { return l + double(xorshift64()) / double(uint64_t(-1)) * (r - l); } } class point2d { public: int x, y; point2d() : x(0), y(0) {}; point2d(int x_, int y_) : x(x_), y(y_) {}; point2d& operator+=(const point2d& p) { x += p.x; y += p.y; return (*this); } point2d& operator-=(const point2d& p) { x -= p.x; y -= p.y; return (*this); } point2d operator+(const point2d& p) const { return point2d(*this) += p; } point2d operator-(const point2d& p) const { return point2d(*this) -= p; } int norm() const { return x * x + y * y; } }; class answer_container { public: vector stations; vector > path; answer_container() : stations(vector()), path(vector >()) {}; answer_container(const std::vector& stations_, const std::vector >& path_) : stations(stations_), path(path_) {}; }; class quadratic { public: // ax^2 + bx + c int a, b, c; quadratic() : a(0), b(0), c(0) {}; quadratic(int a_, int b_, int c_) : a(a_), b(b_), c(c_) {}; quadratic& operator+=(const quadratic& q) { a += q.a; b += q.b; c += q.c; return (*this); } quadratic& operator-=(const quadratic& q) { a -= q.a; b -= q.b; c -= q.c; return (*this); } quadratic operator+(const quadratic& q) const { return quadratic(*this) += q; } quadratic operator-(const quadratic& q) const { return quadratic(*this) -= q; } int get(int x) const { return a * x * x + b * x + c; } int extreme_point() const { if (a == 0) { assert(b == 0); return 0; } auto divide_and_round = [](int x, int y) { if (y > 0) { x = x * 2 + y; y *= 2; } else { x = x * -2 - y; y *= -2; } return (x >= 0 ? x / y : -((-x + y - 1) / y)); }; return divide_and_round(-b, 2 * a); } int extreme_value() const { int x = extreme_point(); return a * x * x + b * x + c; } }; answer_container solve(int N, int M, const std::vector& P, const double TLE) { const int INF = 1012345678; int true_best_cost = INF; vector true_best_perm; vector true_best_qopt; int initial_t = clock(); int loops = 0; while (double(clock() - initial_t) / CLOCKS_PER_SEC < TLE) { loops += 1; vector perm(2 * N + 1); for (int i = 0; i <= N; i++) { perm[2 * i] = i % N; if (i % N != 0) { swap(perm[2 * i], perm[myrandom::next_int(1, i + 1) * 2]); } } for (int i = 0; i < N; i++) { perm[2 * i + 1] = myrandom::next_int(-1, M); } vector qbasex(N); vector qbasey(N); for (int i = 0; i < N; i++) { qbasex[i] = quadratic(1, -2 * P[i].x, P[i].x * P[i].x); qbasey[i] = quadratic(1, -2 * P[i].y, P[i].y * P[i].y); } vector qx(M, quadratic(0, 0, 0)); vector qy(M, quadratic(0, 0, 0)); for (int i = 0; i < N; i++) { if (perm[2 * i + 1] != -1) { qx[perm[2 * i + 1]] += qbasex[perm[2 * i]] + qbasex[perm[2 * i + 2]]; qy[perm[2 * i + 1]] += qbasey[perm[2 * i]] + qbasey[perm[2 * i + 2]]; } } vector qopt(M); for (int i = 0; i < M; i++) { qopt[i].x = qx[i].extreme_point(); qopt[i].y = qy[i].extreme_point(); } vector qext(M); for (int i = 0; i < M; i++) { qext[i] = qx[i].get(qopt[i].x) + qy[i].get(qopt[i].y); } int cost1 = 0, cost2 = 0; for (int i = 0; i < N; i++) { if (perm[2 * i + 1] == -1) { cost1 += (P[perm[2 * i]] - P[perm[2 * i + 2]]).norm(); } } for (int i = 0; i < M; i++) { cost2 += qx[i].extreme_value(); cost2 += qy[i].extreme_value(); } int best_cost = cost1 * 25 + cost2 * 5; vector best_perm = perm; vector best_qopt = qopt; int iteration = 0; int cont = 0; while (double(clock() - initial_t) / CLOCKS_PER_SEC < TLE) { iteration += 1; cont += 1; int vl, vr; do { vl = myrandom::next_int(1, N + 1); vr = myrandom::next_int(1, N + 1); } while(vr - vl <= 0 || (iteration >= 50000 && ((P[perm[2 * vl - 2]] - P[perm[2 * vr - 2]]).norm() >= 810000 || (P[perm[2 * vl]] - P[perm[2 * vr]]).norm() >= 810000))); int idl = perm[2 * vl - 1]; int idr = perm[2 * vr - 1]; int va = perm[2 * vl - 2]; int vb = perm[2 * vl]; int vc = perm[2 * vr - 2]; int vd = perm[2 * vr]; int lstate = ((P[va] - P[vc]).norm() <= (iteration < 50000 ? 360000 : 202500) ? myrandom::next_int(0, 2) : 1); int rstate = ((P[vb] - P[vd]).norm() <= (iteration < 50000 ? 360000 : 202500) ? myrandom::next_int(0, 2) : 1); int next_idl = -1; if (lstate == 1) { int opt_idl = INF; for (int i = 0; i < M; i++) { int c1 = (P[va] - qopt[i]).norm(); int c2 = (P[vc] - qopt[i]).norm(); if (opt_idl > c1 + c2) { opt_idl = c1 + c2; next_idl = i; } } } int next_idr = -1; if (rstate == 1) { int opt_idr = INF; for (int i = 0; i < M; i++) { int c1 = (P[vb] - qopt[i]).norm(); int c2 = (P[vd] - qopt[i]).norm(); if (opt_idr > c1 + c2) { opt_idr = c1 + c2; next_idr = i; } } } vector sqx = qx; vector sqy = qy; int next_cost1 = cost1; if (idl != -1) { sqx[idl] -= (qbasex[va] + qbasex[vb]); sqy[idl] -= (qbasey[va] + qbasey[vb]); } else { next_cost1 -= (P[va] - P[vb]).norm(); } if (idr != -1) { sqx[idr] -= (qbasex[vc] + qbasex[vd]); sqy[idr] -= (qbasey[vc] + qbasey[vd]); } else { next_cost1 -= (P[vc] - P[vd]).norm(); } if (next_idl != -1) { sqx[next_idl] += (qbasex[va] + qbasex[vc]); sqy[next_idl] += (qbasey[va] + qbasey[vc]); } else { next_cost1 += (P[va] - P[vc]).norm(); } if (next_idr != -1) { sqx[next_idr] += (qbasex[vb] + qbasex[vd]); sqy[next_idr] += (qbasey[vb] + qbasey[vd]); } else { next_cost1 += (P[vb] - P[vd]).norm(); } int bit = (1 << idl) | (1 << idr) | (1 << next_idl) | (1 << next_idr); int next_cost2 = cost2; for (int i = 0; i < M; i++) { if (((bit >> i) & 1) == 1) { next_cost2 -= qext[i]; next_cost2 += sqx[i].extreme_value(); next_cost2 += sqy[i].extreme_value(); } } int cost = cost1 * 25 + cost2 * 5; int next_cost = next_cost1 * 25 + next_cost2 * 5; if (cost >= next_cost) { qx = sqx; qy = sqy; cost1 = next_cost1; cost2 = next_cost2; reverse(perm.begin() + (2 * vl), perm.begin() + (2 * vr - 1)); perm[2 * vl - 1] = next_idl; perm[2 * vr - 1] = next_idr; for (int i = 0; i < M; i++) { qopt[i].x = qx[i].extreme_point(); qopt[i].y = qy[i].extreme_point(); qext[i] = qx[i].get(qopt[i].x) + qy[i].get(qopt[i].y); } if (best_cost > next_cost) { best_cost = next_cost; best_perm = perm; best_qopt = qopt; cont = 0; } } if ((cont == N * N && best_cost >= true_best_cost + 150000) || cont == 6 * N * N) { break; } } if (true_best_cost > best_cost) { true_best_cost = best_cost; true_best_perm = best_perm; true_best_qopt = best_qopt; } cerr << "Loop #" << loops << ": Iterations = " << iteration << ", Cost = " << best_cost << endl; } vector > path; for (int i = 0; i <= 2 * N; i++) { if (i % 2 == 0) { path.push_back(make_pair(1, true_best_perm[i])); } else if (true_best_perm[i] != -1) { path.push_back(make_pair(2, true_best_perm[i])); } } return answer_container(true_best_qopt, path); } int calculate_cost(int N, int M, const vector& P, const answer_container& answer) { int cost = 0; for (int i = 0; i < int(answer.path.size()) - 1; i++) { pair u1 = answer.path[i]; pair u2 = answer.path[i + 1]; point2d p1 = (u1.first == 1 ? P[u1.second] : answer.stations[u1.second]); point2d p2 = (u2.first == 1 ? P[u2.second] : answer.stations[u2.second]); cost += (p1 - p2).norm() * (u1.first == 1 ? 5 : 1) * (u2.first == 1 ? 5 : 1); } return cost; } int main() { int N, M; cin >> N >> M; vector P(N); for (int i = 0; i < N; i++) { cin >> P[i].x >> P[i].y; } answer_container answer = solve(N, M, P, 0.94); for (int i = 0; i < M; i++) { cout << answer.stations[i].x << ' ' << answer.stations[i].y << endl; } cout << answer.path.size() << endl; for (int i = 0; i < int(answer.path.size()); i++) { cout << answer.path[i].first << ' ' << answer.path[i].second + 1 << endl; } int cost = calculate_cost(N, M, P, answer); cerr << "Final Cost = " << cost << endl; cerr << "Final Score = " << 1.0e+9 / (1.0e+3 + sqrt(cost)) << endl; return 0; }