Haskell的クイックソート

この記事はHaskell Advent Calendar 2022の4日目の記事です。

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

Wikipediaのおはなし

一時期、Wikipedia.jaの「クイックソート」の項に、Haskellでの「実装例」が示されていました。（⇒ Wikipedia）
「ちょっとこのコードはどうだろう」と思って様子見していたら、「適切でない実装例を削除します。（他モジュールの利用、節々の計算重複、関数型言語で表現する意味を感じられないアイディア、etc.）」という編集コメントと共にざっくり削除されてしまいました。実装例はC言語だけになって、それさえ書いてあれば用は足りるか、という感じになりました。

という昔話をしようと現状を確認したら、C言語に加えて、SchemeとPythonまで書いてあります。いやもうPython実装があるならC実装なくてよくね？ていうか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))))))

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

(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:

(define (qsort f ls)

qsort という関数を定義しているらしい。Haskellなら

qsort f ls = ...

行2,3:

  (if (null? ls)
      '()

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

qsort _ []     = []

行4:

      (let ((x (car ls)) (xs (cdr ls)))

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

qsort f (x:xs) = ...

行5:

        (let ((before (qsort f (filter (lambda (y) ((compose not f) x y)) xs)))

比較関数fにnotを合成して逆にして、xsをxとの逆の比較でfilterした結果をqsortに再帰的に渡した結果をbeforeに束縛するらしい。Haskellなら

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

行6:

              (after  (qsort f (filter (lambda (y) (f x y)) xs))))

同じく、fでxと比較してfilterした結果をqsortしてafterに束縛する。Haskellなら

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

行7:

          (append before (cons x after))))))

beforeと「afterにxをコンスしたもの」を連結したら答えらしい。Haskellなら

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

まとめると

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 を使うことにする。
仕切り直して最初の定義を示す。

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

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

３つに分割する

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

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

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 の初期値は空リストでなくこれにしてしまえばいい。

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要素とすれば、「クイックソートを打ち切る条件」を「ピボットの候補を確保できない場合」と共用できる。

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に掲載してはいけない。