点 \(A,B\) について \(|AB|\) を線分 \(AB\) の長さと定義します。

実は、\(\displaystyle \sum_{i=1}^{N}|P_i Q_{R_i}|\) の最小値を達成するような \(R\) がこの問題の解になります。

証明

まず、以下の命題が成り立つことを示します。

問題の条件を満たさない任意の \(R\) について、\(\displaystyle \sum_{i=1}^{N}|P_i Q_{R'_i}| < \displaystyle \sum_{i=1}^{N}|P_i Q_{R_i}|\) を満たすような \(R'\) が存在する。

問題の条件を満たさない \(R\) において、線分 \(P_a Q_{R_a}\) と線分 \(P_b Q_{R_b}\) \((1 \leq a < b \leq N)\)が交わっているとします。このとき、\(R'_k= \left \{ \begin{array}{ll} R_k &(k \neq a \land k \neq b) \\ R_a &(k = b) \\ R_b &(k=a) \end{array} \right\} \) となる \(R'\) は \(\displaystyle \sum_{i=1}^{N}|P_i Q_{R'_i}| < \displaystyle \sum_{i=1}^{N}|P_i Q_{R_i}|\) を満たします。 これは線分 \(P_a Q_{R_a}\) と線分 \(P_b Q_{R_b}\) の交点を \(X\) としたとき、三角形の成立条件より \(|P_aX|+|Q_{R_b}X|>|P_aQ_{R_b}|\)　及び \(|P_bX|+|Q_{R_a}X|>|P_bQ_{R_a}|\) が言え、\(|P_aQ_{R_a}|+|P_bQ_{R_b}|>|P_aQ_{R_b}|+|P_bQ_{R_a}|\) と示せるため成り立ちます。

また、 \(\displaystyle \sum_{i=1}^{N}|P_i Q_{R_i}|\) は有界であるため最小値が存在します。 よって、先ほどの命題より \(\displaystyle \sum_{i=1}^{N}|P_i Q_{R_i}|\) の最小値を満たす \(R\) は問題文の条件を満たすことが示せます。 (証明終わり)

実際に最小値を達成するような \(R\) を求めるためには最小重み完全マッチング問題を解く必要があります。 最小重み完全マッチング問題はハンガリアン法や最小費用流を用いて高速に求めることができます。最小費用流問題について初めて知ったという方はdrkenさんの記事 などを参照してください。また、今回の場合は同様の問題がLibrary Checkerにあります。(問題)

AtCodoer Library(ACL) には最小費用流ライブラリがあるため、それを用いることで簡潔に実装出来ます。ただし、ここで注意点があり、ACLの最小費用流では整数型のみを保証していて、コストの型を浮動小数点数にすることは厳密には未定義です。そのため、実際のコストをある定数 \(C\) との積を取って四捨五入した値を最小費用流で用いるコストとする必要があります。この \(C\) は今回の制約下では少なくとも \(1.413 \times 10^{12}\) 以上のときでないと正当性を示せません。座標の絶対値が \(M\) 以下であるとき、\(C\) のオーダーはすくなくとも\(O(M^3)\) である必要があります。

しかし、この問題ではACLの最小費用流でコストの型をdouble型やlong double型としてもACが取れるようになっています。

#include "atcoder/mincostflow"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
static constexpr long double C = 1.5e12; //4/(sqrt(4999^2+5000^2)+sqrt(4997^2+4998^2)-2*sqrt(4998^2+4999^2))=1.413*10^12
int main() {
    int n;
    cin >> n;
    vector<ll> xb(n), yb(n), xw(n), yw(n);
    for (int i = 0; i < n; ++i) {
        cin >> xb[i] >> yb[i];
    }
    for (int i = 0; i < n; ++i) {
        cin >> xw[i] >> yw[i];
    }
    auto dist = [&](int i, int j) -> long double {
        return sqrtl((xb[i] - xw[j]) * (xb[i] - xw[j]) + (yb[i] - yw[j]) * (yb[i] - yw[j]));
    };
    atcoder::mcf_graph<ll, ll> g(2 * n + 2);
    int s = 2 * n, t = 2 * n + 1;
    for (int i = 0; i < n; ++i) {
        g.add_edge(s, i, 1, 0);
        g.add_edge(i + n, t, 1, 0);
        for (int j = 0; j < n; ++j) {
            g.add_edge(i, j + n, 1, (ll)(dist(i, j) * C + 0.5));
        }
    }
    auto [flow, cost] = g.flow(s, t);
    int edge_id = 0;
    vector<int> ans(n);
    for (int i = 0; i < n; ++i) {
        edge_id += 2;
        for (int j = 0; j < n; ++j) {
            auto e = g.get_edge(edge_id);
            if (e.flow == 1)
                ans[i] = j + 1;
            edge_id++;
        }
    }
    for (int i = 0; i < n; ++i) {
        cout << ans[i] << " \n"[i == n - 1];
    }
}