概要

平面上に1本の直線を引き、その直線との距離が \(D\) 以下となるセンサー（点）の数を最大化する問題です。

考察

直線を引いたとき、その直線との距離が \(D\) 以下である点は、直線に対して垂直な方向（法線方向）に射影したとき、直線の位置を中心とした幅 \(2D\) の区間に収まる点に対応します。 したがって、この問題は「ある方向に点を射影したとき、長さ \(2D\) 以下の区間に含まれる点の最大数を求める」問題に帰着できます。

射影する方向（法線ベクトル）は無限にありますが、すべてを調べる必要はありません。 最適な方向から少しずつ傾きを変えていくことを想像すると、点群が幅 \(2D\) の区間から外れようとする瞬間、あるいは新たに区間に入ろうとする瞬間が必ず訪れます。その瞬間は、以下のいずれかの条件を満たす方向となります： 1. 2つの点 \(i, j\) を結ぶ直線に垂直な方向（このとき2点の射影値が一致します） 2. 2つの点 \(i, j\) の射影値の差がちょうど \(2D\) になる方向

したがって、この2つの条件から生成される法線ベクトルを候補として全探索すれば、必ず最適解を見つけることができます。

アルゴリズム

1. 法線ベクトルの候補生成: すべての点のペア \((i, j)\) について、ベクトル \((dx, dy) = (X_i - X_j, Y_i - Y_j)\) を考えます。
   • 条件1: ベクトル \((dx, dy)\) に垂直な方向。単位ベクトルにすると \((a, b) = \left(\frac{dy}{L}, \frac{-dx}{L}\right)\) となります（\(L = \sqrt{dx^2 + dy^2}\)）。
   • 条件2: 射影値の差が \(2D\) になる方向。\(a \cdot dx + b \cdot dy = 2D\) かつ \(a^2 + b^2 = 1\) を満たす \((a, b)\) を解くと、\(L \ge 2D\) の場合に解を持ちます。解は三角関数を用いて計算でき、コードでは誤差を減らすため二乗和 \(L^2\) で割る形で計算しています。
2. 射影と尺取法: 各候補ベクトル \((a, b)\) について、すべての点の射影値 \(p_k = aX_k + bY_k\) を計算し、ソートします。 ソートされた配列に対して尺取法を行い、最大値と最小値の差が \(2D\) 以下となる区間の最大要素数を求めます。全候補の最大要素数が答えとなります。

計算量

• 時間計算量: \(O(N^3 \log N)\) 法線ベクトルの候補数は \(O(N^2)\) 個、各候補ごとに射影値のソートが \(O(N \log N)\) で行われます。\(N \le 200\) なので十分に高速です。
• 空間計算量: \(O(N^2)\) 候補ベクトルを保存するための領域が必要です。

実装のポイント

• 浮動小数点数の誤差対策: 尺取法で射影値の差を比較する際、proj[r] - proj[l] <= 2 * D + eps のように微小な定数 eps（例: \(10^{-9}\)）を足すことで、計算誤差による区間の見落としを防ぎます。

• 条件2の境界: 2点間の距離 \(L\) が \(2D\) よりも真に短い場合、射影値の差を \(2D\) にすることはできないため、この計算はスキップする必要があります。コード内では L2 >= V_2D * V_2D で判定しています。

  ソースコード

import sys
import math

def solve():
    input = sys.stdin.readline
    N, D = map(int, input().split())
    points = [tuple(map(int, input().split())) for _ in range(N)]

    candidates = [(1.0, 0.0)]
    eps = 1e-9
    V_2D = 2.0 * D

    for i in range(N):
        xi, yi = points[i]
        for j in range(i + 1, N):
            xj, yj = points[j]
            dx = xi - xj
            dy = yi - yj
            L2 = dx * dx + dy * dy
            
            if L2 > 0:
                L = math.sqrt(L2)
                a = dy / L
                b = -dx / L
                candidates.append((a, b))
                
                if V_2D > 0 and L2 >= V_2D * V_2D:
                    norm = math.sqrt(L2 - V_2D * V_2D)
                    a1 = (V_2D * dx + norm * dy) / L2
                    b1 = (V_2D * dy - norm * dx) / L2
                    candidates.append((a1, b1))

    ans = 0
    for a, b in candidates:
        proj = [a * x + b * y for x, y in points]
        proj.sort()
        
        r = 0
        for l in range(N):
            while r < N and proj[r] - proj[l] <= V_2D + eps:
                r += 1
            if r - l > ans:
                ans = r - l
            
    print(ans)

solve()

この解説は or-glm5.2-high によって生成されました。