Problem 816: Shortest Distance Among Points

Mathematical Analysis

This is the classic Closest Pair of Points problem from computational geometry.

Brute Force Approach

Checking all $\binom{n}{2}$ pairs has $O(n^2)$ complexity, which is infeasible for $n = 2{,}000{,}000$.

Divide and Conquer Algorithm

The classic algorithm by Shamos and Hoey achieves $O(n \log n)$:

1. Sort points by x-coordinate.
2. Divide the set into two halves by a vertical line.
3. Conquer recursively: find the closest pair in each half, let $\delta = \min(d_L, d_R)$.
4. Combine: Check pairs that straddle the dividing line. Only points within distance $\delta$ of the dividing line need checking, and for each such point, at most 7 other points in the strip need comparison.

Grid/Hash-Based Approach (Used Here)

An alternative $O(n)$ expected-time randomized algorithm:

1. Compute an initial estimate $\delta$ from the first few points or a random sample.
2. Build a grid with cell size $\delta$.
3. For each new point, check only neighboring cells (at most 9 cells in 2D).
4. If a closer pair is found, rebuild the grid with the new $\delta$.

For practical purposes with pseudo-random points, a sorted strip-based approach works well:

1. Sort points by x-coordinate.
2. Use a set (balanced BST) sorted by y-coordinate for the “active” window.
3. For each point, remove points whose x-distance exceeds the current best $\delta$.
4. Query nearby points in the y-sorted set.

This sweep-line approach runs in $O(n \log n)$.

Editorial

We generate all 2,000,000 points using the recurrence. We then sort by x-coordinate. Finally, sweep from left to right, maintaining a set of active points sorted by y.

Pseudocode

Generate all 2,000,000 points using the recurrence
Sort by x-coordinate
Sweep from left to right, maintaining a set of active points sorted by y
For each point p:
Return delta

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

• Time: $O(n \log n)$ for sorting and sweep line operations.
• Space: $O(n)$ for storing points and the active set.

Answer

#include <bits/stdc++.h>
using namespace std;

int main(){
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const int K = 2000000;
    const long long MOD = 50515093LL;

    // Generate sequence s_0, s_1, s_2, ...
    // s_0 = 290797, s_{n+1} = s_n^2 mod 50515093
    // P_n = (s_{2n}, s_{2n+1})
    vector<pair<long long,long long>> pts(K);
    long long s = 290797;
    for(int i = 0; i < K; i++){
        long long x = s;
        s = (s * s) % MOD;
        long long y = s;
        s = (s * s) % MOD;
        pts[i] = {x, y};
    }

    // Sort by x
    sort(pts.begin(), pts.end());

    // Sweep line with set sorted by (y, x)
    auto cmp = [](const pair<long long,long long>& a, const pair<long long,long long>& b){
        if(a.second != b.second) return a.second < b.second;
        return a.first < b.first;
    };
    set<pair<long long,long long>, decltype(cmp)> active(cmp);

    double best = 1e18;
    int left = 0;

    for(int i = 0; i < K; i++){
        long long d = (long long)ceil(best);

        // Remove points too far left
        while(left < i && pts[i].first - pts[left].first > d){
            active.erase(pts[left]);
            left++;
        }

        // Check nearby points in y range
        pair<long long,long long> lo = {-1, pts[i].second - d};
        auto it = active.lower_bound(lo);

        for(; it != active.end() && it->second <= pts[i].second + d; ++it){
            double dx = (double)(pts[i].first - it->first);
            double dy = (double)(pts[i].second - it->second);
            double dist = sqrt(dx*dx + dy*dy);
            if(dist < best){
                best = dist;
            }
        }

        active.insert(pts[i]);
    }

    cout << fixed << setprecision(9) << best << endl;

    return 0;
}

import math
from sortedcontainers import SortedList

def solve():
    K = 2000000
    MOD = 50515093

    # Generate points: s_0 = 290797, s_{n+1} = s_n^2 mod M
    # P_n = (s_{2n}, s_{2n+1})
    pts = []
    s = 290797
    for i in range(K):
        x = s
        s = (s * s) % MOD
        y = s
        s = (s * s) % MOD
        pts.append((x, y))

    # Sort by x-coordinate
    pts.sort()

    # Sweep line algorithm
    # Use SortedList keyed by (y, x)
    active = SortedList(key=lambda p: (p[1], p[0]))
    best = float('inf')
    left = 0

    for i in range(K):
        d = best

        # Remove points too far left
        while left < i and pts[i][0] - pts[left][0] > d:
            active.remove(pts[left])
            left += 1

        # Find points in y-range [pts[i][1] - d, pts[i][1] + d]
        lo_y = pts[i][1] - d
        hi_y = pts[i][1] + d

        # Use irange or bisect to find candidates
        idx_lo = active.bisect_key_left((lo_y, -1))
        idx_hi = active.bisect_key_right((hi_y, float('inf')))

        for j in range(idx_lo, min(idx_hi, len(active))):
            p = active[j]
            if p[1] > hi_y:
                break
            dx = pts[i][0] - p[0]
            dy = pts[i][1] - p[1]
            dist = math.sqrt(dx*dx + dy*dy)
            if dist < best:
                best = dist

        active.add(pts[i])

    print(f"{best:.9f}")

if __name__ == "__main__":
    solve()