M.Hiroi's Home Page

お気楽 Haskell プログラミング入門

応用編 : 配列のソート (1)

Copyright (C) 2013-2021 Makoto Hiroi
All rights reserved.

はじめに

今回は配列の簡単な例題として、一次元配列をソートするプログラムを作ってみましょう。ソートはある規則に従ってデータを順番に並べることです。たとえば、データが整数であれば、昇順もしくは降順に並べます。ソートは昔から研究されている分野で、優秀なアルゴリズムが確立しています。

今回作成するプログラムは、配列の添字を 0 から始まる整数 (Int) とし、データも整数 (Int) とします。この場合、boxed タイプの配列 IOArray よりも unboxed タイプの配列 IOUArray のほうが高速になります。プログラムの説明は IOArray で行いますが、実行結果は IOArray と IOUArray の両方を示すことにします。

データは乱数のほかにも、山型データ (中央のデータがいちばん大きく端にいくほど小さいデータになる)、ソート済みのデータ (昇順)、逆順にソートされたデータ (降順) の 4 種類を用意します。プログラムは GHC (ghc -O) でコンパイルして、実行時間を比較してみましょう。

●boxed と unboxed の違い

最初に、boxed と unboxed の違いについて簡単に説明します。Haskell の場合、boxed タイプの配列はデータを直接格納しているわけではありません。データ自体は他のメモリ領域に格納しておいて、そのメモリ領域の場所 (アドレス) を配列に格納している、と考えてください。実際には、boxed タイプの配列は遅延評価されるので、もっと複雑な処理になります。

これに対し、unboxed タイプの配列はデータを配列のメモリ領域に直接格納します。格納できるデータ型は限られていますが、メモリの消費量は boxed タイプよりも少なくなります。また、unboxed タイプの配列は正格評価され、配列のアクセスも boxed タイプより高速に行うことができます。

IOArray の unboxed タイプが IOUArray です。モジュール Data.Array.IO をインポートすれば IOUArray を利用することができます。IOUArray の基本的な操作関数は IOArray と同じです。

●ソートの安定性

次に、ソートの安定性について説明します。安定 (stable) なソートとは、ソートキーが等しい場合、入力された順番が崩れないソートのことです。逆に、不安定なソートとは、入力された順番とは異なる結果になるソートのことです。次の図を見てください。

  元のデータ      安定なソート    不安定なソート
 ------------------------------------------------
  (123, 'abc')    (123, 'abc')    (789, 'abc')
  (456, 'def')    (789, 'abc')    (123, 'abc')
  (789, 'abc')    (456, 'def')    (456, 'def')

              図 : ソートの安定性

2 番目の要素をキーにソートした場合、1 行目と 3 行目はソートキーが等しいですね。安定なソートであればソート結果が 1, 3, 2 の順番になり、入力時の位置関係が保たれています。不安定なソートはソート結果が、たとえば 3, 1, 2 の順番となり、入力時の位置関係が崩れます。

実用的なアプリケーションを作成する場合、ソートの安定性が重要になる場合があります。一般に、単純なソート(遅いソート)は安定で、複雑なソート(高速なソート)は安定ではない、という傾向があります。

●テストプログラム

次に、共通で使用するプログラムを作ります。

リスト : 共通で使用するプログラム

-- 要素の比較
compItem :: Ord a => IOArray Int a -> Int -> Int -> IO Ordering
compItem buff i j = liftM2 (compare) (readArray buff i) (readArray buff j)

-- 要素の交換
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

-- イテレータ
iter_ :: (Int, Int) -> (Int -> IO ()) -> IO ()
iter_ (low, high) fn
  | low > high = return ()
  | otherwise  = do
      fn low
      iter_ (low + 1, high) fn

iterR_ :: (Int, Int) -> (Int -> IO ()) -> IO ()
iterR_ (low, high) fn
  | low > high = return ()
  | otherwise  = do
      fn high
      iterR_ (low, high - 1) fn

compItem と swapItem は「ヒープ」で作成したプログラムと同じです。iter_ は low から high までの整数を引数の関数 fn に渡して呼び出します。逆に、iterR_ は high, high - 1, ..., low までの整数を fn に渡して呼び出します。

テストプログラムは次のようになります。

リスト : テストプログラム

test :: (IOArray Int Int -> IO ()) -> Int -> IO ()
test sort n = do
  let m = n `div` 2
      check i ary 
        | i == n    = return ()
        | otherwise = do
            t <- compItem ary (i - 1) i
            if t == GT then error "test error"
            else check (i + 1) ary
  a <- newListArray (0, n - 1) (take n (randoms (mkStdGen 11) :: [Int])) :: IO (IOArray Int Int)
  b <- newListArray (0, n - 1) [1..n] :: IO (IOArray Int Int)
  c <- newListArray (0, n - 1) [n,n-1..1] :: IO (IOArray Int Int)
  d <- newListArray (0, n - 1) ([1..m] ++ [m,m-1..1]) :: IO (IOArray Int Int)
  x1 <- getCurrentTime
  sort a
  x2 <- getCurrentTime
  check 1 a
  print (diffUTCTime x2 x1)
  x3 <- getCurrentTime
  sort b
  x4 <- getCurrentTime
  check 1 b
  print (diffUTCTime x4 x3)
  x5 <- getCurrentTime
  sort c
  x6 <- getCurrentTime
  check 1 c
  print (diffUTCTime x6 x5)
  x7 <- getCurrentTime
  sort d
  x8 <- getCurrentTime
  check 1 d
  print (diffUTCTime x8 x7)

