22.4 quicksort และ quickselect

merge sort แบ่งง่าย (ตัดกลางเฉยๆ) แต่รวมยาก (ต้องมีอาเรย์พัก) quicksort กลับกัน ทำงานหนักตอนแบ่ง โดยย้ายตัวที่น้อยไปไว้ซ้ายและตัวที่มากไปไว้ขวา พอแบ่งเสร็จก็ไม่ต้องรวมอะไรอีก และไม่ต้องใช้อาเรย์พัก ราคาที่ต้องจ่ายคือถ้าแบ่งไม่ดี มันช้าลงเป็น O(N²) หน้านี้จะทำให้มันช้าให้เห็นจริงกับข้อมูลที่เรียงมาแล้ว แก้ด้วยการสุ่ม แล้วต่อยอดเป็น quickselect ที่หาตัวที่เล็กเป็นอันดับ k ได้โดยไม่ต้องเรียงทั้งหมด ซึ่งก็คือสิ่งที่ nth_element ทำ

แบ่งรอบตัวแบ่ง

เลือกค่าหนึ่งในช่วงมาเป็น ตัวแบ่ง (pivot) แล้วจัดช่วงใหม่ให้ทุกตัวที่น้อยกว่าตัวแบ่งอยู่ทางซ้าย ตัวแบ่งอยู่ตรงกลาง ตัวที่เหลืออยู่ทางขวา ขั้นนี้เรียกว่า partition พอเสร็จ ตัวแบ่งอยู่ในช่องที่มันต้องอยู่ตอนเรียงเสร็จแล้วพอดี เพราะมีตัวที่น้อยกว่ามันอยู่ทางซ้ายครบทุกตัว

วิธีที่เขียนง่ายที่สุดใช้ตัวสุดท้ายของช่วงเป็นตัวแบ่ง แล้วใช้ตัวแปร store จำว่าตัวที่น้อยกว่าตัวแบ่งตัวถัดไปจะต้องไปอยู่ช่องไหน ไล่ i จากซ้ายไปขวา เจอตัวที่น้อยกว่าตัวแบ่งก็สลับไปไว้ที่ store แล้วขยับ store ไปหนึ่งช่อง ทุกช่องก่อน store จึงน้อยกว่าตัวแบ่งเสมอ จบแล้วสลับตัวแบ่งมาไว้ที่ store คะแนนชุดเดียวกับหน้า 22.2 ตัวแบ่งคือ 55

ห้ารอบแรกสลับกับตัวเอง เพราะยังไม่เจอตัวที่ไม่น้อยกว่า 55 เลย พอเจอ 82 แล้ว store ค้างอยู่ที่ช่องของ 82 ตัวที่น้อยกว่าตัวต่อไปคือ 10 จึงสลับกับ 82 จบแล้ว 55 อยู่ช่อง 6 ทางซ้ายมีหกตัวที่น้อยกว่าทั้งหมด ทางขวามี 82 ตัวเดียว partition ไล่ช่วงครั้งเดียว ใช้เวลา O(n)

quicksort

หลัง partition ตัวแบ่งเข้าที่แล้ว เหลือเรียงทางซ้ายกับทางขวาของมันแยกกัน ซึ่งก็ใช้ quicksort เหมือนเดิม และไม่ต้องรวมอะไร เพราะทุกตัวทางซ้ายน้อยกว่าทุกตัวทางขวาอยู่แล้ว ฟังก์ชันที่ทำ partition ตั้งชื่อว่า place_pivot เพราะ partition เป็นชื่อฟังก์ชันที่มีอยู่แล้วใน STL

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

int place_pivot(vector<int> &a, int lo, int hi)
{
    int pivot = a[hi - 1];
    int store = lo;
    for (int i = lo; i < hi - 1; i++) {
        if (a[i] < pivot) {
            swap(a[i], a[store]);
            store++;
        }
    }
    swap(a[store], a[hi - 1]); // (1)!
    return store;
}

