N 個のデータのうち小さい方から数えて k 番目の値を求める方法を「選択アルゴリズム (Selection Algorithm)」と呼びます。簡単な方法はクイックソートなどでソートすることですが、データ数 N が多くなるとクイックソートでも N * log2 N に比例する時間が必要になります。ところが、クイックソートを変形することで、データ数 N に比例する程度の時間で k 番目のデータを求めることができます。これを「クイックセレクト (quick select)」と呼びます。
クイックソートはある値（枢軸）を基準にして、要素をそれより大きいものと小さいものの 2 つに分割していくことでソートを行います。2 つに分けた各々のグループを同様に分割して、さらに 2 つのグループに分けます。分割を続けていくと最後は要素がひとつになりソートが完了します。
クイックセレクトも枢軸を基準に二分割していきますが、枢軸の決め方と二分割したグループの扱い方が異なるのです。
３番目のデータを求める
↓
table[９５３７６４２８] ３番目に格納されている７を
↑ ↑ 枢軸にして二分割する
Ｌ Ｒ
table[２５３４６７９８] Ｒ＞３なので前半部分０－４に対して
↑↑ 同様の操作を行う
ＲＬ
↓
table[２５３４６７９８] 今度は４を枢軸として二分割する
↑ ↑
Ｌ Ｒ
table[２３５４６７９８] Ｌ＜３なので後半部分２－４に対して
↑↑ 同様の操作を行う
ＲＬ
図 : クイックセレクトのアルゴリズム
上図を見てください。配列 table の中から 3 番目のデータを求めます。この場合、枢軸は table[3] に格納されている 7 となります。この値を基準に、クイックソートと同様に配列を二分割すると、0 から R の前半部分に枢軸より小さい値が集まり、L から最後までの後半部分には枢軸より大きい値が集まります。
もし、R が 3 より大きいのであれば、求めるデータは前半部分に含まれることになります。逆に L が 3 より小さいのであれば、後半部分に含まれます。上図の場合は、R は 4 なので前半部分に求めるデータがあるのです。この部分に対して、同様の処理を繰り返します。
今度も table[3] に格納されている要素を枢軸にします。この枢軸に対して配列を二分割します。この処理を繰り返していくと、最後にはグループに含まれる要素がひとつになって分割不能になります。このとき、table[3] に格納されたデータが答となるのです。
このように、クイックセレクトは分割したグループの片方に対してのみ操作を行うので、クイックソートよりも高速に処理することが可能になります。
プログラムは次のようになります。
リスト : クイックセレクト
search1 :: Ord a => IOArray Int a -> a -> Int -> IO Int
search1 buff pv i = do
x <- readArray buff i
if pv > x
then search1 buff pv (i + 1)
else return i
search2 :: Ord a => IOArray Int a -> a -> Int -> IO Int
search2 buff pv j = do
x <- readArray buff j
if pv < x
then search2 buff pv (j - 1)
else return j
quickPartition :: Ord a => IOArray Int a -> a -> Int -> Int -> IO (Int, Int)
quickPartition buff pv low high = do
i <- search1 buff pv low
j <- search2 buff pv high
if i < j
then do swapItem buff i j
quickPartition buff pv (i + 1) (j - 1)
else return (i, j)
quickSelect :: Ord a => IOArray Int a -> Int -> IO a
quickSelect buff k = do
(start, end) <- getBounds buff
qselect start end
where
qselect start end
| start >= end = readArray buff k
| otherwise = do
pv <- readArray buff k
(i, j) <- quickPartition buff pv start end
if j < k
then qselect (j + 1) end
else qselect start (i - 1)
引数 k が求める順番を表し、局所関数 qselect の引数 start と end が分割する区間の両端を表します。関数 search1, search2, partition は「配列のソート (2)」で作成したクイックソートと同じです。
最初に、配列 table の k 番目に位置するデータを枢軸 pivot に選び、partition で区間を二分割します。分割が終わって j が k より小さい場合は後半部分を分割します。区間 (j + 1, end) に対して qselect を再帰呼び出しします。そうでなければ前半部分を分割します。区間 (start, i - 1) に対して qselect を再帰呼び出しします。start が end 以上になったならば buff の k 番目の値を返します。
それでは、どの程度効果があるのか拙作のページ「配列のソート (2)」で作成したクイックソート quickSort と比較してみましょう。データ数を 5000, 10000, 20000, 40000 増やしていき、その中央値を求めることにします。テストプログラムは簡単なので説明は割愛します。詳細はプログラムリストをお読みください。
インタプリタ ghci 上で実行したところ、実行結果は次のようになりました。
ghci> test 5000 159050967812799482 0.030705092s 159050967812799482 0.146068966s ghci> test 10000 159755821740510529 0.050423603s 159755821740510529 0.307371815s ghci> test 20000 42943112069563441 0.115262209s 42943112069563441 0.649589248s ghci> test 40000 88201593115835707 0.214998141s 88201593115835707 1.349078594s
select がとても速いことがわかります。
ご参考までに、リストで動作するプログラムも作ってみました。この場合、枢軸はリストの先頭要素とし、分割した前半のリストの長さ a で場合分けします。a が k と等しい場合、枢軸が k 番目の要素になります。a が k よりも大きい場合、k 番目の要素は前半のリストにあります。それ以外の場合、k 番目の要素は後半のリストにあります。興味のある方はプログラムリストをお読みくださいませ。
--
-- quickselect.hs : クイックセレクト
--
-- Copyright (C) 2013-2021 Makoto Hiroi
--
import Data.Array.IO
import Control.Monad
import Data.Time
import System.Random
-- 要素の交換
swapItem :: IOArray Int a -> Int -> Int -> IO ()
swapItem buff i j = do
a <- readArray buff i
b <- readArray buff j
writeArray buff i b
writeArray buff j a
-- ソート
search1 :: Ord a => IOArray Int a -> a -> Int -> IO Int
search1 buff pv i = do
x <- readArray buff i
if pv > x
then search1 buff pv (i + 1)
else return i
search2 :: Ord a => IOArray Int a -> a -> Int -> IO Int
search2 buff pv j = do
x <- readArray buff j
if pv < x
then search2 buff pv (j - 1)
else return j
quickPartition :: Ord a => IOArray Int a -> a -> Int -> Int -> IO (Int, Int)
quickPartition buff pv low high = do
i <- search1 buff pv low
j <- search2 buff pv high
if i < j
then do swapItem buff i j
quickPartition buff pv (i + 1) (j - 1)
else return (i, j)
quickSort :: Ord a => IOArray Int a -> IO ()
quickSort buff = do
(low, high) <- getBounds buff
qsort low high
where
qsort low high = do
pv <- readArray buff ((low + high) `div` 2)
(i, j) <- quickPartition buff pv low high
when (low < i - 1) $ qsort low (i - 1)
when (high > j + 1) $ qsort (j + 1) high
-- 選択
quickSelect :: Ord a => IOArray Int a -> Int -> IO a
quickSelect buff k = do
(start, end) <- getBounds buff
qselect start end
where
qselect start end
| start >= end = readArray buff k
| otherwise = do
pv <- readArray buff k
(i, j) <- quickPartition buff pv start end
if j < k
then qselect (j + 1) end
else qselect start (i - 1)
-- リスト版
quickSort' :: Ord a => [a] -> [a]
quickSort' [] = []
quickSort' (x:xs) = quickSort' [y | y <- xs, y < x]
++ [x]
++ quickSort' [y | y <- xs, y >= x]
quickSelect' :: Ord a => [a] -> Int -> a
quickSelect' xs k =
iter xs k
where iter [x] 0 = x
iter (x:xs) k = if a == k then x
else if a > k then iter ys k
else iter zs (k - a - 1)
where ys = [y | y <- xs, y < x]
a = length ys
zs = [z | z <- xs, z >= x]
test :: Int -> IO ()
test n = do
a <- newListArray (0, n - 1) (take n (randoms (mkStdGen 11) :: [Int])) :: IO (IOArray Int Int)
b <- newListArray (0, n - 1) (take n (randoms (mkStdGen 11) :: [Int])) :: IO (IOArray Int Int)
x1 <- getCurrentTime
y1 <- quickSelect a (n `div` 2)
x2 <- getCurrentTime
print y1
print (diffUTCTime x2 x1)
x3 <- getCurrentTime
quickSort b
y2 <- readArray b (n `div` 2)
x4 <- getCurrentTime
print y2
print (diffUTCTime x4 x3)
test' :: Int -> IO ()
test' n = do
let a = take n (randoms (mkStdGen 11) :: [Int])
b = take n (randoms (mkStdGen 11) :: [Int])
x1 <- getCurrentTime
print (quickSelect' a (n `div` 2))
x2 <- getCurrentTime
print (diffUTCTime x2 x1)
x3 <- getCurrentTime
let c = quickSort' b
print (c !! (n `div` 2))
x4 <- getCurrentTime
print (diffUTCTime x4 x3)