main :: IO ()
main = do
  let xs1 = [2000, 4000, 8000, 16000]
      xs2 = [40000, 80000, 160000, 320000]
  (x:_) <- getArgs
  case x of
    "bubleSort"   -> mapM_ (test bubleSort) xs1
    "bubleSort'"  -> mapM_ (test bubleSort') xs1
    "bubleSort''" -> mapM_ (test bubleSort'') xs1
    "selectSort"  -> mapM_ (test selectSort) xs1
    "insertSort"  -> mapM_ (test insertSort) xs1
    "shellSort"   -> mapM_ (test shellSort) xs2
    "shellSort'"  -> mapM_ (test shellSort') xs2
    "combSort"    -> mapM_ (test combSort) xs2

ファイル名を sort.hs とすると、ghc -O sort でコンパイルします。あとはコマンドラインで sort bubleSort のようにテストするソートを指定してください。実行時間が表示されます。

なお、プログラムの実行時間は、筆者のコーディング、実行したマシン、使用するプログラミング言語(またはコンパイラ)などの環境に大きく依存しています。また、これらの環境だけではなく、データの種類によっても実行時間は大きく左右されます。たとえば、高速なソートアルゴリズムとして有名な「クイックソート」の場合、データによっては実行時間が要素数の 2 乗に比例する最悪のケースが存在します。どのようなデータに対しても最速を誇るソートアルゴリズムは存在しません。興味のある方は、いろいろなデータをご自分の環境で試してみてください。

●バブルソート

最初は、単純なソートから説明しましょう。バブルソート (buble sort) は泡がぶくぶくと浮いてくるように、いちばん小さいデータが後ろから前に浮かび上がってくるアルゴリズムです。

隣接する 2 つのデータを比較して、順序が逆であれば入れ換えます。これを順番に後ろから前に行っていけば、いちばん小さなデータは頂上に浮かび上がっているというわけです。先頭が決まったならば、残りのデータに対して同じことを行えば、2 番目には残りのデータの中でいちばん小さいものが浮かび上がってきます。これをデータ数だけ繰り返せばソートが完了します。

 9 5 3 7 6 4 8   交換しない
           ~~~
 9 5 3 7 6 4 8   交換する
         ~~~
 9 5 3 7 4 6 8   交換する
       ~~~
 9 5 3 4 7 6 8   交換しない
     ~~~
 9 5 3 4 7 6 8   交換する
   ~~~
 9 3 5 4 7 6 8   交換する
 ~~~
 3 9 5 4 7 6 8   いちばん小さいデータが決定する
 +               残りのデータに対して同様な操作を行う

        図 : バブルソート

これをそのままプログラムすると次のようになります。

リスト : バブルソート

bubleSort :: Ord a => IOArray Int a -> IO ()
bubleSort buff = do
  (low, high) <- getBounds buff
  iter_
    (low, high - 1)
    (\i -> iterR_
             (i + 1, high)
             (\j -> do t <- compItem buff j (j - 1)
                       when (t == LT) $ swapItem buff j (j - 1)))

bubleSort のプログラムは簡単です。iter_ で (データの個数 - 1) 回だけラムダ式を呼び出します。ラムダ式の中では、iterR_ で buff の後ろから前に向かって、確定していないデータを比較していき、もしも順番が逆になっていたら交換します。

それでは実行してみましょう。

ghci> a <- newListArray (0,9) [5,6,4,7,3,8,2,9,1,0] :: IO (IOArray Int Int)
ghci> bubleSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> b <- newListArray (0,9) [9,8..0] :: IO (IOArray Int Int)
ghci> bubleSort b
ghci> getElems b
[0,1,2,3,4,5,6,7,8,9]
ghci> c <- newListArray (0,9) [0..9] :: IO (IOArray Int Int)
ghci> bubleSort c
ghci> getElems c
[0,1,2,3,4,5,6,7,8,9]
  表 : buble sort の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.056  0.038  0.044  0.039
 4000 : 0.187  0.132  0.165  0.146
 8000 : 0.758  0.520  0.638  0.581
16000 : 3.107  2.070  2.643  2.405

 [IOUArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.042  0.031  0.034  0.033
 4000 : 0.154  0.114  0.126  0.127
 8000 : 0.606  0.450  0.487  0.473
16000 : 2.408  1.775  1.954  1.871

バブルソートは簡単なアルゴリズムですが、データ数が多くなると時間がかかります。データ数を N とすると実行時間は N の 2 乗に比例します。バブルソートは安定なソートですが遅いアルゴリズムなのです。

●バブルソートの改良

バブルソートは改良すべきポイントがひとつあります。次の例を見てください。

   4 5 6 7 8 9 3   交換する
             ~~~
   ・・省略・・

   3 4 5 6 7 8 9   3 を決定
   +

   3 4 5 6 7 8 9   交換しない
   +         ~~~

   ・・省略・・    一度もデータを交換しない

   3 4 5 6 7 8 9   4 を決定 => ソート終了
   + +

この例の場合、3 を先頭に移動すればソートが完了しますね。 バブルソートの場合、データの交換がまったく行われなければ、 ソートは完了したと判断することができます。 この例では 3 を先頭に移動したあと、さらにバブルソートを続けますが、 データの交換は行われずに次のデータ 4 が決まりますね。 ここでソートを終了してよいのです。この改良を行うと昇順データは極めて高速にソートできるようになります。

プログラムは次のようになります。

リスト : バブルソートの改良

-- IORef を使う方法
bubleSort' :: Ord a => IOArray Int a -> IO ()
bubleSort' buff = do
  (low, high) <- getBounds buff
  loop low high
  where loop low high
          | low >= high = return ()
          | otherwise = do
              f <- newIORef False
              iterR_
                (low + 1, high)
                (\j -> do t <- compItem buff j (j - 1)
                          when (t == LT) $ do 
                            swapItem buff j (j - 1)
                            writeIORef f True)
              x <- readIORef f
              when x $ loop (low + 1) high

-- IORef を使わない方法
bubleSort'' :: Ord a => IOArray Int a -> IO ()
bubleSort'' buff = do
  (low, high) <- getBounds buff
  loop low high
  where loop low high
          | low >= high = return ()
          | otherwise = do
              t <- loop1 (low + 1) high False
              when t $ loop (low + 1) high
        loop1 i j done
          | i > j     = return done
          | otherwise = do
              t <- compItem buff j (j - 1)
              if t == LT
                then do
                  swapItem buff j (j - 1)
                  loop1 i (j - 1) True
                else loop1 i (j - 1) done

bubleSort' は newIORef で書き換え可能な変数を生成し、データの交換が行われたら変数の値を True に書き換えます。bubleSort'' は newIORef で生成した変数のかわりに、局所関数 loop1 の第 3 引数を使います。データの交換が行われたら第 3 引数を True にして loop1 を再帰呼び出しします。あとは loop1 の返り値をチェックして、False であれば loop の再帰呼び出しを行わずに処理を終了します。

それでは実行結果を示します。

  表 : buble sort' の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.050  0.000  0.046  0.031
 4000 : 0.184  0.000  0.169  0.118
 8000 : 0.740  0.000  0.677  0.459
16000 : 3.423  0.000  2.724  1.833

 [IOUArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.056  0.000  0.037  0.027
 4000 : 0.153  0.000  0.133  0.094
 8000 : 0.589  0.000  0.526  0.367
16000 : 2.315  0.000  2.072  1.444
  表 : buble sort'' の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.093  0.000  0.037  0.025
 4000 : 0.158  0.000  0.140  0.095
 8000 : 0.628  0.000  0.557  0.378
16000 : 2.572  0.000  2.229  1.515

 [IOUArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.059  0.000  0.028  0.020
 4000 : 0.118  0.000  0.109  0.074
 8000 : 0.470  0.000  0.414  0.292
16000 : 1.871  0.000  1.650  1.140

どちらの関数も昇順のデータに対しては極めて高速に動作します。他のデータの場合、IORef を使った bubleSort' ではそれほど速くなりませんでした。bubleSort'' は改良の効果により bubleSort よりも速くなりました。Haskell は純粋な関数型言語なので、IORef の処理は時間がかかるのかもしれません。

●選択ソート

選択ソート (selection sort) は、ソートしていないデータの中から最小値(または最大値)を見つけ、それを先頭のデータと交換する、という手順を繰り返すことでソートを行います。最初は、すべてのデータの中から最小値を探し、それを配列の先頭 buff[0] と交換します。次は、buff[1] 以降のデータの中から最小値を探し、それを buff[1] と交換します。これを繰り返すことでソートすることができます。

 [9 5 3 7 6 4 8]   3 と 9 を交換する
  +   +

 3 [5 9 7 6 4 8]   5 と 4 を交換する
    +       +

 3 4 [9 7 6 5 8]   9 と 5 を交換する
      +     +

 3 4 5 [7 6 9 8]   7 と 6 を交換する
        + +

 3 4 5 6 [7 9 8]   7 と 7 を交換する
          +

 3 4 5 6 7 [9 8]   9 と 8 を交換してソート終了
            + +

        図 : 選択ソート

このように、選択ソートは単純でわかりやすいアルゴリズムです。プログラムは次のようになります。

リスト : 選択ソート

selectMin buff low high = iter buff (low + 1) low
  where iter buff i n
          | i > high  = return n
          | otherwise = do
              r <- compItem buff i n
              if r == LT
                then iter buff (i + 1) i
                else iter buff (i + 1) n

selectSort :: Ord a => IOArray Int a -> IO ()
selectSort buff = do
  (low, high) <- getBounds buff
  iter_ (low, high - 1) (\i -> do x <- selectMin buff i high
                                  swapItem buff i x)

選択ソートは配列の区間 (low, high) から最小値を求める関数 seletMin を作ると簡単です。実際の処理は局所関数 iter で行います。第 3 引数 n が最小値の位置を表します。compareItem で i 番目と n 番目の要素を比較し、i 番目の要素が小さい場合は第 4 引数を i に変更します。

selectSort は low から high までの中から最小値を selectMin で選び、それと low 番目の要素を swapItem で交換します。これを 0 番目からから high - 1 番目まで繰り返します。

それでは実行結果を示します。

ghci> a <- newListArray (0,9) [5,6,4,7,3,8,2,9,1,0] :: IO (IOArray Int Int)
ghci> selectSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> b <- newListArray (0,9) [9,8..0] :: IO (IOArray Int Int)
ghci> selectSort b
ghci> getElems b
[0,1,2,3,4,5,6,7,8,9]
ghci> c <- newListArray (0,9) [0..9] :: IO (IOArray Int Int)
ghci> selectSort c
ghci> getElems c
[0,1,2,3,4,5,6,7,8,9]
  表 : select sort の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.041  0.030  0.031  0.036
 4000 : 0.130  0.120  0.119  0.121
 8000 : 0.488  0.468  0.470  0.473
16000 : 1.924  1.873  1.858  1.871

 [IOUArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.030  0.030  0.028  0.029
 4000 : 0.104  0.103  0.105  0.104
 8000 : 0.410  0.403  0.401  0.403
16000 : 1.599  1.598  1.596  1.493

選択ソートも実行時間はデータの個数の 2 乗に比例します。選択ソートは遅いソートですが安定です。

●単純挿入ソート

単純挿入ソートはソート済みの配列に新しいデータを挿入していくことでソートします。最初は先頭のデータひとつがソート済みと考え、2 番目のデータをそこに挿入することからスタートします。データを挿入するので、そこにあるデータをどかさないといけません。そこで、挿入位置を決めるため後ろから順番に比較するとき、いっしょにデータの移動も行うことにします。

 [9] 5 3 7 6 4 8    5 を取り出す

 [9] * 3 7 6 4 8    5 を[9]の中に挿入する

 [5 9] 3 7 6 4 8    9 をひとつずらして先頭に 5 を挿入

 [5 9] * 7 6 4 8    3 を取り出して[5 9]の中に挿入する

 [3 5 9] 7 6 4 8    先頭に 3 を挿入

 [3 5 9] * 6 4 8    7 を取り出して[3 5 9] に挿入

 [3 5 7 9] 6 4 8    9 を動かして 7 を挿入
                    残りの要素も同様に行う

        図 : 単純挿入ソート

プログラムは次のようになります。

リスト : 単純挿入ソート

search_move :: Ord a => IOArray Int a -> a -> Int -> Int -> IO Int
search_move buff x i g
  | i < 0 = return (i + g)
  | otherwise = do
      y <- readArray buff i
      if x < y
        then do writeArray buff (i + g) y
                search_move buff x (i - g) g
        else return (i + g)

insertElement :: Ord a => IOArray Int a -> Int -> Int -> IO ()
insertElement buff i gap = do
  tmp <- readArray buff i
  pos <- search_move buff tmp (i - gap) gap
  writeArray buff pos tmp

insertSort :: Ord a => IOArray Int a -> IO ()
insertSort buff = do
  (low, high) <- getBounds buff
  iter_ (low + 1, high) (\i -> insertElement buff i 1)

関数 insertElement は i 番目の要素を 0 から i 番目の間に挿入します。単純挿入ソートの場合、引数 gap には 1 を指定します。シェルソートで再利用するため、引数に gap を用意しました。最初に、i 番目の要素を取り出して変数 tmp にセットします。あとは i - 1 番目から先頭に向かって挿入する位置を探します。挿入位置の探索は関数 search_move で行います。このとき、データの移動も同時に行っていることに注意してください。関数 insertSort は low + 1 番目の要素から順番にデータを挿入していくだけです。

それでは実行結果を示します。

ghci> a <- newListArray (0,9) [5,6,4,7,3,8,2,9,1,0] :: IO (IOArray Int Int)
ghci> insertSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> b <- newListArray (0,9) [9,8..0] :: IO (IOArray Int Int)
ghci> insertSort b
ghci> getElems b
[0,1,2,3,4,5,6,7,8,9]
ghci> c <- newListArray (0,9) [0..9] :: IO (IOArray Int Int)
ghci> insertSort c
ghci> getElems c
[0,1,2,3,4,5,6,7,8,9]
  表 : insert sort の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.009  0.000  0.018  0.009
 4000 : 0.035  0.000  0.073  0.036
 8000 : 0.135  0.000  0.269  0.137
16000 : 0.565  0.000  1.072  0.537

 [IOUArray]
 個数    乱数   昇順   逆順   山型
----------------------------------
 2000 : 0.008  0.000  0.015  0.008
 4000 : 0.026  0.000  0.056  0.026
 8000 : 0.102  0.000  0.222  0.102
16000 : 0.407  0.000  0.811  0.461

挿入ソートも実行時間はデータの個数の 2 乗に比例します。挿入ソートは遅いソートですが安定です。

ところで、単純挿入ソートは興味深い特性をもっています。昇順データの結果をみると、とても速いことがわかります。データがソートされていれば、挿入位置を探索する search_move はすぐに終了するため、最初の iter_ の繰り返し回数でソートが完了することになります。したがって、与えられたデータが大まかにでもソートされていれば、search_move で繰り返す回数が少なくなり、単純挿入ソートでも高速にソートすることができます。

選択ソートと単純挿入ソートを比較すると、データの比較回数は単純挿入ソートのほうが少なくなり(平均すると約半分になる)、データの移動回数は選択ソートの方が少なくなります。今回は単純挿入ソートの方が速くなりましたが、使用するプログラミング言語やデータの種類によっては、結果が異なる場合もあるでしょう。興味のある方はいろいろ試してみてください。

●シェルソート

シェルソート (shell sort) は単純挿入ソートの改良版ともいえる方法です。最初は遠く離れた要素間でソートを開始し、徐々に間隔を狭めていきます。最後は隣り合った要素間でソートします。つまり、単純挿入ソートと同じになります。

間隔が大きいときは要素の個数が少ないので、単純なアルゴリズムでもソートにかかる時間は少なくてすみます。間隔が小さくなると要素の個数は多くなりますが、大まかにソートされているので単純挿入ソートでも高速にソートすることが可能です。

 9 5 3 7 6 4 2 8  最初の状態

 9       6        間隔を 4 で分割する
   5       4
     3       8
       7       2

 6       9        各群をソートする
   4       5
     3       8
       2       7

 6   3   9   8    間隔を 2 で分割する
   4   2   5   7

 3   6   8   9    各群をソートする
   2   4   5   7

 3 2 6 4 8 5 9 7  間隔を 1 で分割する(単純挿入ソートと同じ)

 2 3 4 5 6 7 8 9  ソート完了

        図 : シェルソート

プログラムは次のようになります。

リスト : シェルソート

shellSort :: Ord a => IOArray Int a -> IO ()
shellSort buff = do
  (low, high) <- getBounds buff
  loop1 ((high - low + 1) `div` 2) high
  where loop1 0   _    = return ()
        loop1 gap high = do iter_ (gap, high) (\i -> insertElement buff i gap)
                            loop1 (gap `div` 2) high

局所関数 loop1 で間隔 (gap) を徐々に狭めていきます。ここでは単純に 2 で割っていくことにしました。iter_ では、各群を並列にソートしていることに注意してください。あとは insertElement で i 番目のデータを取り出して、適切な位置に挿入します。このとき、探索は隣り合った要素ではなく gap 離れた要素を比較することに注意してください。最後に gap は 1 になるので、単純挿入ソートと同じになりソートが完了します。

それでは実行結果を示します。

ghci> a <- newListArray (0,9) [5,6,4,7,3,8,2,9,1,0] :: IO (IOArray Int Int)
ghci> shellSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> a <- newListArray (0,9) [0..9] :: IO (IOArray Int Int)
ghci> shellSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> a <- newListArray (0,9) [9,8..0] :: IO (IOArray Int Int)
ghci> shellSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
  表 : shell sort の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
-----------------------------------
 40000 : 0.040  0.016  0.020  0.017
 80000 : 0.102  0.035  0.040  0.036
160000 : 0.221  0.068  0.087  0.076
320000 : 0.577  0.144  0.178  0.160

 [IOUArray]
 個数    乱数   昇順   逆順   山型
-----------------------------------
 40000 : 0.054  0.013  0.015  0.014
 80000 : 0.054  0.027  0.031  0.031
160000 : 0.117  0.057  0.071  0.061
320000 : 0.253  0.114  0.133  0.125

データ数は 20 倍になっていることに注意してください。結果を見ればお分かりのように、シェルソートは速いですね。gap を常に奇数になるようにすると、シェルソートの実行速度はデータの個数 N の 1.5 乗に比例します。また、Knuth によると、gap の値に次の数列を用いると、シェルソートは N の 1.25 乗に比例するそうです。

gap = ..., 121, 40, 13, 4, 1

この数列は 3 倍して 1 を加えることで得られる数列を逆にしたものです。これをプログラムすると、次のようになります。

リスト : シェルソートの改良版

shellSort' :: Ord a => IOArray Int a -> IO ()
shellSort' buff = do
  (low, high) <- getBounds buff
  loop1 (getGap (high - low + 1) 1) high
  where getGap size gap
          | gap < size `div` 9 = getGap size (gap * 3 + 1)
          | otherwise          = gap
        loop1 0   _    = return ()
        loop1 gap high = do iter_ (gap, high) (\i -> insertElement buff i gap)
                            loop1 (gap `div` 3) high
  表 : shell sort (改良版) の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
-----------------------------------
 40000 : 0.053  0.010  0.011  0.012
 80000 : 0.072  0.023  0.027  0.025
160000 : 0.162  0.043  0.056  0.054
320000 : 0.419  0.100  0.115  0.115

 [IOUArray]
 個数    乱数   昇順   逆順   山型
-----------------------------------
 40000 : 0.029  0.009  0.010  0.009
 80000 : 0.040  0.017  0.020  0.019
160000 : 0.090  0.035  0.043  0.040
320000 : 0.195  0.073  0.086  0.086

改良の効果は大きくて、けっこう速くなっていることがわかります。シェルソートは実装が簡単で、極端に要素数が大きくなければ十分実用になるソートです。なお、シェルソートは不安定なソートです。

●コムソート

コムソート (comb sort) はバブルソートの改良版といえる方法です。コームソートと呼ばれることもあります。今野氏の日記 (http://www.bsdclub.org/~motoyuki/d/d200005c.html#26-1) より引用します。

ソートの様子は間隔があいた櫛 (comb) ですいていくのに似ている。
  • 最初に間隔が (要素数 / 1.3) の櫛で一度すく。
  • 櫛の間隔を 1 / 1.3 しながら繰り返す。
  • 櫛の間隔が 1 になったら bubble sort
という感じである。 1.3 というのは数多くのデータを sort して実験的に得た数値とのこと。

出典は 『日経バイト 1991 年 11 月号』 とのことです。コムソートの考え方はシェルソートとよく似ています。シェルソートは挿入ソートの改良版といえる方法で、最初は遠く離れた要素間でソートを開始し、徐々に間隔を狭めていきます。そして、最後には隣り合った要素間でソートします。

シェルソートは実装が簡単で、要素数が極端に多くなければ十分実用になるソートですが、平均すればあとで説明するクイックソートよりも遅くなります。コムソートはシェルソートと同様の手法で改良を行っているようですが、それでどの程度の速さになるのかとても興味があります。

さっそくプログラムを作ってみましょう。次のリストを見てください。

リスト:コムソート

combSort :: Ord a => IOArray Int a -> IO ()
combSort buff = do
  (low, high) <- getBounds buff
  loop (high - low + 1) high
  where
    loop gap high
      | gap < 2 = bubleSort'' buff
      | otherwise = do
          let gap1 = if gap == 9 || gap == 10 then 11 else gap
          iterR_ (gap, high) (\i -> do t <- compItem buff (i - gap) i
                                       when (t == GT) $ swapItem buff (i - gap) i)
          loop (gap1 * 10 `div` 13) high

gap の間隔でデータを比較していきます。データの順序が逆ならば交換するところはバブルソートと同じです。gap を 1 / 1.3 ずつ狭めていき、gap が 1 になったらバブルソートと同じになります。

バブルソートの場合、データの交換が行われないときはソートが完了しています。櫛をすくことでデータはおおまかにソートされるため、バブルソートの途中でソートが完了する可能性は極めて高くなるでしょう。これがコムソートの原理です。

それから、間隔 gap が 9, 10 の場合には強制的に 11 にすることで、速度が向上するように改良したコムソートを comb sort 11 というそうです。今回のプログラムは comb sort 11 を使っています。

それでは実行結果を示します。

ghci> a <- newListArray (0,9) [5,6,4,7,3,8,2,9,1,0] :: IO (IOArray Int Int)
ghci> combSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> a <- newListArray (0,9) [0..9] :: IO (IOArray Int Int)
ghci> combSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
ghci> a <- newListArray (0,9) [9,8..0] :: IO (IOArray Int Int)
ghci> combSort a
ghci> getElems a
[0,1,2,3,4,5,6,7,8,9]
  表 : comb sort の結果 (単位 : 秒)

 [IOArray]
 個数    乱数   昇順   逆順   山型
-----------------------------------
 40000 : 0.047  0.022  0.026  0.027
 80000 : 0.103  0.047  0.058  0.061
160000 : 0.287  0.096  0.118  0.127
320000 : 0.819  0.203  0.249  0.261

 [IOUArray]
 個数    乱数   昇順   逆順   山型
-----------------------------------
 40000 : 0.029  0.019  0.021  0.022
 80000 : 0.054  0.038  0.041  0.042
160000 : 0.115  0.079  0.084  0.087
320000 : 0.245  0.169  0.180  0.186

バブルソートやシェーカーソートと比べると、コムソートの方が断然速いですね。しかしながら、シェルソート改良版よりは少し遅くなりました。コムソートは手続き型言語のほうが簡潔にプログラムを記述することができます。Haskwell でも同じようにプログラムしたほうが速くなるかもしれません。興味のある方は試してみてください。

今回はここまでです。次回はもっと高速なソートアルゴリズムを取り上げる予定です。


●プログラムリスト1

--
-- sort.hs : 配列のソート
--
--           Copyright (C) 2013-2021 Makoto Hiroi
--
import Data.Array.IO
import Data.IORef
import Control.Monad
import Data.Time
import System.Random
import System.Environment

-- 要素の比較
compItem :: Ord a => IOArray Int a -> Int -> Int -> IO Ordering
compItem buff i j = liftM2 (compare) (readArray buff i) (readArray buff j)

-- 要素の交換
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

-- イテレータ
iter_ :: (Int, Int) -> (Int -> IO ()) -> IO ()
iter_ (low, high) fn
  | low > high = return ()
  | otherwise  = do
      fn low
      iter_ (low + 1, high) fn

iterR_ :: (Int, Int) -> (Int -> IO ()) -> IO ()
iterR_ (low, high) fn
  | low > high = return ()
  | otherwise  = do
      fn high
      iterR_ (low, high - 1) fn

-- バブルソート
bubleSort :: Ord a => IOArray Int a -> IO ()
bubleSort buff = do
  (low, high) <- getBounds buff
  iter_
    (low, high - 1)
    (\i -> iterR_
             (i + 1, high)
             (\j -> do t <- compItem buff j (j - 1)
                       when (t == LT) $ swapItem buff j (j - 1)))

bubleSort' :: Ord a => IOArray Int a -> IO ()
bubleSort' buff = do
  (low, high) <- getBounds buff
  loop low high
  where loop low high
          | low >= high = return ()
          | otherwise = do
              f <- newIORef False
              iterR_
                (low + 1, high)
                (\j -> do t <- compItem buff j (j - 1)
                          when (t == LT) $ do 
                            swapItem buff j (j - 1)
                            writeIORef f True)
              x <- readIORef f
              when x $ loop (low + 1) high

bubleSort'' :: Ord a => IOArray Int a -> IO ()
bubleSort'' buff = do
  (low, high) <- getBounds buff
  loop low high
  where loop low high
          | low >= high = return ()
          | otherwise = do
              t <- loop1 (low + 1) high False
              when t $ loop (low + 1) high
        loop1 i j done
          | i > j     = return done
          | otherwise = do
              t <- compItem buff j (j - 1)
              if t == LT
                then do
                  swapItem buff j (j - 1)
                  loop1 i (j - 1) True
                else loop1 i (j - 1) done

-- 選択ソート

-- 最小値の選択
selectMin buff low high = iter buff (low + 1) low
  where iter buff i n
          | i > high  = return n
          | otherwise = do
              r <- compItem buff i n
              if r == LT
                then iter buff (i + 1) i
                else iter buff (i + 1) n

selectSort buff = do
  (low, high) <- getBounds buff
  iter_ (low, high - 1) (\i -> do x <- selectMin buff i high
                                  swapItem buff i x)

-- 単純挿入ソート
search_move :: Ord a => IOArray Int a -> a -> Int -> Int -> IO Int
search_move buff x i g
  | i < 0 = return (i + g)
  | otherwise = do
      y <- readArray buff i
      if x < y
        then do writeArray buff (i + g) y
                search_move buff x (i - g) g
        else return (i + g)

insertElement :: Ord a => IOArray Int a -> Int -> Int -> IO ()
insertElement buff i gap = do
  tmp <- readArray buff i
  pos <- search_move buff tmp (i - gap) gap
  writeArray buff pos tmp

insertSort :: Ord a => IOArray Int a -> IO ()
insertSort buff = do
  (low, high) <- getBounds buff
  iter_ (low + 1, high) (\i -> insertElement buff i 1)

-- シェルソート
shellSort :: Ord a => IOArray Int a -> IO ()
shellSort buff = do
  (low, high) <- getBounds buff
  loop1 ((high - low + 1) `div` 2) high
  where loop1 0   _    = return ()
        loop1 gap high = do iter_ (gap, high) (\i -> insertElement buff i gap)
                            loop1 (gap `div` 2) high

shellSort' :: Ord a => IOArray Int a -> IO ()
shellSort' buff = do
  (low, high) <- getBounds buff
  loop1 (getGap (high - low + 1) 1) high
  where getGap size gap
          | gap < size `div` 9 = getGap size (gap * 3 + 1)
          | otherwise          = gap
        loop1 0   _    = return ()
        loop1 gap high = do iter_ (gap, high) (\i -> insertElement buff i gap)
                            loop1 (gap `div` 3) high

-- コムソート
combSort :: Ord a => IOArray Int a -> IO ()
combSort buff = do
  (low, high) <- getBounds buff
  loop (high - low + 1) high
  where
    loop gap high
      | gap < 2 = bubleSort' buff
      | otherwise = do
          let gap1 = if gap == 9 || gap == 10 then 11 else gap
          iterR_ (gap, high) (\i -> do t <- compItem buff (i - gap) i
                                       when (t == GT) $ swapItem buff (i - gap) i)
          loop (gap1 * 10 `div` 13) high

-- 時間計測
test :: (IOArray Int Int -> IO ()) -> Int -> IO ()
test sort n = do
  let m = n `div` 2
      check i ary 
        | i == n    = return ()
        | otherwise = do
            t <- compItem ary (i - 1) i
            if t == GT then error "test error"
            else check (i + 1) ary
  a <- newListArray (0, n - 1) (take n (randoms (mkStdGen 11) :: [Int])) :: IO (IOArray Int Int)
  b <- newListArray (0, n - 1) [1..n] :: IO (IOArray Int Int)
  c <- newListArray (0, n - 1) [n,n-1..1] :: IO (IOArray Int Int)
  d <- newListArray (0, n - 1) ([1..m] ++ [m,m-1..1]) :: IO (IOArray Int Int)
  x1 <- getCurrentTime
  sort a
  x2 <- getCurrentTime
  check 1 a
  print (diffUTCTime x2 x1)
  x3 <- getCurrentTime
  sort b
  x4 <- getCurrentTime
  check 1 b
  print (diffUTCTime x4 x3)
  x5 <- getCurrentTime
  sort c
  x6 <- getCurrentTime
  check 1 c
  print (diffUTCTime x6 x5)
  x7 <- getCurrentTime
  sort d
  x8 <- getCurrentTime
  check 1 d
  print (diffUTCTime x8 x7)

main :: IO ()
main = do
  let xs1 = [2000, 4000, 8000, 16000]
      xs2 = [40000, 80000, 160000, 320000]
  (x:_) <- getArgs
  case x of
    "bubleSort"   -> mapM_ (test bubleSort) xs1
    "bubleSort'"  -> mapM_ (test bubleSort') xs1
    "bubleSort''" -> mapM_ (test bubleSort'') xs1
    "selectSort"  -> mapM_ (test selectSort) xs1
    "insertSort"  -> mapM_ (test insertSort) xs1
    "shellSort"   -> mapM_ (test shellSort) xs2
    "shellSort'"  -> mapM_ (test shellSort') xs2
    "combSort"    -> mapM_ (test combSort) xs2

●プログラムリスト2

--
-- sortu.hs : 配列のソート (unboxed type)
--
--            Copyright (C) 2013-2021 Makoto Hiroi
--
import Data.Array.IO
import Data.IORef
import Control.Monad
import Data.Time
import System.Random
import System.Environment

-- 要素の比較
compItem :: IOUArray Int Int -> Int -> Int -> IO Ordering
compItem buff i j = liftM2 (compare) (readArray buff i) (readArray buff j)

-- 要素の交換
swapItem :: IOUArray Int Int -> Int -> Int -> IO ()
swapItem buff i j = do
  a <- readArray buff i
  b <- readArray buff j
  writeArray buff i b
  writeArray buff j a

-- イテレータ
iter_ :: (Int, Int) -> (Int -> IO ()) -> IO ()
iter_ (low, high) fn
  | low > high = return ()
  | otherwise  = do
      fn low
      iter_ (low + 1, high) fn

iterR_ :: (Int, Int) -> (Int -> IO ()) -> IO ()
iterR_ (low, high) fn
  | low > high = return ()
  | otherwise  = do
      fn high
      iterR_ (low, high - 1) fn

-- バブルソート
bubleSort :: IOUArray Int Int -> IO ()
bubleSort buff = do
  (low, high) <- getBounds buff
  iter_
    (low, high - 1)
    (\i -> iterR_
             (i + 1, high)
             (\j -> do t <- compItem buff j (j - 1)
                       when (t == LT) $ swapItem buff j (j - 1)))

bubleSort' :: IOUArray Int Int -> IO ()
bubleSort' buff = do
  (low, high) <- getBounds buff
  loop low high
  where loop low high
          | low >= high = return ()
          | otherwise = do
              f <- newIORef False
              iterR_
                (low + 1, high)
                (\j -> do t <- compItem buff j (j - 1)
                          when (t == LT) $ do 
                            swapItem buff j (j - 1)
                            writeIORef f True)
              x <- readIORef f
              when x $ loop (low + 1) high

bubleSort'' :: IOUArray Int Int -> IO ()
bubleSort'' buff = do
  (low, high) <- getBounds buff
  loop low high
  where loop low high
          | low >= high = return ()
          | otherwise = do
              t <- loop1 (low + 1) high False
              when t $ loop (low + 1) high
        loop1 i j done
          | i > j     = return done
          | otherwise = do
              t <- compItem buff j (j - 1)
              if t == LT
                then do
                  swapItem buff j (j - 1)
                  loop1 i (j - 1) True
                else loop1 i (j - 1) done


-- 選択ソート

-- 最小値の選択
selectMin buff low high = iter buff (low + 1) low
  where iter buff i n
          | i > high  = return n
          | otherwise = do
              r <- compItem buff i n
              if r == LT
                then iter buff (i + 1) i
                else iter buff (i + 1) n

selectSort :: IOUArray Int Int -> IO ()
selectSort buff = do
  (low, high) <- getBounds buff
  iter_ (low, high - 1) (\i -> do x <- selectMin buff i high
                                  swapItem buff i x)


-- 単純挿入ソート
search_move :: IOUArray Int Int -> Int -> Int -> Int -> IO Int
search_move buff x i g
  | i < 0 = return (i + g)
  | otherwise = do
      y <- readArray buff i
      if x < y
        then do writeArray buff (i + g) y
                search_move buff x (i - g) g
        else return (i + g)

insertElement :: IOUArray Int Int -> Int -> Int -> IO ()
insertElement buff i gap = do
  tmp <- readArray buff i
  pos <- search_move buff tmp (i - gap) gap
  writeArray buff pos tmp

insertSort :: IOUArray Int Int -> IO ()
insertSort buff = do
  (low, high) <- getBounds buff
  iter_ (low + 1, high) (\i -> insertElement buff i 1)

-- シェルソート
shellSort :: IOUArray Int Int -> IO ()
shellSort buff = do
  (low, high) <- getBounds buff
  loop1 ((high - low + 1) `div` 2) high
  where loop1 0   _    = return ()
        loop1 gap high = do iter_ (gap, high) (\i -> insertElement buff i gap)
                            loop1 (gap `div` 2) high

shellSort' :: IOUArray Int Int -> IO ()
shellSort' buff = do
  (low, high) <- getBounds buff
  loop1 (getGap (high - low + 1) 1) high
  where getGap size gap
          | gap < size `div` 9 = getGap size (gap * 3 + 1)
          | otherwise          = gap
        loop1 0   _    = return ()
        loop1 gap high = do iter_ (gap, high) (\i -> insertElement buff i gap)
                            loop1 (gap `div` 3) high


-- コムソート
combSort :: IOUArray Int Int -> IO ()
combSort buff = do
  (low, high) <- getBounds buff
  loop (high - low + 1) high
  where
    loop gap high
      | gap < 2 = bubleSort'' buff
      | otherwise = do
          let gap1 = if gap == 9 || gap == 10 then 11 else gap
          iterR_ (gap, high) (\i -> do t <- compItem buff (i - gap) i
                                       when (t == GT) $ swapItem buff (i - gap) i)
          loop (gap1 * 10 `div` 13) high


test :: (IOUArray Int Int -> IO ()) -> Int -> IO ()
test sort n = do
  let m = n `div` 2
      check i ary 
        | i == n    = return ()
        | otherwise = do
            t <- compItem ary (i - 1) i
            if t == GT then error "test error"
            else check (i + 1) ary
  a <- newListArray (0, n - 1) (take n (randoms (mkStdGen 11) :: [Int])) :: IO (IOUArray Int Int)
  b <- newListArray (0, n - 1) [1..n] :: IO (IOUArray Int Int)
  c <- newListArray (0, n - 1) [n,n-1..1] :: IO (IOUArray Int Int)
  d <- newListArray (0, n - 1) ([1..m] ++ [m,m-1..1]) :: IO (IOUArray Int Int)
  x1 <- getCurrentTime
  sort a
  x2 <- getCurrentTime
  check 1 a
  print (diffUTCTime x2 x1)
  x3 <- getCurrentTime
  sort b
  x4 <- getCurrentTime
  check 1 b
  print (diffUTCTime x4 x3)
  x5 <- getCurrentTime
  sort c
  x6 <- getCurrentTime
  check 1 c
  print (diffUTCTime x6 x5)
  x7 <- getCurrentTime
  sort d
  x8 <- getCurrentTime
  check 1 d
  print (diffUTCTime x8 x7)

main :: IO ()
main = do
  let xs1 = [2000, 4000, 8000, 16000]
      xs2 = [40000, 80000, 160000, 320000]
  (x:_) <- getArgs
  case x of
    "bubleSort"   -> mapM_ (test bubleSort) xs1
    "bubleSort'"  -> mapM_ (test bubleSort') xs1
    "bubleSort''" -> mapM_ (test bubleSort'') xs1
    "selectSort"  -> mapM_ (test selectSort) xs1
    "insertSort"  -> mapM_ (test insertSort) xs1
    "shellSort"   -> mapM_ (test shellSort) xs2
    "shellSort'"  -> mapM_ (test shellSort') xs2
    "combSort"    -> mapM_ (test combSort) xs2

初版 2013 年 5 月 19 日
改訂 2021 年 8 月 8 日