void quick_sort(vector<int> &a, int lo, int hi)
{
    if (hi - lo <= 1) {
        return;
    }
    int p = place_pivot(a, lo, hi);
    cout << "[" << lo << ", " << hi << ") pivot " << a[p] << " -> index " << p << ":";
    for (int t = lo; t < hi; t++) {
        cout << ' ' << a[t];
    }
    cout << '\n';
    quick_sort(a, lo, p); // (2)!
    quick_sort(a, p + 1, hi);
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n;
    cin >> n;
    vector<int> a(n);
    for (int i = 0; i < n; i++) {
        cin >> a[i];
    }
    quick_sort(a, 0, n);
    for (int i = 0; i < n; i++) {
        cout << a[i] << (i == n - 1 ? '\n' : ' ');
    }
    return 0;
}

1. swap ใน STL สลับค่าของตัวแปรสองตัวที่ส่งเข้าไป ไม่ต้องเขียนเองแบบหน้า 13.2
2. ช่อง p เข้าที่แล้ว จึงไม่อยู่ในช่วงทั้งสองฝั่ง ทุกครั้งที่เรียกต่อ ช่วงจึงสั้นลงอย่างน้อยหนึ่งช่อง

8
38 27 43 3 9 82 10 55

[0, 8) pivot 55 -> index 6: 38 27 43 3 9 10 55 82
[0, 6) pivot 10 -> index 2: 3 9 10 38 27 43
[0, 2) pivot 9 -> index 1: 3 9
[3, 6) pivot 43 -> index 5: 38 27 43
[3, 5) pivot 27 -> index 3: 27 38
3 9 10 27 38 43 55 82

ต่างจาก merge sort ที่งานอยู่ในขากลับ quicksort ทำงานในขาลง พิมพ์ก่อนแล้วค่อยเรียกตัวเอง และช่วงที่ได้ไม่เท่ากัน ครั้งแรกตัวแบ่ง 55 แบ่ง 8 ตัวเป็น 6 กับ 1 ไม่ใช่ 4 กับ 4 ถ้าตัวแบ่งตกตรงกลางทุกครั้ง ต้นไม้จะเหมือน merge sort คือ log₂ N ชั้น ชั้นละ O(N) รวม O(N log N) แต่ถ้าตัวแบ่งเป็นตัวมากที่สุดทุกครั้ง แต่ละชั้นตัดออกได้แค่ตัวเดียว เป็นการแบ่งไม่ครึ่งแบบในหน้า 22.1 ต้นไม้ลึก N ชั้น และเทียบรวม (N - 1) + (N - 2) + ... + 1 ราว N² / 2 ครั้ง

ข้อมูลที่เรียงแล้วคือกรณีแย่ที่สุด

ตัวแบ่งเป็นตัวมากที่สุดทุกครั้งเกิดเมื่อไร คำตอบคือเมื่อข้อมูลเรียงมาแล้ว ตัวสุดท้ายของทุกช่วงคือตัวมากที่สุดของช่วงนั้นพอดี ข้อมูลที่เรียงแล้วหรือเกือบเรียงเจอบ่อยมากในชุดทดสอบ เช่น ชุดที่คนตั้งโจทย์ตั้งใจใส่ไว้ดักวิธีนี้ วัด quicksort ข้างบน (ตัดส่วนที่พิมพ์ออก) แบบเดียวกับหน้า 22.1 กับข้อมูลสุ่มจาก mt19937 และข้อมูล 0 1 2 ... N - 1 ที่เรียงแล้ว

