22.5 คู่จุดที่ใกล้ที่สุด

สามหน้าที่ผ่านมาใช้ divide and conquer กับข้อมูลที่เรียงเป็นแถว หน้านี้เป็นตัวอย่างขั้นสูงบนระนาบสองมิติ มีจุด N จุด หาคู่ที่อยู่ใกล้กันที่สุด เช่น เสาสัญญาณที่ตั้งใกล้กันเกินจะรบกวนกัน ลองทุกคู่เป็น O(N²) ซึ่งที่ N เป็นแสนต้องคำนวณห้าพันล้านคู่ วิธีแบ่งครึ่งในหน้านี้ลดเหลือ O(N log² N) ความยากอยู่ที่ขั้นรวม ต้องหาคู่ที่ข้ามเส้นแบ่งให้ได้โดยไม่ลองทุกคู่ข้ามฝั่ง ซึ่งเป็นความคิดเดียวกับการนับคู่ข้ามครึ่งในหน้า 22.3 แต่ใช้เรขาคณิตช่วยตัดคู่ที่ไม่มีหวังทิ้ง

โจทย์และวิธีลองทุกคู่

เสาสัญญาณ N ต้นอยู่ที่พิกัดจำนวนเต็ม (x, y) ระยะระหว่างสองต้นคือ √((x₁ - x₂)² + (y₁ - y₂)²) ตามหน้า 12.5 เพื่อไม่ต้องยุ่งกับทศนิยม ให้เทียบกันด้วย ระยะกำลังสอง ซึ่งเป็นจำนวนเต็ม ระยะจริงยิ่งน้อย ระยะกำลังสองก็ยิ่งน้อย ลำดับจึงเหมือนกัน วิธีตรงไปตรงมาคือลูปซ้อนเทียบทุกคู่ i < j เป็น O(N²) เขียนได้ในไม่กี่บรรทัด และควรเขียนไว้เทียบคำตอบเสมอ

long long best = LLONG_MAX;
for (int i = 0; i < n; i++) {
    for (int j = i + 1; j < n; j++) {
        best = min(best, dist2(towers[i], towers[j]));
    }
}

แบ่งด้วยเส้นตั้ง

เรียงเสาตามพิกัด x แล้วแบ่งที่ต้นตรงกลาง ครึ่งซ้ายกับครึ่งขวามีจำนวนเท่าๆ กัน หาคู่ที่ใกล้ที่สุดของแต่ละครึ่งด้วยการเรียกตัวเอง ได้ระยะกำลังสองที่ดีที่สุดของสองฝั่ง ให้ d เป็นค่าที่น้อยกว่า เหลือแต่คู่ที่ข้ามเส้นแบ่ง ซึ่งการเรียกตัวเองมองไม่เห็น

คู่ข้ามเส้นที่จะชนะ d ได้ ต้องห่างกันตามแนว x น้อยกว่า √d ทั้งสองต้นจึงต้องอยู่ห่างจากเส้นแบ่งไม่ถึง √d ต้นที่อยู่ไกลกว่านั้นตัดทิ้งได้เลย เหลือเฉพาะต้นใน แถบ กว้าง √d สองข้างของเส้น ในโค้ดเทียบแบบกำลังสองทั้งหมด คือเก็บต้นที่ (x - x_กลาง)² < d

ตัวอย่าง เสา 8 ต้น เรียงตาม x แล้วครึ่งซ้ายคือ (2, 2) (5, 30) (8, 27) (21, 27) ครึ่งขวาคือ (22, 13) (23, 28) (28, 21) (34, 15) เส้นแบ่งอยู่ที่ x = 22 ซึ่งเป็นต้นแรกของครึ่งขวา

