結果

問題 No.5007 Steiner Space Travel
ユーザー square1001square1001
提出日時 2022-07-30 17:59:06
言語 C++14
(gcc 12.3.0 + boost 1.83.0)
結果
AC  
実行時間 969 ms / 1,000 ms
コード長 8,875 bytes
コンパイル時間 1,279 ms
実行使用メモリ 6,952 KB
スコア 8,822,957
最終ジャッジ日時 2022-07-30 17:59:45
合計ジャッジ時間 32,543 ms
ジャッジサーバーID
（参考情報） 		judge14 / judge16
純コード判定しない問題か言語
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 945 ms
4,900 KB
testcase_01 AC 945 ms
4,904 KB
testcase_02 AC 946 ms
6,948 KB
testcase_03 AC 946 ms
4,904 KB
testcase_04 AC 945 ms
4,904 KB
testcase_05 AC 945 ms
4,904 KB
testcase_06 AC 946 ms
4,908 KB
testcase_07 AC 945 ms
6,952 KB
testcase_08 AC 968 ms
6,948 KB
testcase_09 AC 947 ms
4,900 KB
testcase_10 AC 945 ms
4,900 KB
testcase_11 AC 947 ms
6,948 KB
testcase_12 AC 968 ms
6,952 KB
testcase_13 AC 969 ms
4,900 KB
testcase_14 AC 957 ms
4,904 KB
testcase_15 AC 953 ms
4,908 KB
testcase_16 AC 947 ms
4,900 KB
testcase_17 AC 945 ms
4,904 KB
testcase_18 AC 968 ms
4,900 KB
testcase_19 AC 947 ms
6,948 KB
testcase_20 AC 962 ms
4,900 KB
testcase_21 AC 946 ms
4,904 KB
testcase_22 AC 945 ms
6,948 KB
testcase_23 AC 947 ms
4,900 KB
testcase_24 AC 945 ms
4,904 KB
testcase_25 AC 968 ms
4,904 KB
testcase_26 AC 946 ms
4,900 KB
testcase_27 AC 969 ms
4,900 KB
testcase_28 AC 968 ms
4,904 KB
testcase_29 AC 947 ms
4,904 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <cmath>
#include <ctime>
#include <vector>
#include <cassert>
#include <iostream>
#include <algorithm>
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<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 + 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<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();
		}
		vector<double> 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<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 >= 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<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 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<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;
}
0