ข้อมูลสุ่มเร็วพอๆ กับ merge sort ในหน้า 22.2 แต่ข้อมูลที่เรียงแล้ว N เพิ่มสองเท่า เวลาเพิ่มราวสี่เท่า คือโตแบบ N² แค่สองหมื่นตัวก็ช้ากว่าการเรียงล้านตัวสุ่มแล้ว และแบบ Code::Blocks ที่สี่หมื่นตัวยังไม่ทันได้ช้า ก็พังด้วย stack overflow ก่อน เพราะลึกหนึ่งชั้นต่อหนึ่งตัวแบบหน้า 9.4 แบบ -O2 ยังไม่พังที่แสนตัว แต่ก็ใช้เกือบสี่วินาที และห้ามพึ่งว่า -O2 จะรอดตามที่หน้า 9.1 เตือนไว้

สุ่มตัวแบ่ง

ปัญหาไม่ได้อยู่ที่ข้อมูล แต่อยู่ที่เราเลือกตัวแบ่งแบบที่คนอื่นเดาได้ ถ้าสุ่มช่องที่จะเป็นตัวแบ่ง ข้อมูลเรียงมาแบบไหนก็ไม่มีผล แก้แค่ต้นฟังก์ชัน place_pivot สุ่มช่อง r ในช่วง แล้วสลับมาไว้ท้ายช่วงก่อน ส่วนที่เหลือเหมือนเดิมทุกบรรทัด

mt19937 pick(2569); // (1)!

int place_pivot(vector<int> &a, int lo, int hi)
{
    int r = lo + pick() % (hi - lo);
    swap(a[r], a[hi - 1]);
    int pivot = a[hi - 1];
    int store = lo;
    for (int i = lo; i < hi - 1; i++) {
        if (a[i] < pivot) {
            swap(a[i], a[store]);
            store++;
        }
    }
    swap(a[store], a[hi - 1]);
    return store;
}

1. ตัวสุ่มประกาศนอกฟังก์ชัน ถ้าประกาศข้างในจะเริ่มใหม่จาก seed เดิมทุกครั้งที่เรียก ได้เลขชุดเดิมซ้ำๆ อย่าใช้ rand() ซึ่งบน Windows ได้เลขไม่เกิน 32767 สุ่มช่องในช่วงที่ยาวกว่านั้นไม่ได้ครบ

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

sort ของ C++ ข้างในทำอะไร

sort ของ g++ ใช้วิธีที่เรียกว่า introsort คือ quicksort ที่เลือกตัวแบ่งจากค่ากลางของสามตัว และคอยนับความลึกไว้ ถ้าลึกเกินราว 2 log₂ N ชั้น แปลว่ากำลังเจอกรณีแย่ ก็เปลี่ยนไปใช้ heap sort ซึ่งเป็น O(N log N) เสมอ ส่วนช่วงที่เหลือไม่เกิน 16 ตัวใช้ insertion sort จาก 13.3 เพราะเร็วกว่ากับข้อมูลน้อย sort จึงไม่มีวันช้าเป็น O(N²) heap sort อาศัยโครงสร้างที่ชื่อ heap ซึ่งเป็นที่มาของ priority_queue ในหน้า 20.5 และจะเขียนเองในบทที่ 26 สรุปวิธีเรียงทั้งหมดที่เจอมาได้ตามตาราง

quickselect หาอันดับที่ k

นักกีฬาวิ่ง 100 เมตร 9 คน เวลาเป็นร้อยละของวินาที (1342 คือ 13.42 วินาที) อยากรู้เวลาของคนที่เร็วเป็นอันดับ 3 เพื่อใช้เป็นเกณฑ์เข้ารอบ เรียงทั้งหมดก็ได้คำตอบ แต่เราไม่ได้สนใจลำดับของคนอื่นเลย partition บอกอะไรมากกว่าที่คิด หลังแบ่งรอบตัวแบ่งที่ไปอยู่ช่อง p ตัวแบ่งเป็นตัวที่เล็กเป็นอันดับ p (นับจาก 0) แน่นอน

• ถ้า p เท่ากับอันดับที่ต้องการ ได้คำตอบเลย
• ถ้าอันดับที่ต้องการน้อยกว่า p คำตอบอยู่ทางซ้าย ทิ้งทางขวาทั้งหมด
• ถ้ามากกว่า คำตอบอยู่ทางขวา ทิ้งทางซ้ายทั้งหมด

