# Pattern-defeating Quicksort 閱讀筆記
###### tags: `Algorithm`
完整的 python 實作可參考 [Github 連結](https://github.com/Chang-Chia-Chi/pdqsort) [10]
## Introduction
在論文撰寫的時間點，最常使用的複合(hybrid)排序算法應該是內省排序 (intro sort) [2]，內省排序由插入排序 (insertion sort)、堆排序 (heap sort) 以及快速排序 (quick sort) 組合使用而成，當遞迴深度過深時採用堆排序，若分區的尺寸低於某個大小時則改用插入排序。
:::success
堆排序:
確保算法的時間複雜度能維持在 O(NlogN)

插入排序:
節省函式遞迴呼叫的成本，降低記憶體的使用量 (省去 logX 的遞迴棧幀 (stack frame)，以 C++ 為例 X = 16 [2]) 並且提昇快取的 locality
:notes:
:::

Pattern-defeating Quicksort (下簡稱 **pdqsort**) 受到內省排序的啟發，保有快速排序節省記憶體空間且在多數情況高效的特性，另外在某些特定情況下，時間複雜度可壓縮到線性時間，這也是為什麼叫做 pattern-defeating 的原因。當然無可避免地，快速排序不穩定的性質也顯現在 pdqsort 上，因此不適用於需要穩定排序的場景。

## Overview
- 同樣由快速排序、插入排序以及堆排序組合而成
- 使用三向分段法 (tripartite partitioning) 搭配無分支 (branchless) 分段 Block partitioning [3] (註：為求簡單，本筆記參考 [4] 僅使用三向分段)
- 基於 partition badness 作為遞迴深度的依據，論文引入 bad partition 的概念，當偵測到某次分段狀況不佳時，會交換 (swap) 部份元素，以引入更高的隨機性，但當 bad partition 發生次數過多時，代表數據很有可能本身就不是均勻分佈的，將轉換為推排序
- 採用 median-of-3 ninther [5] 選擇鉚點 (pivot) $ninther(A) = med3(med3(A[1 ...\frac{1}{3}n]), med3(A[\frac{1}{3}n ... \frac{2}{3}n]), med3(A[\frac{2}{3}n ... n]))$
實作如下：

    def _choosePivot(cls, left, right, size, arr):
        def sort3(a, b, c):
            if arr[b] < arr[a]: a, b = b, a
            if arr[c] < arr[b]: b, c = c, b
            if arr[b] < arr[a]: a, b = b, a
            return b

        half = size // 2
        if size > PDQSort.NINTHER_THRESHOLD:
          a = sort3(left + 1, left + half - 1, right - 1)
          b = sort3(left, left + half, right)
          c = sort3(left + 2, left + half + 2, right - 2)
          return sort3(a, b, c)
        else:
          return sort3(left + half, left, right)

:::info
無分支的比較函式，通常有兩種方式 [3]
1. 條件轉移 (a < b)? a : b (condition move，x86 指令 CMOVcc)
2. 將布林值轉換為整數 1 or 0 (x86 指令 SETcc)
:::

## Dutch national flag problem (荷蘭國旗問題)
快速排序在基數越小 (low cardinality) 時，時間複雜度越容易退化成 $O(n^2)$，而要妥善處理資料中的相同數據 (equal comparing elements)，我們需要使用三向分段，而三向分段非常類似荷蘭國旗問題

荷蘭國旗問題 [6] 是 Dijkstra 提出的演算法問題，該問題其實是三色球的排序與分類，只是剛好 Dijkstra 的家鄉荷蘭，國旗也是三色的，就如此命名了。如果單純用快速排序，時間複雜度很高機率會是最糟糕的 $O(n^2)$，而改用如合併排序這種穩定的排序方法，Dijkstra 老人家好像也沒有必要提出這個問題，所以要思考的是如何在 $O(n)$ 時間完成三色球的排序

假設我們將三色球分別編號成 `0` `1` `2`，我們要做的事情是將 `0` 放到 `1` 的左邊，將 `2` 放到 `1` 的右邊，這不就是以 `1` 為鉚點的分段嗎？只是這次我們會需要 3 個指針：
- `left`： 指向陣列左方的元素
- `curr`： 指向目前迭代的元素
- `right`： 指向陣列右方的元素

當 `curr` 指向的元素為 `0` 時，代表我們應該要把該元素放到陣列左側，因此交換 `left` & `curr`，同時移動兩個指針 ; 同理，當 `curr` 指向 `2` 時，我們交換 `right` & `curr`，但是注意這時候我們僅移動右指針，因為交換之後 `curr` 有可能指向 `0` ; 最後如果 `curr` 指向 `1` 的話，我們直接移動 `curr` 指針即可，因為此時在 `curr` 左邊的元素理應都小於等於 `1`。實作可參考下方程式碼

def sortColors(nums):
    curr, left, right = 0, 0, len(nums)-1
    while curr <= right:
        if nums[curr] == 0:
            nums[curr], nums[left] = nums[left], nums[curr]
            curr += 1
            left += 1
        elif nums[curr] == 2:
            nums[curr], nums[right] = nums[right], nums[curr]
            right -= 1
        else:
            curr += 1

:::warning
順帶一題，LeetCode 75 題 [7] 就是在考荷蘭國旗問題 :lower_left_paintbrush:
:::

回到我們的問題上，對於快速排序來說，會造成效率低落的原因之一，就是有許多重複的元素造成分段不平衡，導致某一側的遞迴深度變得非常深，而且每個遞迴呼叫並沒有排序到多少元素。我們雖然不能直接套用荷蘭國旗問題的算法，但可以借鑒這個算法的思想：在分段的過程中，將與鉚點相同的元素集中在陣列的中間，如下圖所示


一個直覺的作法是，我們先將與鉚點相同的元素換到左側或右側，最後再將這些元素換回中間，但這個方法的缺點是，除了與鉚點比較大小之外，每個元素在交換前，都需確認是否與鉚點相同，導致額外的比較成本

為了解決這個問題，論文使用兩個不同的分段函式 **partition_left** & **partition_right**，**partition_left** 將小於等於鉚點的元素，放到鉚點的左側，而 **partition_right** 將大於等於鉚點的元素，放到鉚點的右側。假設我們剛開始進行排序，並且找到鉚點 $p$，剩餘元素統稱為 $\alpha$，整個分段的步驟如下：
1. 對 $\alpha$ 執行 **partition_right**，將大於等於 $p$ 的元素放在右側
2. 找到新的鉚點 $q$
3. 如果 $p\,!= q$，則對 $\beta$ 執行 **partition_right**，假如往右遞迴的過程中，$p$ 都不等於 $q$，代表 $p$ 右側的元素都大於 $p$
4. 如果 $p = q$，改執行 **partition_left**，我們知道 $∀x ∈ β : x ≥ p$，因此 $∄x ∈ β : x < q$，當呼叫 **partition_left** 時，會將小於等於 $q$ 的元素丟到左側，但從上述的推論可以得知， $\beta'$ 的元素一定不小於 $q$，所以分段後 $\beta'$ 左側的元素一定等於 $p, q$，右側的元素一定大於 $p, q$
:::success
簡單來說，論文採用的三向分段方法，會呼叫 N 次 **partition_right** + 1 or 0 次 **partition_left**。**partition_left** 不會往下遞迴呼叫，因為該區間的元素都等於 $q$。好處就是省去比較以及紀錄元素是否與鉚點相同的成本，並且將時間複雜度壓在 $O(nk)$，詳情請參閱後續的證明
:notes:
:::

參考論文及 [4]，**partition_right** 程式碼實作如下，使用雙指針的技巧，從陣列左側及右側向中間夾，當 `i` 指向的元素大於等於鉚點且 `j` 指向的元素小於鉚點時進行交換。**partition_left** 的實作方式近乎雷同就不放上來了

def partitionRight(cls, left, right, pivot_idx, array):
    """
    Partitions array within [left, right] with array[pivotIdx].
    Elements larger or equal to pivot will be put on the right-hand
    side of the array.

    @params left: left bound of partition
    @params right: right bound of partition
    @params pivot_idx: index of pivot in array
    @params array: array to be sorted

    @return (next_pivot_idx, no_swap): (next pivot index, true if range [left, right] of array is already partitioned)
    """
    pivot = array[pivot_idx]
    array[left], array[pivot_idx] = array[pivot_idx], array[left]

    i, j = left + 1, right
    while i <= j and array[i]  < pivot: i += 1
    while i <= j and array[j] >= pivot: j -= 1

    no_swap = i > j

    while i <= j:
        array[i], array[j] = array[j], array[i]
        while array[i]  < pivot: i += 1
        while array[j] >= pivot: j -= 1

    array[left], array[j] = array[j], array[left]
    return j, no_swap

:::warning
**no_swap** 的功用會在稍後的章節講解
:::

### An $O(nk)$ worst case of pdqsort with k distinct elements
*註：$O(nk)$ 為單純使用快速排序的上限
#### Lemma 1. 任何子序列*的 predecessor 都是其 ancestor 的鉚點
如果某個子序列是其父序列的右子序列，則其 predecessor (子序列的前一個元素) 會是父序列的鉚點，但如果是左子序列的話，其 predecessor 會是父序列的 predecessor


註：只考慮被傳入遞迴呼叫的子序列，代表該子序列不是最左子序列，一定存在 predecessor

#### Lemma 2. 當某個 distinct value $v$ 被選為鉚點，一定不同於 predecessor
假設 $v$ 等於 predecessor，但從 Lemma 1 我們知道 predecessor 是其祖先序列的鉚點，與 Lemma 1 矛盾，因此 $v$ 一定不同於 predecessor
:::info
由 Lemma 2 可知，當 $v$ 第一次被選做鉚點時，其一定搭配 **partition_right** 進行分段
:::

#### Lemma 3. 直到與 $v$ 相同的元素再次被選為鉚點，所有相同於 $v$ 的元素 $x$ 都位於 $v$ 的右側
由於我們以 $v$ 執行 **partition_right**，任何其他分段 $w$ < $v$ 與 $x$ 都無關。如果以 $w$ > $v$ 進行分段，也不可能將 $x$ 放到 $w$ 的右側，因為 $x ≥ w > v$ 與分段的定義產生矛盾。由此得證所有 $x=v$ 的元素都會被放在 $v$ 的右側

#### Lemma 4. 當我們第 2 次使用與 $v$ 相同的元素為鉚點時，所有 $x=v$ 的元素都已經在正確的位置上
當 $x=v$ 被選做鉚點時，由 Lemma 3 可知 $v$ 一定是其 predecessor，所以會使用 **partition_left** 進行分段，所有與 $x$ ($v$) 相同的元素會被分在其左側，代表在這一次分段步驟中，所有與 $v$ 相同的元素已經處理完成並擺放在正確的位置上

#### Theorem 1. 當陣列中有 k 種不同的元素，pdqsort 的時間複雜度為 $O(nk)$
Lemma 4 說明所有 distinct value 可以被選做鉚點至多兩次，而且兩次分段後所有與該鉚點的元素都已被擺在正確的位置上。每一次分段時間複雜度 $O(n)$，假如陣列有 $k$ 種不同的元素，則最差的時間複雜度為 $O(nk)$

## Other novel techniques
### Preventing quicksort’s $O(n^2)$ worst case
**pdqsort** 以 **bad partition** 發生的次數來決定切換成堆排序的時機。當鉚點將陣列平均切分時，我們稱之為完美分段 (perfect partition)，鉚點的百分位數 $p=0.5$，當 $p$ 偏離 0.5 越多代表這次的分段品質越差。但怎樣才稱之為 **bad** 呢？我們可以從資訊熵的角度出發


論文採用 $p=0.125$ 作為 **bad partition** 的閥值，一來是 0.125 在任何平台上，都可以用簡單的位元運算處理，若我們將上式畫圖可以發現，當 $p$ 約等於 0.125 時，效率約為完美分段的 1/2，這也說明為什麼快速排序在多數情況下，有不錯的表現


:::success
作者在文中表明，使用 **bad partition** 作為判斷基準比遞迴深度要好上許多。原因是很多時候，內省排序雖然一開始會有比較差的分段，但分了幾次之後通常這種情況會被打破，如果使用遞迴深度為基準，很可能在元素分佈其實適合使用快速排序的狀況下，採用效率較差的堆排序
:::

不免俗地，還是要證明一下 **pdqsort** 最差的時間複雜度為 $O(nlogn)$

#### Lemma 5. bad partitions 的花費時間至多 $O(nlogn)$
若某個遞迴子樹遇到 $logn$ 次 **bad partitions**，將轉換為堆排序。每一層的時間複雜度最差為 $O(n)$，因此時間複雜度限制在 $O(nlogn)$

#### Lemma 6. good partitions 的花費時間至多 $O(nlogn)$
假若每次分段都剛好形成 **最差** 的 good partition，則其遞迴式可表達如下


由 Akra-Bazzi [9] 定理，$Θ(T (n, p)) = Θ(n log n)$

#### Theorem 2. **pdqsort** 的時間複雜度為 $O(nlogn)$
由 Lemma 5 & Lemma 6 可以證明 **pdqsort** 最糟的時間複雜度為 $O(nlogn)$

:::info
在 Python 中，可以使用 `bit_length()` 來計算 **bad partition** 發生次數的上限

limit = len(arr).bit_length()

另外要特別注意，這個數值並非全局變數，而是每個遞迴子樹的局部變數
:::
### Introducing fresh pivot candidates on bad partitions
作者發現，假如前一次的分段是 **bad partition**，有很高機率下一次算法依然會挑到不好的鉚點，作者將這個現象稱之為自相似性 **self similarity**。論文提出的解決方案是將頭、尾的元素，跟 1/4、3/4 位置的元素進行交換

def _shufflePattern(cls, left, right, pivot_idx, arr):
    left_size, right_size = pivot_idx - left, right - (pivot_idx + 1)
    cls._shuffle(left, pivot_idx, left_size, arr)
    cls._shuffle(pivot_idx + 1, right, right_size, arr)

def _shuffle(cls, start, end, size, arr):
if size >= PDQSort.INSERTION_SORT_THRESHOLD:
    arr[start], arr[start + size//4] = arr[start + size//4], arr[start]
    arr[end], arr[end - size//4] = arr[end - size//4], arr[end]

    if size > PDQSort.NINTHER_THRESHOLD:
        arr[start + 1], arr[start + (size//4 + 1)] = arr[start + (size//4 + 1)], arr[start + 1]
        arr[start + 2], arr[start + (size//4 + 2)] = arr[start + (size//4 + 2)], arr[start + 2]
        arr[end - 1], arr[end - size//4 - 1] = arr[end - size//4 - 1], arr[end - 1]
        arr[end - 2], arr[end - size//4 - 2] = arr[end - size//4 - 2], arr[end - 2]

經過實測發現，洗牌策略確實能明顯提昇分段的品質，測試的程式碼可參考 [10]，以 [10] 的測試案例來說，大概可以從 $p = 0.1$ 提昇至 $p = 0.25$
:::success
那為什麼這個簡單的邏輯有用呢？個人推測當發生 **bad partition** 時，表示目前 median of 3 的分佈，讓鉚點很接近陣列的某一側，又因為分段的不平衡，導致下一次仍然挑選到不好的鉚點。如果要打破這種情況，那我們改變 median of 3 的分佈就好啦，而從遠離鉚點的位置引入不同元素，就是一個暴力、簡單而且有效的方法
:::

## Previously known techniques
### Insertion sort
所有快速排序的最佳化，都會在陣列長度小於某個數值時 (16 ~ 32)，轉為其他穩定的排序算法，通常採用插入排序。參考論文的實作，使用移動元素的方式排序，減少不必要的元素交換。此外，本文採用 **binary search**，以 $O(logn)$ 時間找到插入點，降低元素大小的比較次數

class InsertionSort:
    @classmethod
    def sort(cls, left, right, array):
        for i in range(left + 1, right + 1):
            if array[i] < array[i - 1]:
                tmp = array[i]
                insert_idx = bisect.bisect_right(array, array[i], left, i)

            j = i
            while j >= insert_idx:
                array[j] = array[j - 1]
                j -= 1
            array[insert_idx] = tmp

### Optimistic insertion sorting on a swapless partition
當我們在某次分段發現除了鉚點之外，不用交換任何其他元素，我們稱這次分段為 **swapless**。如果 **swapless** 及 **good partition** 同時發生，代表這次的分段有可能已經將陣列排序的差不多，此時我們轉用部份插入排序 (partial insertion sort) 對左右兩側子序列進行排序。假若我們妥善選擇鉚點，則部份插入排序在升序、降序或升序陣列尾巴附帶一些任意元素的情況下，能夠在 $O(n)$ 時間完成排序動作

:::warning
**partial** 的意思是，我們會追蹤這次排序交換了多少次元素，如果超過某個閥值 (論文取 8)，則直接結束插入排序的動作。若低於閥值的話，代表我們已經排序好左右子序列，不用繼續往下遞迴
:::

部份插入排序的實作與插入排序雷同，只多了追蹤閥值的部份

@classmethod
def partial_sort(cls, left, right, array):
    limit = 0
    for i in range(left + 1, right + 1):
        if array[i] < array[i - 1]:
            tmp = array[i]
            insert_idx = bisect.bisect_right(array, array[i], left, i)

            j = i
            while j >= insert_idx:
                array[j] = array[j - 1]
                j -= 1
            array[insert_idx] = tmp
            limit += i - insert_idx

        if limit > InsertionSort.PARTIAL_INSERTION_SORT_LIMIt: return False
    return True

### Block partitioning
本筆記的實作並未採用，細節請自行參閱 [3]

## Reference
[1. Pattern-defeating Quicksort](https://arxiv.org/pdf/2106.05123.pdf)
[2. 內省排序- 維基百科，自由的百科全書](https://zh.wikipedia.org/zh-tw/%E5%86%85%E7%9C%81%E6%8E%92%E5%BA%8F)
[3. BlockQuicksort: How Branch Mispredictions don’t affect Quicksort](https://arxiv.org/pdf/1604.06697.pdf)
[4. zhangyunhao116/pdqsort](https://github.com/zhangyunhao116/pdqsort/tree/master)
[5. Median - Wikipedia](https://en.wikipedia.org/wiki/Median#Ninther)
[6. Dutch national flag problem - Wikipedia](https://en.wikipedia.org/wiki/Dutch_national_flag_problem)
[7. Leetcode - Sort Colors](https://leetcode.com/problems/sort-colors/)
[8. Go1.19で採用された Pattern-defeating Quicksort の紹介](https://speakerdeck.com/po3rin/go1-dot-19decai-yong-sareta-pattern-defeating-quicksort-falseshao-jie)
[9. Notes on Better Master Theorems for Divide-and-Conquer Recurrences](https://courses.csail.mit.edu/6.046/spring04/handouts/akrabazzi.pdf)
[10. Chang-Chia-Chi/pdqsort](https://github.com/Chang-Chia-Chi/pdqsort)