結果
| 問題 | No.5007 Steiner Space Travel |
| ユーザー |
 square1001
|
| 提出日時 | 2022-07-30 17:44:27 |
| 言語 | C++14 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 973 ms / 1,000 ms |
| コード長 | 8,430 bytes |
| コンパイル時間 | 1,226 ms |
| 実行使用メモリ | 5,164 KB |
| スコア | 8,786,873 |
| 最終ジャッジ日時 | 2022-07-30 17:45:42 |
| 合計ジャッジ時間 | 32,599 ms |
|
ジャッジサーバーID （参考情報） |
judge11 / judge15 |
| 純コード判定しない問題か言語 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 30 |
ソースコード
#include <cmath>#include <ctime>#include <vector>#include <cassert>#include <iostream>#include <algorithm>using namespace std;namespace myrandom {uint64_t seed = 123456789ULL;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<point2d> stations;vector<pair<int, int> > path;answer_container() : stations(vector<point2d>()), path(vector<pair<int, int> >()) {};answer_container(const std::vector<point2d>& stations_, const std::vector<pair<int, int> >& path_) : stations(stations_), path(path_) {};};class quadratic {public:// ax^2 + bx + cint 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 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<point2d>& P, const double TLE) {const int INF = 1012345678;int true_best_cost = INF;vector<int> true_best_perm;vector<point2d> true_best_qopt;int initial_t = clock();int loops = 0;while (double(clock() - initial_t) / CLOCKS_PER_SEC < TLE) {loops += 1;vector<int> 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<quadratic> qbasex(N);vector<quadratic> 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<quadratic> qx(M, quadratic(0, 0, 0));vector<quadratic> 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<point2d> qopt(M);for (int i = 0; i < M; i++) {qopt[i].x = qx[i].extreme_point();qopt[i].y = qy[i].extreme_point();}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<int> best_perm = perm;vector<point2d> 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 >= 30000 && ((P[perm[2 * vl - 2]] - P[perm[2 * vr - 2]]).norm() >= 640000 || (P[perm[2 * vl]] -P[perm[2 * vr]]).norm() >= 640000)));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 < 30000 ? 300000 : 150000) ? myrandom::next_int(0, 2) : 1);int rstate = ((P[vb] - P[vd]).norm() <= (iteration < 30000 ? 300000 : 150000) ? 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<quadratic> sqx = qx;vector<quadratic> 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 next_cost2 = 0;for (int i = 0; i < M; 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();}if (best_cost > next_cost) {best_cost = next_cost;best_perm = perm;best_qopt = qopt;cont = 0;}}if (cont == 3 * 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<pair<int, int> > 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<point2d>& P, const answer_container& answer) {int cost = 0;for (int i = 0; i < int(answer.path.size()) - 1; i++) {pair<int, int> u1 = answer.path[i];pair<int, int> 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<point2d> 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;}