เหลือส่วนเดียวทุกครั้ง เป็น decrease and conquer แบบ binary search ในหน้า 22.1 จึงเขียนเป็นลูปได้ไม่ต้องเรียกตัวเอง ใช้ place_pivot แบบสุ่มตัวแบ่ง

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

mt19937 pick(2569);

int place_pivot(vector<int> &a, int lo, int hi)
{
    int r = lo + pick() % (hi - lo);
    swap(a[r], a[hi - 1]);
    int pivot = a[hi - 1];
    int store = lo;
    for (int i = lo; i < hi - 1; i++) {
        if (a[i] < pivot) {
            swap(a[i], a[store]);
            store++;
        }
    }
    swap(a[store], a[hi - 1]);
    return store;
}

int kth_smallest(vector<int> &a, int k) // (1)!
{
    int lo = 0, hi = a.size();
    while (true) {
        int p = place_pivot(a, lo, hi);
        cout << "[" << lo << ", " << hi << ") pivot " << a[p] << " -> index " << p << '\n';
        if (p == k) {
            return a[p];
        }
        if (k < p) {
            hi = p;
        } else {
            lo = p + 1;
        }
    }
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n, k;
    cin >> n >> k;
    vector<int> times(n);
    for (int i = 0; i < n; i++) {
        cin >> times[i];
    }
    int answer = kth_smallest(times, k - 1); // (2)!
    cout << "rank " << k << " time = " << answer << '\n';
    return 0;
}

1. k นับจาก 0 ลูปจบแน่นอน เพราะ k อยู่ในช่วง [lo, hi) เสมอ และช่วงสั้นลงอย่างน้อยหนึ่งทุกรอบ ถ้าช่วงเหลือช่องเดียว p ก็ต้องเท่ากับ k
2. โจทย์นับอันดับจาก 1 ส่วนช่องนับจาก 0 จึงลบหนึ่ง

9 3
1342 1287 1405 1199 1350 1263 1311 1228 1376

[0, 9) pivot 1350 -> index 6
[0, 6) pivot 1311 -> index 4
[0, 4) pivot 1228 -> index 1
[2, 4) pivot 1263 -> index 2
rank 3 time = 1263

รอบแรกสุ่มได้ 1350 ไปอยู่ช่อง 6 แปลว่า 1350 เร็วเป็นอันดับ 7 อันดับ 3 (ช่อง 2) จึงอยู่ทางซ้าย ทิ้งสองตัวทางขวาไปเลย ทำต่อจนตัวแบ่งตกช่อง 2 พอดี ได้ 12.63 วินาที เรียงด้วยมือจะได้ 1199 1228 1263 ... ตรงกัน mt19937 ให้เลขชุดเดียวกันทุกเครื่องเมื่อ seed เท่ากัน ผลจึงเหมือนกันทุกครั้งที่รัน

ถ้าตัวแบ่งตกแถวกลางๆ ช่วงจะหดราวครึ่งหนึ่งทุกรอบ งานรวมคือ N + N/2 + N/4 + ... ไม่เกิน 2N เป็น O(N) โดยเฉลี่ย ตามแถวที่สองของตารางในหน้า 22.1 ต่างจากการเรียงที่ต้อง O(N log N) เพราะทิ้งได้ทีละครึ่ง ไม่ต้องทำทั้งสองฝั่ง

nth_element