ครึ่งซ้ายได้ 18 จาก (5, 30) กับ (8, 27) ครึ่งขวาได้ 72 จาก (28, 21) กับ (34, 15) d จึงเป็น 18 แถบคือต้นที่ห่างจาก x = 22 ไม่ถึง √18 ราว 4.24 เหลือแค่ (21, 27) (22, 13) (23, 28) และคู่ (21, 27) กับ (23, 28) ห่างกันกำลังสองเท่ากับ 2² + 1² = 5 ชนะทั้งสองฝั่ง คำตอบคือ 5

ในแถบเทียบแค่ไม่กี่ต้น

แถบอาจยังมีเสาเยอะ ถ้าเสาเกือบทั้งหมดอยู่ใกล้เส้นแบ่ง ลองทุกคู่ในแถบก็กลับไปเป็น O(N²) เคล็ดคือเรียงต้นในแถบตามพิกัด y แล้วเทียบแต่ละต้นกับต้นถัดไปข้างบนเท่านั้น หยุดทันทีที่ผลต่างของ y ยกกำลังสองไม่น้อยกว่า d เพราะต้นที่สูงขึ้นไปอีกยิ่งห่างตามแนว y มากกว่านั้น

ทำไมถึงเทียบแค่ไม่กี่ต้น ต้นที่ต้องเทียบกับต้น p อยู่ในกรอบสี่เหลี่ยมกว้าง 2√d (แถบทั้งสองฝั่ง) สูง √d เหนือ p ลองตัดกรอบนี้เป็นช่องสี่เหลี่ยมจัตุรัสด้านละ √d / 2 ได้ 8 ช่อง ฝั่งละ 4 ช่อง แต่ละช่องมีเสาได้ไม่เกินหนึ่งต้น เพราะเสาสองต้นในช่องเดียวกันอยู่ฝั่งเดียวกันของเส้นแบ่ง และห่างกันไม่เกินเส้นทแยงของช่องคือ √d / √2 ซึ่งน้อยกว่า √d แต่ d คือค่าที่น้อยกว่าของคำตอบสองฝั่ง จึงไม่มากกว่าคำตอบของฝั่งไหนเลย ทุกคู่ที่อยู่ฝั่งเดียวกันจึงห่างกันอย่างน้อย √d ขัดกัน กรอบนี้จึงมีเสาไม่เกิน 8 ต้นรวม p เอง แต่ละต้นในแถบจึงเทียบกับต้นถัดไปไม่เกิน 7 ต้น

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

vector<pair<int, int>> towers;
long long checks = 0;

long long dist2(pair<int, int> p, pair<int, int> q)
{
    long long dx = p.first - q.first; // (1)!
    long long dy = p.second - q.second;
    return dx * dx + dy * dy;
}

long long closest(int lo, int hi)
{
    long long best = LLONG_MAX;
    if (hi - lo <= 3) { // (2)!
        for (int i = lo; i < hi; i++) {
            for (int j = i + 1; j < hi; j++) {
                checks++;
                best = min(best, dist2(towers[i], towers[j]));
            }
        }
        return best;
    }
    int mid = (lo + hi) / 2;
    int mid_x = towers[mid].first;
    best = min(closest(lo, mid), closest(mid, hi));
    vector<pair<int, int>> strip;
    for (int i = lo; i < hi; i++) {
        long long dx = towers[i].first - mid_x;
        if (dx * dx < best) {
            strip.push_back(towers[i]);
        }
    }
    sort(strip.begin(), strip.end(), [](pair<int, int> p, pair<int, int> q) { // (3)!
        return p.second < q.second;
    });
    for (int i = 0; i < (int)strip.size(); i++) {
        for (int j = i + 1; j < (int)strip.size(); j++) {
            long long dy = strip[j].second - strip[i].second;
            if (dy * dy >= best) { // (4)!
                break;
            }
            checks++;
            best = min(best, dist2(strip[i], strip[j]));
        }
    }
    return best;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n;
    cin >> n;
    towers.resize(n);
    for (int i = 0; i < n; i++) {
        cin >> towers[i].first >> towers[i].second;
    }
    sort(towers.begin(), towers.end()); // (5)!
    long long best = closest(0, n);
    cout << "closest squared distance = " << best << '\n';
    cout << "checks = " << checks << '\n';
    return 0;
}

