結果

問題 No.5007 Steiner Space Travel
ユーザー BinomialSheepBinomialSheep
提出日時 2022-08-04 01:01:50
言語 C++17
(gcc 12.3.0 + boost 1.83.0)
結果
AC  
実行時間 983 ms / 1,000 ms
コード長 10,957 bytes
コンパイル時間 3,232 ms
実行使用メモリ 4,324 KB
スコア 8,879,279
最終ジャッジ日時 2022-08-04 01:02:26
合計ジャッジ時間 35,347 ms
ジャッジサーバーID
（参考情報） 		judge14 / judge12
純コード判定しない問題か言語
このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示
入力 結果 実行時間
実行使用メモリ
testcase_00 AC 981 ms
4,168 KB
testcase_01 AC 981 ms
4,180 KB
testcase_02 AC 981 ms
4,216 KB
testcase_03 AC 981 ms
4,116 KB
testcase_04 AC 981 ms
4,172 KB
testcase_05 AC 982 ms
4,180 KB
testcase_06 AC 982 ms
4,084 KB
testcase_07 AC 981 ms
4,172 KB
testcase_08 AC 981 ms
4,176 KB
testcase_09 AC 980 ms
4,184 KB
testcase_10 AC 981 ms
4,064 KB
testcase_11 AC 981 ms
4,084 KB
testcase_12 AC 982 ms
4,104 KB
testcase_13 AC 980 ms
4,216 KB
testcase_14 AC 981 ms
4,192 KB
testcase_15 AC 982 ms
4,272 KB
testcase_16 AC 981 ms
4,324 KB
testcase_17 AC 981 ms
4,028 KB
testcase_18 AC 982 ms
4,172 KB
testcase_19 AC 981 ms
4,032 KB
testcase_20 AC 981 ms
3,996 KB
testcase_21 AC 982 ms
4,116 KB
testcase_22 AC 981 ms
4,184 KB
testcase_23 AC 981 ms
4,320 KB
testcase_24 AC 980 ms
3,984 KB
testcase_25 AC 981 ms
4,180 KB
testcase_26 AC 981 ms
4,064 KB
testcase_27 AC 983 ms
4,176 KB
testcase_28 AC 981 ms
4,128 KB
testcase_29 AC 981 ms
3,988 KB
権限があれば一括ダウンロードができます

ソースコード

diff # 

#pragma GCC target("avx2")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")

#include <bits/stdc++.h>

using namespace std;
using namespace chrono;
using ll = long long;
using vi = vector<int>;
using vl = vector<long long>;
using vs = vector<string>;
using vc = vector<char>;
using vb = vector<bool>;
using vpii = vector<pair<int, int>>;
using vvi = vector<vector<int>>;
using vvl = vector<vector<long long>>;
using pii = pair<int, int>;
#define rep(i, n) for (int i = 0; i < (int)(n); i++)
#define all(x) (x).begin(), (x).end()
#define INFTY (1 << 30)

template <typename T>
inline bool chmax(T &a, T b) {
  return ((a < b) ? (a = b, true) : (false));
}
template <typename T>
inline bool chmin(T &a, T b) {
  return ((a > b) ? (a = b, true) : (false));
}

// 焼きなまし法の参考にしたページ
// https://shindannin.hatenadiary.com/entry/2021/03/06/115415

// 0以上UINT_MAX(0xffffffff)以下の整数をとる乱数 xorshift
// https://ja.wikipedia.org/wiki/Xorshift
static uint32_t randXor() {
  static uint32_t x = 123456789;
  static uint32_t y = 362436069;
  static uint32_t z = 521288629;
  static uint32_t w = 88675123;
  uint32_t t;

  t = x ^ (x << 11);
  x = y;
  y = z;
  z = w;
  return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
}

// 0以上1未満の小数をとる乱数
static double rand01() { return (randXor() + 0.5) * (1.0 / UINT_MAX); }

int N, M;
vi a, b;
vpii cd;
// α
const int A = 5;
const int MAP_SIZE = 1000;

// ワーシャルフロイドで求めた各惑星間の最短コスト
vvl G;