STL มี quickselect ให้ใช้ชื่อ nth_element(begin, nth, end) หลังเรียก ช่องที่ nth ชี้จะมีค่าเดียวกับที่ช่องนั้นจะมีถ้าเรียงทั้งหมด ทุกช่องก่อนหน้าไม่มากกว่า และทุกช่องหลังจากนั้นไม่น้อยกว่า แต่ภายในแต่ละฝั่งไม่เรียง ถ้าต้องการ k ตัวที่เล็กที่สุดแบบเรียงด้วย ให้ sort เฉพาะ k ช่องแรกต่อ ใช้เวลาเพิ่ม O(k log k) เท่านั้น

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

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n, k;
    cin >> n >> k;
    vector<int> times(n);
    for (int i = 0; i < n; i++) {
        cin >> times[i];
    }
    nth_element(times.begin(), times.begin() + (k - 1), times.end()); // (1)!
    cout << "rank " << k << " time = " << times[k - 1] << '\n';
    sort(times.begin(), times.begin() + k); // (2)!
    cout << "fastest " << k << ":";
    for (int i = 0; i < k; i++) {
        cout << ' ' << times[i];
    }
    cout << '\n';
    return 0;
}

1. อันดับ k นับจาก 1 อยู่ช่อง k - 1
2. เรียงเฉพาะ k ช่องแรก ก่อนบรรทัดนี้ k ช่องแรกเป็น k ตัวที่เร็วที่สุดอยู่แล้ว แต่ลำดับข้างในไม่แน่นอน

9 3
1342 1287 1405 1199 1350 1263 1311 1228 1376

rank 3 time = 1263
fastest 3: 1199 1228 1263

วัดการหาค่ามัธยฐาน (k = N / 2) ของเลขสุ่มจาก mt19937 สามวิธี วัดแบบเดียวกับหน้า 22.1 และวัดรอบเดียวกับตารางในหน้า 22.2 แถว sort ที่ล้านตัวจึงได้เท่ากัน

ทั้งสามวิธีได้ค่ามัธยฐานเท่ากัน nth_element เร็วกว่า sort ราวสิบห้าเท่าที่สิบล้านตัว ของที่เขียนเองก็เร็วระดับเดียวกัน แต่ขึ้นกับดวงของการสุ่มตัวแบ่ง ตารางใช้ pick(2569) ซึ่งได้ตัวแบ่งชุดเดิมทุกครั้ง ลองเปลี่ยน seed เป็น 1 ถึง 10 ที่สิบล้านตัวแบบ -O2 ได้ตั้งแต่ 0.043 ถึง 0.066 วินาที เมื่อโจทย์ต้องการแค่ค่ามัธยฐาน อันดับที่ k หรือ k ตัวที่เล็กที่สุด ใช้ nth_element ได้เลย

ลองทำให้พัง

แบบที่ 1 ข้อมูลซ้ำกันเยอะ

เอา quicksort แบบสุ่มตัวแบ่งไปเรียงข้อมูลที่ทุกตัวเท่ากัน เช่น คะแนนเต็มทั้งห้อง วัดแบบเดียวกับตารางข้างบน

การสุ่มไม่ช่วยเลย เพราะ partition เทียบด้วย a[i] < pivot ตัวที่เท่ากับตัวแบ่งไปอยู่ทางขวาทั้งหมด สุ่มได้ตัวไหนก็เท่ากันหมด ทางซ้ายจึงว่างทุกครั้ง แบ่งไม่ครึ่งแบบเดียวกับข้อมูลที่เรียงแล้ว ข้อมูลที่มีค่าไม่กี่แบบ เช่น เกรด 0 ถึง 4 ของนักเรียนแสนคน ก็เจอปัญหาเดียวกัน

แก้ด้วยการแบ่งเป็นสามกลุ่ม น้อยกว่า เท่ากับ และมากกว่าตัวแบ่ง แล้วเรียกตัวเองเฉพาะกลุ่มน้อยกว่ากับกลุ่มมากกว่า กลุ่มที่เท่ากันเข้าที่แล้วทั้งกลุ่ม ใช้ตัวชี้สามตัว ทุกช่องก่อน lt น้อยกว่าตัวแบ่ง ช่อง lt ถึงก่อน i เท่ากับตัวแบ่ง ช่องตั้งแต่ gt มากกว่าตัวแบ่ง ส่วน i ถึงก่อน gt คือช่องที่ยังไม่ได้ดู