หน้านี้พิกัดอยู่ระหว่าง 0 ถึง 10⁹ ผลต่างจึงยังอยู่ใน int แต่ยกกำลังสองแล้วถึง 10¹⁸ ต้องเก็บเป็น long long ก่อนคูณ ถ้าพิกัดติดลบได้ ผลต่างเองก็เกิน int ดูแบบฝึกหัดข้อ 8
กรณีฐานคือเหลือไม่เกินสามต้น เทียบทุกคู่ตรงๆ ใช้สามแทนหนึ่งเพื่อไม่ให้มีช่วงที่เหลือต้นเดียว ซึ่งไม่มีคู่ให้เทียบ ช่วงที่ยาวสี่ขึ้นไปแบ่งแล้วได้อย่างน้อยฝั่งละสองต้นเสมอ
เรียงแถบตาม y ด้วย lambda แบบหน้า 16.3
ต้นถัดไปทุกต้นสูงกว่านี้ ห่างตามแนว y มากกว่านี้ จึงหยุดได้ทันที ไม่ใช่แค่ข้าม
pair เรียงตาม x ก่อน ถ้า x เท่ากันเรียงตาม y ตามหน้า 16.4

closest squared distance = 5
checks = 10

วิธีลองทุกคู่ต้องเทียบ 8 × 7 / 2 = 28 คู่ วิธีนี้เทียบ 10 คู่ ข้อมูลน้อยแบบนี้ยังต่างกันไม่มาก แต่ตัวนับ checks โตเกือบเป็นเส้นตรงตาม N ตามตารางข้างล่าง

เวลาที่ใช้

แต่ละชั้นของต้นไม้ สร้างแถบใช้ O(n) เรียงแถบใช้ O(n log n) และเทียบในแถบใช้ O(n) เพราะแต่ละต้นเทียบไม่เกิน 7 ต้น มี log₂ N ชั้น รวมเป็น O(N log² N) ตามแถวที่ห้าของตารางในหน้า 22.1 วัดแบบเดียวกับหน้า 22.1 พิกัดสุ่มจาก mt19937 ระหว่าง 0 ถึง 10⁹ จับเวลารวมการเรียงตาม x ด้วย

วิธี	N	Code::Blocks (`-Wall`)	`-O2`	checks
ลองทุกคู่ O(N²)	10,000	0.410 วินาที	0.074 วินาที	ราว 5 × 10⁷
ลองทุกคู่ O(N²)	100,000	27.906 วินาที (5 ครั้งอยู่ระหว่าง 24.2 ถึง 36.7)	4.936 วินาที	ราว 5 × 10⁹
แบ่งครึ่ง O(N log² N)	100,000	0.097 วินาที	0.027 วินาที	142,609
แบ่งครึ่ง O(N log² N)	1,000,000	1.429 วินาที	0.364 วินาที	1,216,103

ทั้งสองวิธีได้ระยะกำลังสอง 28528400 เท่ากันที่ N เป็นแสน แต่วิธีแบ่งครึ่งเทียบแค่ราว 1.4 คู่ต่อต้น น้อยกว่าขีด 7 ที่พิสูจน์ไว้มาก เพราะข้อมูลสุ่มมักมีเสาในแถบไม่กี่ต้น ส่วนเวลาที่เหลือส่วนใหญ่คือการเรียง ถ้าอยากได้ O(N log N) จริงๆ ให้เรียงตาม y ระหว่างทางขากลับแบบ merge sort ในหน้า 22.2 แทนการ sort แถบใหม่ทุกชั้น แต่ N ถึงสองแสนแบบนี้ก็ทันแล้ว

ลองทำให้พัง

แบบที่ 1 ลืมคู่ที่ข้ามเส้น

คิดว่าการเรียกตัวเองกับสองครึ่งครอบคลุมทุกคู่แล้ว จึงคืนค่าที่น้อยกว่าของสองฝั่งเลย

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