// 解説の「貪欲による解法2」を採用
// https://yukicoder.me/problems/no/5007/editorial
vi terryInit() {
  const ll INF = 1e18;
  G.resize(N);
  rep(i, N) rep(j, N) G[i].push_back(INF);

  auto calDist = [](int x1, int y1, int x2, int y2) {
    return (ll)((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
  };

  rep(i, N) rep(j, i) {
    G[i][j] = G[j][i] = A * A * calDist(a[i], b[i], a[j], b[j]);
  }
  vvl dirctG = G;

  rep(i, N) G[i][i] = 0;
  rep(k, N) rep(i, N) {
    if (G[i][k] == INF) continue;
    rep(j, N) {
      if (G[k][j] == INF) continue;
      chmin(G[i][j], G[i][k] + G[k][j]);
    }
  }
  vi ret;

  // 経路復元できるダイクストラ
  auto dijkstra = [&](int now, int nexV) {
    vl dist(N, INFTY);
    vi prev(N, -1);
    priority_queue<pii, vpii, greater<pii>> pq;
    pq.emplace(0, now);
    dist[now] = 0;
    while (!pq.empty()) {
      auto p = pq.top();
      pq.pop();
      int from = p.second;
      if (dist[from] < p.first) continue;
      rep(to, N) {
        if (chmin(dist[to], dist[from] + dirctG[from][to])) {
          prev[to] = from;
          pq.emplace(dist[to], to);
        }
      }
    }

    // 経路復元する
    int current = nexV;
    stack<int> nowToNexRev;
    while (current != now) {
      nowToNexRev.push(current + 1);
      current = prev[current];
    }
    while (!nowToNexRev.empty()) {
      ret.push_back(nowToNexRev.top());
      nowToNexRev.pop();
    }
  };

  vi color(N);
  int now = 0;
  color[0] = 1;
  ret.push_back(1);

  rep(i, N - 1) {
    // 未訪問で1番近い頂点を見つける
    ll nexD = INF;
    int nexV = -1;
    rep(j, N) {
      if (color[j]) continue;
      if (chmin(nexD, G[now][j])) nexV = j;
    }
    // ダイクストラして経路復元
    dijkstra(now, nexV);
    //
    now = nexV;
    color[nexV] = 1;
  }
  // 最終地点から1までもダイクストラ
  dijkstra(now, 0);

  return ret;
}

// 宇宙ステーションの初期配置
vpii initStation() {
  vpii ret;
  ret.emplace_back(300, 300);
  ret.emplace_back(300, 500);
  ret.emplace_back(300, 700);
  ret.emplace_back(500, 300);
  ret.emplace_back(500, 700);
  ret.emplace_back(700, 300);
  ret.emplace_back(700, 500);
  ret.emplace_back(700, 700);
  return ret;
}

ll calcDist(int prev, int next) {
  --prev, --next;
  if (prev < N && next < N) {
    return G[prev][next];
  }
  auto dist = [](int x1, int y1, int x2, int y2) {
    return (ll)((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
  };
  int x1 = prev < N ? a[prev] : cd[prev - N].first;
  int y1 = prev < N ? b[prev] : cd[prev - N].second;
  int x2 = next < N ? a[next] : cd[next - N].first;
  int y2 = next < N ? b[next] : cd[next - N].second;
  ll ret = dist(x1, y1, x2, y2);
  if (prev < N || next < N) ret *= 5;
  return ret;
}

const static double START_TEMP = 1e4;  // 開始時の温度
const static double END_TEMP = 10;     // 終了時の温度
const double MAX_DELTA = 400;
const double MIN_DELTA = 10;

ll annealing(double endTime, vi &tr, ll currentScore) {
  /* 焼きなまし法 */
  auto startClock = system_clock::now();

  double time = 0.0;  // 経過時間（秒）
  double progressRatio = 0;
  double temp = START_TEMP;

  // ループ回数
  int cnt = 0;
  do {
    cnt++;
    // if (cnt == 30) break;
    //    進捗。開始時が0.0で、終了時が1.0
    if (cnt % 2000) {
      // 時間更新
      time =
          ((double)duration_cast<microseconds>(system_clock::now() - startClock)
               .count() *
           1e-6);
      // 温度更新
      progressRatio = time / endTime;
      temp = START_TEMP + (END_TEMP - START_TEMP) * progressRatio;
    }

    // 近傍の割合
    int neighType = randXor() % 10;

    if (neighType < 2) {
      /* 近傍1：ステーションを適当な位置に挿入 */
      int stationId = N + randXor() % M + 1;
      // 最初と最後は1固定なので最初や最後には入れない
      int idx = 1 + randXor() % ((int)tr.size() - 1);
      int prev = tr[idx - 1];
      int next = tr[idx];
      ll oldScore = calcDist(prev, next);
      ll newScore = calcDist(prev, stationId) + calcDist(stationId, next);
      // 解の更新
      ll deltaScore = oldScore - newScore;
      const double probability = exp((double)deltaScore / temp);
      if (probability >= rand01()) {
        currentScore -= deltaScore;
        // 挿入する
        tr.insert(tr.begin() + idx, stationId);
      }
    } else if (neighType < 4) {
      /* 近傍2：適当な位置のステーションを削除 */
      int len = (int)tr.size();
      // ステーションを探す（適当に上限を10としている）
      int idx = 0;
      int trial = 0;
      while (tr[idx] <= N && trial < 10) {
        idx = randXor() % len;
        trial++;
      }
      if (tr[idx] <= N) continue;
      // ステーションを指すidxが見つかった場合
      int stationId = tr[idx];
      int prev = tr[idx - 1];
      int next = tr[idx + 1];
      ll oldScore = calcDist(prev, stationId) + calcDist(stationId, next);
      ll newScore = calcDist(prev, next);
      // 解の更新
      ll deltaScore = oldScore - newScore;
      const double probability = exp((double)deltaScore / temp);
      if (probability >= rand01()) {
        currentScore -= deltaScore;
        // 挿入する
        tr.erase(tr.begin() + idx);
      }
    } else if (neighType < 5) {
      /* 近傍3：ステーションの移動 */
      int stationId = N + randXor() % M + 1;
      // 暫定解からステーションを削除
      vi tmpTr;
      rep(i, (int)tr.size()) if (tr[i] != stationId) tmpTr.push_back(tr[i]);
      // ステーションを移動（移動範囲も徐々に小さくしてる！）
      pii oldP = cd[stationId - N - 1];
      pii newP = oldP;
      int d = (int)(MAX_DELTA * (endTime - time) + MIN_DELTA * time);
      int mn = max(0, newP.first - d);
      newP.first = mn + randXor() % (min(MAP_SIZE, newP.first + d) - mn);
      mn = max(0, newP.second - d);
      newP.second = mn + randXor() % (min(MAP_SIZE, newP.second + d) - mn);
      cd[stationId - N - 1] = newP;
      // 各辺でステーションを使うかgreedyに判定
      vi newTr;
      ll newScore = 0;
      rep(i, (int)tmpTr.size() - 1) {
        newTr.push_back(tmpTr[i]);
        ll oldDist = calcDist(tmpTr[i], tmpTr[i + 1]);
        ll newDist =
            calcDist(tmpTr[i], stationId) + calcDist(stationId, tmpTr[i + 1]);
        if (newDist < oldDist) {
          newScore += newDist;
          newTr.push_back(stationId);
        } else {
          newScore += oldDist;
        }
      }
      newTr.push_back(tmpTr.back());
      // 解の更新
      ll deltaScore = currentScore - newScore;
      const double probability = exp((double)deltaScore / temp);
      if (probability >= rand01()) {
        currentScore = newScore;
        tr.resize(newTr.size());
        tr = newTr;
      } else {
        cd[stationId - N - 1] = oldP;
      }
    } else {
      /* 近傍4：2-opt */
      // 最初と最後は1固定なので入れ替えない
      int l = (int)tr.size();
      int from = 1 + randXor() % (l - 2);
      int to = 1 + randXor() % (l - 2);
      if (from == to) {
        to += 1;
        if (to == l - 1) to = 2;
      }
      if (from > to) swap(from, to);
      int i0 = tr[from - 1];
      int i1 = tr[from];
      int i2 = tr[to - 1];
      int i3 = tr[to];
      ll d0123 = calcDist(i0, i1) + calcDist(i2, i3);
      ll d0213 = calcDist(i0, i2) + calcDist(i3, i1);
      ll newScore = currentScore + d0213 - d0123;
      // 解の更新
      ll deltaScore = currentScore - newScore;
      const double probability = exp((double)deltaScore / temp);
      if (probability >= rand01()) {
        currentScore = newScore;
        // 繋ぎ直す
        reverse(tr.begin() + from, tr.begin() + to);
      }
    }

    // trを綺麗にする
    // if (cnt % 1000) {
    //   vi cleanTr;
    //   cleanTr.push_back(tr[0]);
    //   rep(i, (int)tr.size() - 1) {
    //     if (tr[i + 1] != cleanTr.back()) cleanTr.push_back(tr[i + 1]);
    //   }
    //   tr.resize((int)cleanTr.size());
    //   tr = cleanTr;
    // }

  } while (time < endTime);

  // 戻り値
  return currentScore;
}

int main() {
  ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  /* input */
  cin >> N >> M;
  a.resize(N);
  b.resize(N);
  rep(i, N) cin >> a[i] >> b[i];

  /* solve */
  // 解説の「貪欲法による解法（その2）」を再現
  // ワーフロして経路復元ダイクストラ
  vi initTr = terryInit();

  // debug(calcScore(tr));
  // auto maxScore = calcScore(tr);
  ll currentScore = 0;
  rep(i, (int)initTr.size() - 1) currentScore +=
      calcDist(initTr[i], initTr[i + 1]);

  ll bestScore = currentScore;
  vi bestTr;
  vpii bestCd(8);
  // 3回焼きなます
  cd = initStation();
  rep(i, 3) {
    vi thisTr = initTr;
    ll thisScore = annealing(0.322, thisTr, currentScore);
    if (chmin(bestScore, thisScore)) {
      bestTr = thisTr;
      bestCd = cd;
    }
  }

  /* output */
  int V = (int)bestTr.size();
  // assert(bestTr[0] == 1);
  // assert(bestTr[V - 1] == 1);

  rep(i, M) cout << bestCd[i].first << " " << bestCd[i].second << "\n";

  cout << V << endl;
  rep(i, V) {
    if (bestTr[i] <= N) {
      cout << "1 " << bestTr[i] << "\n";
    } else {
      cout << "2 " << bestTr[i] - N << "\n";
    }
  }
  return 0;
}
0