This problem has various solutions. In this editorial, we explain a solution that uses the working principle of Mo’s algorithm.

Simply put, in the writer’s solution, the tour looks like a “snake” shape.
First of all, this problem asks to find a Hamiltonian cycle, which is a cycle that visits every location exactly once. Once a Hamiltonian cycle is found, the vertex order can be rotated to obtain a cycle that starts at vertex \(1\), so for now we may ignore the starting point. Thus, let us ignore the requirement of the starting point being location \(1\), and interpret it as a problem to find a cycle.
Denote the coordinate range by \(W = 2 \times 10^7\). First, fix a width \(B\) and split the region into vertical rectangular region of width \(B\). There will be about \(W/B\) sub-rectangles. Next, we visit the locations within each sub-rectangle, in order. Within a sub-rectangle, the locations are visited in

• ascending order of \(y\) within the odd-indexed sub-rectangles, and
• descending order of \(y\) within the even-indexed sub-rectangles.

Below is a conceptual figure. (\(x\) axis points right, and \(y\) axis points up.) Although not illustrated in the figure, do not forget that you need to come back from the final point to the first point.

Let us consider the worst-case required time of a tour. It can be roughly estimated as follows:

• Vertical move within a sub-rectangle: at most \(W\) seconds within a sub-rectangle. With \(W/B\) sub-rectangle, the total is \(W^2/B\) seconds.
• Horizontal move within a sub-rectangle: at most \(N\) moves between locations can occur, with a horizontal move by at most \(B\) per move, for a total of \(NB\) seconds.
• Move between adjacent sub-rectangles: occurs about \(W/B\) times, with each costing \(2B\), for a total of \(2W\) seconds.
• Final move from the last to first location: moves at most \(2W\), so it costs \(2W\) seconds.
• Total: \(W^2 / B + NB + 4W\)

By AM-GM (arithmetic mean-geometric mean) inequality, it is optimal to take a value of about \(B = \frac{W}{\sqrt{N}}\), in which case the cost is about \(2 W \sqrt{N} + 4W\).
Under the constraints of the problem (\(N=6\times 10^4, W=2\times 10^7\)), we have \(2 W \sqrt{N} + 4W \simeq 9.87 \times 10^9\), which fits within\(10^{10}\) seconds. The complexity is \(\mathrm{O}(N \log N)\) originating from the sorting, which is fast enough.

We will introduce alternative solutions. In the solution above, the distances sum to \(2W\sqrt{N}\), but if we ignore the constant factor, there are numerous solutions whose distance sums to \(\mathrm{O}(W \sqrt{N})\). They include:

• A solution that is basically the same as above, but the details are worsening the constant factor:
  • \(B\) is far off from \(\frac{W}{\sqrt{N}}\).
  • Instead of snake shape, each block is always scanned in ascending order of \(y\).
• Visit in the order of Hilbert curve or Z-order curve.
• \(\vdots\)

All these solutions have a slightly worse constant factor than the intended solution, so some implementation may lead to exceeding the \(10^{10}\)-second limit. Nevertheless, we may apply heuristic improvements to squeeze it into the limit. For example:

• Implement several solutions and choose one with the minimum traveling distance
• Randomly pick \(B\) until the total fits in \(10^{10}\) seconds
• Keep running a hill-climbing algorithm based on 2-opt within the execution time limit

Finally, we will introduce a fun fact. In the solution explained here, the total distance at the worst case was about \(2W\sqrt{N} + \mathrm{O}(W)\). In fact, noshi91 proved (in the blog (Japanese)) that:

• With an ad-hoc improvement, we can construct an algorithm where the total distance at the worst case is \(\sqrt{2} W \sqrt{N} + \mathrm{O}(W)\).
• Conversely, there is a case where any Hamiltonian cycle has a total distance of \(\sqrt{2} W \sqrt{N}\).
• Hence, the proposed method is optimal up to an \(\mathrm{O}(W)\) difference.

This implies that setting the upper bound in the problem to \(7 \times 10^9\) still has a clean construction, and that construction achieves the theoretical lower bound. The contents are also interesting, so we recommend you to read the blog. (The blog discusses a triangular region, so the constant factor is different by \(\sqrt{2}\), but it is directly applicable to our problem.)

• Sample code (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];
}