この問題は様々な解法があります。今回は Mo’s algorithm 内部で利用される手法を適用する解法を説明します。

想定解は簡単に言うと、領域内を「蛇行」するように移動する手法です。
まず、この問題は ハミルトンサイクル 、すなわち全ての地点を \(1\) 回ずつ通る閉路を求める問題です。ハミルトンサイクルが 1 個求まれば、頂点順を rotate することで頂点 \(1\) を始点とする閉路に変換できるので、始点がどこであるかはいったん考えなくてもよいことになります。つまり、いったん始点が地点 \(1\) である必要があることは忘れて、閉路を探す問題として考えましょう。
座標の範囲を \(W = 2 \times 10^7\) で表します。まず、幅 \(B\) を決めて、領域を \(x\) 軸側が幅 \(B\) になるように短冊状に分割します。短冊状の領域はおよそ \(W/B\) 個生成されます。その後、左の領域から順に地点を訪問していきます。ここで短冊状の領域内では、

• 左から奇数番目の領域では \(y\) 座標 昇順 に、
• 左から偶数番目の領域では \(y\) 座標 降順 になるように訪問します。

以下の図はイメージ図です。(右方向が \(x\) 軸方向、上方向が \(y\) 軸方向) 図には書いていませんが、求めるのはハミルトンサイクルなので最後に終着点から最初の点に戻ることを忘れないようにしましょう。

この解法でかかる時間の最悪ケースを考えてみましょう。大きめに見積もると次のようになります：

• 短冊内の縦方向の移動 : 1 個の短冊で最大 \(W\) 秒移動して、短冊が \(W/B\) 個だから \(W^2/B\) 秒
• 短冊内の横方向の移動 : 地点間の移動が最大 \(N\) 回発生して、1 回の移動で最大横方向に \(B\) 動くので \(NB\) 秒
• 隣り合う短冊をまたぐ移動：\(W/B\) 回程度発生して毎回最大で \(2B\) 動くので、\(2W\) 秒
• 最後の終点から始点へ戻る移動：最大 \(2W\) 動くので、\(2W\) 秒
• 合計 \(W^2 / B + NB + 4W\) 秒

よって相加相乗平均の公式を考えると \(B = \frac{W}{\sqrt{N}}\) 程度の値を取るのが最適とわかり、この時合計で \(2 W \sqrt{N} + 4W\) 秒かかることになります。
今回の \(N=6\times 10^4, W=2\times 10^7\) という制約では \(2 W \sqrt{N} + 4W \simeq 9.87 \times 10^9\) なので \(10^{10}\) 秒に収まります。計算量はソートにかかる \(\mathrm{O}(N \log N)\) で十分高速です。

別解について解説します。今回の問題の想定解は距離の総和が \(2W\sqrt{N}\) となっていますが、定数倍を無視すると、距離の総和が \(\mathrm{O}(W \sqrt{N})\) となる解法は無数にあります。列挙すると次の通りです：

• 想定解とほぼ同じだが細部で定数倍を悪化させている
  • \(B\) の値を \(\frac{W}{\sqrt{N}}\) から大きく外れた値を取る
  • 蛇行せずに毎回 \(y\) 昇順に頂点を訪問する
• Hilbert curve や Z-order curve の順に頂点を訪問する
• マンハッタン距離上の最小全域木を求めて(これは \(\mathrm{O}(N \log N)\) で求まる) DFS 順に頂点を訪問する
• \(\vdots\)

いずれの解法も定数倍が想定解よりわずかに悪いため、実装によっては \(10^{10}\) 秒を超えてしまうようです。しかしその場合でもヒューリスティックな解法を適用するとこの問題を通すことが出来ると考えられます。いくつか例を挙げると次の通りです：

• いくつかの手法を実装して距離最小のものを選ぶ
• \(10^{10}\) 秒に収まるまで \(B\) の値を乱択する
• 2-opt を用いた山登り法を時間いっぱい行う

最後に、少し面白い事実を紹介します。今回紹介した解法は最悪ケースでの距離の総和が \(2W\sqrt{N} + \mathrm{O}(W)\) でしたが、 noshi91 さんにより次の事実が示されました。(ブログ)

• ad-hoc な工夫をすると最悪ケースでの距離の総和が \(\sqrt{2} W \sqrt{N} + \mathrm{O}(W)\) となる手法を構成できる
• また、どのようなハミルトンサイクルを選んでも必ず \(\sqrt{2} W \sqrt{N}\) かかるケースが存在する
• よって提案手法は \(\mathrm{O}(W)\) の差を除いて最適である

つまり、今回の問題で制約を「 \(7 \times 10^9\) 秒以下」にしてもキレイな構築が存在して、かつその構築は理論的な下界を達成していることになります。内容も面白いので是非ブログの方を読んでみてください。(ブログに書かれている内容は三角形領域で議論しているので定数倍が \(\sqrt{2}\) 倍異なりますが、今回の問題にもそのまま適用できます。)

• 実装例(C++)

#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <vector>
using namespace std;

int main() {
  cin.tie(0)->sync_with_stdio(0);
  int N;
  cin >> N;
  vector<int> X(N), Y(N);
  for (int i = 0; i < N; i++) cin >> X[i] >> Y[i];

  int B = 2e7 / sqrt(N);
  vector<int> ord(N);
  iota(begin(ord), end(ord), 0);
  sort(begin(ord), end(ord), [&](int i, int j) {
    if (X[i] / B != X[j] / B) return X[i] / B < X[j] / B;
    if ((X[i] / B) % 2) return Y[i] > Y[j];
    return Y[i] < Y[j];
  });
  int pos0 = 0;
  while (ord[pos0] != 0) pos0++;
  rotate(begin(ord), begin(ord) + pos0, end(ord));

  for (int i = 0; i < N; i++) cout << ord[i] + 1 << " \n"[i + 1 == N];
}