void quick_sort(vector<int> &a, int lo, int hi)
{
    if (hi - lo <= 1) {
        return;
    }
    int pivot = a[lo + pick() % (hi - lo)];
    int lt = lo, i = lo, gt = hi;
    while (i < gt) {
        if (a[i] < pivot) {
            swap(a[i], a[lt]);
            lt++;
            i++;
        } else if (a[i] > pivot) {
            gt--;
            swap(a[i], a[gt]); // (1)!
        } else {
            i++;
        }
    }
    quick_sort(a, lo, lt);
    quick_sort(a, gt, hi);
}

1. ไม่เลื่อน i เพราะค่าที่สลับมาจากช่อง gt ยังไม่เคยถูกดู

ข้อมูลที่เท่ากันหมดล้านตัวใช้ 0.002 วินาทีแบบ Code::Blocks และไม่ถึง 0.001 วินาทีแบบ -O2 เพราะรอบแรกทุกตัวเข้ากลุ่มเท่ากับ แล้วจบเลย ส่วนข้อมูลสุ่มล้านตัวใช้ 0.447 กับ 0.080 วินาที ข้อมูลที่เรียงแล้วล้านตัวใช้ 0.277 กับ 0.035 วินาที ยังเป็น O(N log N) ตามเดิม แต่แบบ Code::Blocks ช้ากว่าแบบสองกลุ่มราวเท่าตัวเพราะเทียบมากกว่า ชุดทดสอบของโจทย์เรียงข้อมูลควรมีทั้งข้อมูลสุ่ม ข้อมูลที่เรียงแล้ว และข้อมูลที่ซ้ำกันหมด

แบบที่ 2 ส่งอันดับให้ nth_element ผิดไปหนึ่ง

ในโปรแกรม nth-element เขียน times.begin() + k และอ่าน times[k] แทน k - 1 ป้อนข้อมูลชุดเดิมได้ rank 3 time = 1287 ซึ่งเป็นเวลาของอันดับ 4 คอมไพล์ผ่านโดยไม่มีคำเตือน และถ้าป้อน k เท่ากับ N คือ 9 จะอ่านช่องที่ 9 ซึ่งอยู่นอก vector บนเครื่องที่เขียนเอกสารนี้ได้ rank 9 time = 0 ทุกครั้งที่ใช้อันดับจากโจทย์ ให้ถามตัวเองว่านับจาก 0 หรือ 1 และทดสอบด้วย k = 1 กับ k = N เสมอ

สรุป

• partition จัดให้ตัวที่น้อยกว่าตัวแบ่งอยู่ซ้าย ตัวแบ่งเข้าที่ของมันพอดี ใช้เวลา O(n) quicksort เรียกตัวเองกับสองฝั่งโดยไม่ต้องรวม
• ตัวแบ่งตัวสุดท้ายกับข้อมูลที่เรียงแล้วคือกรณีแย่ที่สุด O(N²) บนเครื่องที่วัด สี่หมื่นตัวแบบ Code::Blocks พังด้วย stack overflow
• สุ่มตัวแบ่งด้วย mt19937 แล้วแบ่งสองกลุ่ม ได้ O(N log N) โดยเฉลี่ยเมื่อค่าไม่ซ้ำ ไม่ว่าข้อมูลจะเรียงมาแบบไหน ข้อมูลซ้ำเยอะต้องแบ่งสามกลุ่ม
• sort คือ quicksort ที่สลับไปใช้ heap sort เมื่อลึกเกิน จึงเป็น O(N log N) เสมอ ส่วน merge sort ใช้เมื่อต้องการ stable หรือทำงานเพิ่มระหว่างรวม
• quickselect ทิ้งฝั่งที่ไม่มีคำตอบทุกรอบ ได้อันดับที่ k ใน O(N) โดยเฉลี่ย ใน STL คือ nth_element ระวังอันดับนับจาก 0 หรือ 1