Haskell的クイックソート

この記事は[Haskell Advent Calendar 2022](https://qiita.com/advent-calendar/2022/haskell)の4日目の記事です。

何番煎じかもわからないネタですが、お付き合いください。

# Wikipediaのおはなし

一時期、Wikipedia.jaの「クイックソート」の項に、Haskellでの「実装例」が示されていました。[（⇒ Wikipedia）](https://ja.wikipedia.org/w/index.php?title=%E3%82%AF%E3%82%A4%E3%83%83%E3%82%AF%E3%82%BD%E3%83%BC%E3%83%88&oldid=48066057#%E5%AE%9F%E8%A3%85%E4%BE%8B%EF%BC%92)
「ちょっとこのコードはどうだろう」と思って様子見していたら、「適切でない実装例を削除します。（他モジュールの利用、節々の計算重複、関数型言語で表現する意味を感じられないアイディア、etc.）」という編集コメントと共にざっくり削除されてしまいました。実装例はC言語だけになって、それさえ書いてあれば用は足りるか、という感じになりました。

という昔話をしようと現状を確認したら、C言語に加えて、SchemeとPythonまで書いてあります。いやもうPython実装があるならC実装なくてよくね？ていうかScheme実装を見てみると…

```scheme
(require-extension srfi-1)              ; SRFI 1 の呼び出し方は実装依存(場合によってはデフォルト)。これは Chicken Scheme の例。

(define (qsort f ls)
  (if (null? ls)
      '()
      (let ((x (car ls)) (xs (cdr ls)))
        (let ((before
               (qsort f (filter (lambda (y) ; filter は SRFI 1 の手続き
                                  ;; compose は Chicken、Gauche、Guile、Racket 等に備わってる「合成関数を作る」手続き。
                                  ;; 詳細は Paul Graham の On Lisp
                                  ;; http://www.asahi-net.or.jp/~kc7k-nd/onlispjhtml/returningFunctions.html
                                  ;; を参照。
                                  ((compose not f) x y)) xs)))
              (after
               (qsort f (filter (lambda (y) ; filter は SRFI 1 の手続き
                                  (f x y)) xs))))
          (append before (cons x after))))))
```

妙に長くて言い訳くさいコメントを消して改行調整して

```scheme
(define (qsort f ls)
  (if (null? ls)
      '()
      (let ((x (car ls)) (xs (cdr ls)))
        (let ((before (qsort f (filter (lambda (y) ((compose not f) x y)) xs)))
              (after  (qsort f (filter (lambda (y) (f x y)) xs))))
          (append before (cons x after))))))
```

前から読んでいく。
#### 行1:
```scheme
(define (qsort f ls)
```
 `qsort` という関数を定義しているらしい。Haskellなら

```haskell
qsort f ls = ...
```

#### 行2,3:
```scheme
  (if (null? ls)
      '()
```

`ls`が空リストなら空リストを返すらしい。Haskellなら

```haskell
qsort _ []     = []
```

#### 行4:
```scheme
      (let ((x (car ls)) (xs (cdr ls)))
```

`ls`が空リストでないとき、`ls`の`car`と`cdr`に`x`と`xs`という名前を付けるらしい。Haskellなら

```haskell
qsort f (x:xs) = ...
```

#### 行5:
```scheme
        (let ((before (qsort f (filter (lambda (y) ((compose not f) x y)) xs)))
```
比較関数`f`に`not`を合成して逆にして、`xs`を`x`との逆の比較で`filter`した結果を`qsort`に再帰的に渡した結果を`before`に束縛するらしい。Haskellなら

```haskell
    before = qsort f (filter (not . f x) xs)
```

#### 行6:
```scheme
              (after  (qsort f (filter (lambda (y) (f x y)) xs))))
```
同じく、`f`で`x`と比較して`filter`した結果を`qsort`して`after`に束縛する。Haskellなら

```haskell
    after  = qsort f (filter (      f x) xs)
```

#### 行7:
```scheme
          (append before (cons x after))))))
```
`before`と「`after`に`x`をコンスしたもの」を連結したら答えらしい。Haskellなら

```haskell
qsort f (x:xs) = before ++ x : after
```

#### まとめると
```haskell
qsort :: (a -> a -> Bool) -> [a] -> [a]
qsort _ []     = []
qsort f (x:xs) = before ++ x : after
  where
    before = qsort f (filter (not . f x) xs)
    after  = qsort f (filter (      f x) xs)
```
これHaskellの教科書の最初の方に必ず書いてある（ダメな）クイックソートそのものじゃないですか。削除された方のHaskell版はともかく、Schemeのそれを載せるよりもこの教科書的なHaskell版を載せる方がよくないですか…
（少なくとも、C,Python実装がしている、ピボット候補を3点からサンプルしている動作がないことは注記するべきでは。）

# Haskell的効率化

ここから本題。
比較関数のことはとりあえず忘れて、`Ord` の `compare` を使うことにする。
仕切り直して最初の定義を示す。

```haskell
import Data.List (partition)

quickSort :: Ord a => [a] -> [a]
quickSort []     = []
quickSort (p:xs) = quickSort as ++ p : quickSort cs
  where
    (as, cs) = partition (< p) xs
```

リストの先頭をピボットとし、それ未満 `as` とそれ以上 `cs` に分け、再帰的にソートした結果を繋ぎ合わせる、クイックソートの分割統治な側面を説明する究極に簡潔なコード。
（配列の左と右から、ピボット以上とピボット未満を探してスワップする in-place な計算がクイックソートの本質だ、という命令型パラダイム厨はお帰り下さい。）

これがよくないと言われる理由は：

1. `(++)` を多用しているので遅い。
1. ピボットの選択が単純で、最悪計算量に陥りやすい。
1. 特に、`cs` の中に `p` と等しい値が紛れ込む。これがピボットに選ばれるたびに `as` は空になって再帰がスカる。

これらに対処していこう。

## ３つに分割する

`partition` で「ピボット未満、ピボット以上」の2つに分割する代わりに、`compare` を用いて「未満、等しい、より大きい」に3分割してみる。
これで、問題点3の、一度使用したピボットが再度ピボットとして選ばれることがなくなる。

```haskell
quickSort []     = []
quickSort (p:xs) = quickSort as ++ p : bs ++ quickSort cs
  where
    (as, bs, cs) = foldr (step p) ([],[],[]) xs

step p x (as,bs,cs) =
  case compare p x of
    LT -> (as, bs, x:cs)
    EQ -> (as, x:bs, cs)
    GT -> (x:as, bs, cs)
```

## `(++)` をなくす

問題点1について、`(++)` をなくすには、リストを返す計算に対して、その結果の続きに `(++)` されるリストを追加の引数で渡し、これの前に要素を付け足す形に変形することが定番のやり方。

```haskell
quickSort xs = qs xs []  -- qs xs rest は、xsを整列した結果をrestの前に連結する
  where
    qs []     rest = rest
    qs (p:xs) rest = qs as (p : bs ++ qs cs rest)     -- (※1)
      where
        (as, bs, cs) = foldr (step p) ([],[],[]) xs   -- (※2)

step p x (as,bs,cs) = -- 略
```
(※1)の箇所、`as`の整列結果と`p:bs`を連結する`(++)`が消えた。`cs`の整列結果と`rest`を連結する`(++)`は元々いらない。しかし`p:bs`と`cs`の整列結果を連結する`(++)`が残ってしまった。

ここで(※2)の `(as,bs,cs)` を求める計算に注目する。空リストから始めて、ピボットと比較した結果によって3タプルのどこかに付け足すことで、`as`,`cs`にピボット未満、ピボットより大きい値を集める。これらは再帰的に整列する対象となる。
しかし `bs` については、集めたらそれ以上することはない。また、`bs` の後続は `qs cs rest` の結果である。つまり、`foldr` における `bs` の初期値は空リストでなくこれにしてしまえばいい。

```haskell
quickSort xs = qs xs []
  where
    qs []     rest = rest
    qs (p:xs) rest = qs as (p:bs)
      where
        (as, bs, cs) = foldr (step p) ([], qs cs rest, []) xs

step p x (as,bs,cs) = -- 略
```
遅延評価だからできる技。
`cs` が `foldr` の式の両辺に出現する、初心者にはとても見せられない定義になってしまった。

## ピボットを3か所から採取する + 短いときは選択ソートに切り替える

常に一定の位置からピボットを選ぶと、例えば昇順または降順に整列されたリストに対して、最悪計算量になる。これを回避するために、いくつかのピボット候補から、大きくも小さくもない中間の値を選ぶことで、ピボットとしての有効性を高めるという方法で問題点2に対処する。

そのため例えば、リストの先頭 `(!! 0)` に加えて、17番目の値 `(!! 16)` と33番目の値 `(!! 32)` を候補として取り出すことが考えられる。しかしそのためには、入力が33要素あることが必要である。

ところで、クイックソートは要素数が少ないときには効率がよくないため、再帰が進んで要素が減ったときは選択ソートなどに切り替えるという手法が有効であるという。この切り替えの閾値を32要素とすれば、「クイックソートを打ち切る条件」を「ピボットの候補を確保できない場合」と共用できる。

```haskell
quickSort xs = qs xs []
  where
    qs xs rest
      | null xs32 = foldr insert rest xs -- 挿入ソート
      | otherwise = qs as bs
      where
        xs32 = drop (halflimit * 2) xs
        p = med3 (head xs) (xs !! halflimit) (head xs32)
        (as, bs, cs) = foldr (step p) ([], qs cs rest, []) xs

halflimit = 16

-- 3つの値の中央を選ぶ
med3 a b c
  | c <= x    = x   -- c - x - y の順
  | y <= c    = y   -- x - y - c の順
  | otherwise = c   -- x - c - y の順
  where
    (x,y) = minMax a b  -- aとbは x - y の順

-- min と max を同時に
minMax a b = if a <= b then (a,b) else (b,a)

step p x (as,bs,cs) = -- 再掲
  case compare p x of
    LT -> (as, bs, x:cs)
    EQ -> (as, x:bs, cs)
    GT -> (x:as, bs, cs)
```

## おしまい。

こうして、Haskellらしい改造を加えたクイックソートが完成したが、こんなものを「これがHaskellで書くクイックソートでござい」とWikipediaに掲載してはいけない。