vector<pair<int, int>> towers;

long long dist2(pair<int, int> p, pair<int, int> q)
{
    long long dx = p.first - q.first;
    long long dy = p.second - q.second;
    return dx * dx + dy * dy;
}

long long closest(int lo, int hi)
{
    long long best = LLONG_MAX;
    if (hi - lo <= 3) {
        for (int i = lo; i < hi; i++) {
            for (int j = i + 1; j < hi; j++) {
                best = min(best, dist2(towers[i], towers[j]));
            }
        }
        return best;
    }
    int mid = (lo + hi) / 2;
    return min(closest(lo, mid), closest(mid, hi));
}

int main()
{
    int n;
    cin >> n;
    towers.resize(n);
    for (int i = 0; i < n; i++) {
        cin >> towers[i].first >> towers[i].second;
    }
    sort(towers.begin(), towers.end());
    cout << "closest squared distance = " << closest(0, n) << '\n';
    return 0;
}

closest squared distance = 72

ได้ 72 แทนที่จะเป็น 5 และไม่ได้ 18 ของครึ่งซ้ายด้วย เพราะคู่ (5, 30) กับ (8, 27) ก็ข้ามเส้นแบ่งในชั้นที่ลึกลงไป ถ้าไม่ดูคู่ข้ามเส้น คำตอบขึ้นกับว่าคู่ที่ใกล้ที่สุดบังเอิญอยู่ฝั่งเดียวกันในทุกชั้นหรือไม่ ขั้นรวมของ divide and conquer คือส่วนที่ต้องคิดให้ครบที่สุด เหมือนคู่ข้ามครึ่งในหน้า 22.3

แบบที่ 2 ลืมเรียงแถบตาม y

ในโปรแกรม closest-pair ลบบรรทัด sort(strip.begin(), ...) กับ lambda ออก แถบจึงเรียงตาม x ตามที่ได้มาจาก towers ข้อมูลเข้าชุดเดิมได้ closest squared distance = 18 แถบชั้นบนสุดเรียงตาม x คือ (21, 27) (22, 13) (23, 28) ต้นแรกเทียบกับ (22, 13) ผลต่างของ y คือ 14 ยกกำลังสองได้ 196 ไม่น้อยกว่า 18 จึงหยุด ไม่เคยได้เทียบกับ (23, 28) ที่อยู่ถัดไป คำสั่ง break ถูกต้องก็ต่อเมื่อต้นที่ตามมาทุกต้นสูงกว่า ซึ่งจริงเฉพาะเมื่อเรียงตาม y แล้ว ลองสุ่มชุดเล็ก 200 ชุดเทียบกับวิธีลองทุกคู่ ผิด 40 ชุด เขียนวิธีลองทุกคู่ไว้เทียบกับข้อมูลสุ่มเล็กๆ จึงสำคัญกับโจทย์แบบนี้

สรุป

เทียบด้วยระยะกำลังสองที่เป็นจำนวนเต็ม เก็บใน long long ไม่ต้องใช้ทศนิยม
เรียงตาม x แบ่งที่ต้นกลาง หาคำตอบสองฝั่งได้ d แล้วดูเฉพาะแถบที่ (x - x_กลาง)² < d
เรียงแถบตาม y เทียบกับต้นถัดไปจนผลต่างของ y ยกกำลังสองถึง d แต่ละต้นเทียบไม่เกิน 7 ต้น
รวมเป็น O(N log² N) บนเครื่องที่วัด แสนต้นใช้ 0.027 วินาทีแบบ -O2 แทน 4.9 วินาทีของวิธีลองทุกคู่
ขั้นรวมต้องดูคู่ข้ามเส้นแบ่งเสมอ และทดสอบกับวิธีลองทุกคู่ด้วยข้อมูลสุ่มเล็กๆ

กลับไปเช็กรายการ ผ่านบทนี้เมื่อ แล้วไปทำแบบฝึกหัด