乱数
今回は「乱数 (random numbers)」について簡単に説明します。
●乱数とは？
私たちが適当な数字を決める場合、サイコロを使えば 1 から 6 までの数字を簡単に決めることができます。たとえば、サイコロを振って出た目を記録したら次のようになったとしましょう。
5, 2, 1, 2, 6, 3, 4, 3, 1, 5, .....
サイコロの目はイカサマをしないかぎり、出る確率が 1 / 6 で規則性はまったくありません。したがって、数字の出る順番には規則性はなく、まったくでたらめになります。いま 2 が出たから次は 1 が出るとか、3, 4 と続いたから次は 5 が出るなどのように、前に出た数字から次に出る数字を予測することはできません。このように、でたらめに並んだ数列を「乱数列 (random sequence)」といい、乱数列の中のひとつひとつの数字を乱数 (random numbers) といいます。
コンピュータは決められた手順（プログラム）を高速に実行することは得意なのですが、まったくでたらめの数を作れといわれると、とたんに困ってしまいます。そこで、何かしらの数式をプログラムして、それを実行することで乱数を発生させます。厳密にいえば乱数ではありませんが、それを乱数としてみなして使うことにするのです。このような乱数を「疑似乱数 (pseudo-random numbers)」といいます。
●擬似乱数の性質
乱数には重要な性質が 2 つあります。ひとつは「等確率性 (等出現性)」です。たとえば、サイコロの目はイカサマをしないかぎり、出る確率が 1/6 になります。したがって、サイコロを振る回数 (試行回数 n) を増やしていくと、どの目も出現する回数は理論値 (n / 6) に限りなく近づいていきます。このような確率的な性質を「大数の法則」といいます。そして、ある範囲の整数値が等確率で現れる乱数を「整数の一様乱数」といいます。実数の一様乱数は、整数の一様乱数から生成することができます。
もうひとつの重要な性質が「無規則性 (独立性)」です。たとえば、サイコロには規則性がないので、目の出る順番はまったくでたらめになります。いま 2 が出たから次は 1 が出るとか、3, 4 と続いたから次は 5 が出るなどのような規則性はありません。ようするに、以前に出現した数によって次に出現する数が決まるようなことはないのです。これが無規則性で、出現した各数には相関性がないことから、無相関性とか独立性と呼ばれることもあります。
ただし、コンピュータで生成する擬似乱数の場合、無規則性を完全に実現することはできない、ということに注意してください。一般的な擬似乱数のアルゴリズムは、コンピュータのメモリに現在の状態 (数値) を記憶しておきます。乱数を生成するとき、メモリの内容をある手順 (A) で書き換えて新しい状態へ移行します。そして、その状態をある手順 (B) で数値に変換し、それを乱数として出力します。これを繰り返すことで乱数列を生成しています。
ここで、手順 (A) や (B) をどんなに複雑なものにしたとしても、あらかじめ「決められた手順 (アルゴリズム)」を実行していることにかわりはありません。したがって、擬似乱数の規則性をゼロにすることは原理的に不可能というわけです。また、コンピュータにあるメモリは有限なので、擬似乱数は周期性を持つこともわかるでしょう。
これに対し、何らかの物理的現象 (たとえば電子回路の熱雑音など) を利用して乱数を生成する方法があります。このような方法で生成された乱数を「物理乱数」と呼ぶ場合があります。物理乱数は規則性がまったくないので、まさに「真の乱数」といえるでしょう。そのかわり、同じ乱数列を再現するには、発生した乱数を記憶しておかない限り不可能です。つまり、物理乱数には「再現性」がまったくないのです。
擬似乱数は真の乱数ではありませんが、物理乱数にはない「再現性」があります。擬似乱数は、種 (シード) に同じ値を設定することで同じ乱数列を簡単に発生させることができます。シミュレーションやモンテカルロ法などの実験で乱数を用いる場合、この機能はとても役に立ちます。シードに別な値を設定すれば、異なる乱数列を簡単に発生させることができますし、シードに同じ値を与えれば、追試も簡単に行うことができます。
●乱数生成器
Haskell の場合、モジュール System.Random をインポートすると乱数を簡単に生成することができます。System.Random が見つからない場合は、パッケージ random をインストールしてください。stack を使うと、次のコマンドで簡単にインストールすることができます。
$ stack install random
乱数を表す型クラスが Random です。ghci の :info コマンドで Random を調べると次のように表示されます。
Prelude> :m + System.Random
Prelude System.Random> :i Random
class Random a where
randomR :: RandomGen g => (a, a) -> g -> (a, g)
random :: RandomGen g => g -> (a, g)
randomRs :: RandomGen g => (a, a) -> g -> [a]
randoms :: RandomGen g => g -> [a]
randomRIO :: (a, a) -> IO a
randomIO :: IO a
{-# MINIMAL randomR, random #-}
-- Defined in ‘System.Random’
instance Random Word -- Defined in ‘System.Random’
instance Random Integer -- Defined in ‘System.Random’
instance Random Int -- Defined in ‘System.Random’
instance Random Float -- Defined in ‘System.Random’
instance Random Double -- Defined in ‘System.Random’
instance Random Char -- Defined in ‘System.Random’
instance Random Bool -- Defined in ‘System.Random’
RandomGen は「乱数生成器」を表す型クラスです。
Prelude System.Random> :i RandomGen
class RandomGen g where
next :: g -> (Int, g)
genRange :: g -> (Int, Int)
split :: g -> (g, g)
{-# MINIMAL next, split #-}
-- Defined in ‘System.Random’
instance RandomGen StdGen -- Defined in ‘System.Random’
RandomGen は整数の一様乱数を生成します。その範囲は genRange の返り値 (Int, Int) で表します。next は乱数を生成して新しい乱数生成器といっしょにタプルに格納して返します。split は二つの新しい乱数生成器をタプルに格納して返します。なお、実際に乱数を使用するときは next を使うのではなく、型クラス Random に定義されている関数を使います。
Haskell で乱数を生成する場合、何らかの乱数生成器が必要になります。System.Random に用意されている乱数生成器が StdGen です。
Prelude System.Random> :i StdGen
data StdGen
= System.Random.StdGen {-# UNPACK #-}GHC.Int.Int32
{-# UNPACK #-}GHC.Int.Int32
-- Defined in ‘System.Random’
instance Show StdGen -- Defined in ‘System.Random’
instance RandomGen StdGen -- Defined in ‘System.Random’
instance Read StdGen -- Defined in ‘System.Random’
StdGen を生成する関数が mkStdGen です。
mkStdGen :: Int -> StdGen
mkStdGen は整数 (Int) から StdGen を生成します。また、Haskell は起動時にシステム内の適当な値 (たとえば時刻など) を使って StdGen をひとつ生成します。この値は大域変数 getStdGen に格納されています。
getStdGen :: IO StdGen
それでは実際に試してみましょう。
Prelude System.Random> g = mkStdGen 11
Prelude System.Random> g
12 1
Prelude System.Random> genRange g
(0,2147483562)
Prelude System.Random> (n1, g1) = next g
Prelude System.Random> n1
439476
Prelude System.Random> g1
480168 40692
Prelude System.Random> next g1
(377734984,2033573848 1655838864)
Prelude System.Random> mkStdGen 256
257 1
Prelude System.Random> mkStdGen 123456789
123456790 1
Prelude System.Random> a <- getStdGen
Prelude System.Random> a
2094583970 1
このように、StdGen は二つの整数 (Int32) で構成されています。
●乱数の取得
乱数生成器から乱数を取得する基本的な関数が randomR です。
randomR :: (RandomGen g, Random a) => (a, a) -> g -> (a, g)
randomR は生成する乱数の範囲と乱数生成器を受け取り、乱数と新しい乱数生成器を返します。同じ乱数生成器を渡すと同じ乱数しか生成しません。返された乱数生成器を使って新しい乱数を生成します。ご注意くださいませ。
簡単な実行例を示します。
Prelude System.Random> (n1, g1) = randomR (1, 6) (mkStdGen 11)
Prelude System.Random> n1
6
Prelude System.Random> g1
480168 40692
Prelude System.Random> (n2, g2) = randomR (1, 6) g1
Prelude System.Random> n2
4
Prelude System.Random> g2
2033573848 1655838864
Prelude System.Random> (n3, g3) = randomR (1, 6) g2
Prelude System.Random> n3
4
Prelude System.Random> g3
1124268239 2103410263
関数 randomRs は乱数をリストに格納して返します。返り値は遅延リスト (無限リスト) になることに注意してください。
randomRs :: (RandomGen g, Random a) => (a, a) -> g -> [a]
簡単な実行例を示します。
Prelude System.Random> xs = randomRs (1, 6) (mkStdGen 11)
Prelude System.Random> take 16 xs
[6,4,4,6,1,2,6,5,6,5,3,2,2,2,4,4]
Prelude System.Random> ys = randomRs ('a', 'z') (mkStdGen 11)
Prelude System.Random> take 26 ys
"xpzrspdkfuuzlhfxmovwuamjkt"
関数 random, randoms は randomR, randomRs と似ていますが、生成する乱数の範囲の指定がありません。
random :: (RandomGen g, Random a) => g -> (a, g)
randoms :: (RandomGen g, Random a) => g -> [a]
乱数の範囲は指定したデータ型の (最小値, 最大値) となります。Integer の場合は Int と同じになります。実数の場合は 0 以上 1.0 未満になります。
簡単な実行例を示します。
Prelude System.Random> (n1, g1) = random (mkStdGen 11) :: (Integer, StdGen)
Prelude System.Random> n1
5260538044923710387
Prelude System.Random> g1
1124268239 2103410263
Prelude System.Random> (n2, g2) = random g1 :: (Integer, StdGen)
Prelude System.Random> n2
4361398698747678847
Prelude System.Random> g2
1908496843 1780294415
Prelude System.Random> (n3, g3) = random g2 :: (Integer, StdGen)
Prelude System.Random> n3
-8221315287270277529
Prelude System.Random> g3
857054714 1422611300
Prelude System.Random> xs = randoms (mkStdGen 11) :: [Integer]
Prelude System.Random> take 10 xs
[5260538044923710387,4361398698747678847,-8221315287270277529,7278185606566790575,1652507602255180489,
6207436798600535810,2828579254075873640,-3624220446475129215,-866363518444507197,-2767556111145761201]
Prelude System.Random> ys = randoms (mkStdGen 11) :: [Double]
Prelude System.Random> take 5 ys
[3.704593909093601e-2,0.21252549198699444,0.25064753697547537,4.092379058618245e-2,0.4651988391811638]
●大域乱数発生器から乱数を取得する
大域変数 getStdGen から乱数を取得する場合は、関数 getStdRandom, randomIO, randomRIO を使うと簡単です。
getStdRandom :: (StdGen -> (a, StdGen)) -> IO a
randomIO :: Random a => IO a
randomRIO :: Random a => (a, a) -> IO a
getStdRandom は引数に random や randomR を受け取って、返り値を変数に束縛します。このあと、その変数を参照するたびに乱数が生成されます。もっと簡単なのが randomIO と randomRIO です。randomIO は参照するたびに乱数が生成され、randomRIO は呼び出すたびに乱数が生成されます。当然ですが、これらの動作はすべて I/O アクションになります。
簡単な実行例を示します。
Prelude System.Random> rand = getStdRandom random :: IO Integer
Prelude System.Random> rand
8359930581553067021
Prelude System.Random> rand
4962155463747272277
Prelude System.Random> rand
910859405987817638
Prelude System.Random> rand
4249964703329071031
Prelude System.Random> dice = getStdRandom (randomR (1,6))
Prelude System.Random> dice
4
Prelude System.Random> dice
4
Prelude System.Random> dice
3
Prelude System.Random> dice
3
Prelude System.Random> dice
1
Prelude System.Random> dice
5
Prelude System.Random> dice
3
Prelude System.Random> randomIO :: IO Integer
3439421987660092452
Prelude System.Random> randomIO :: IO Integer
5035637377291934312
Prelude System.Random> randomIO :: IO Integer
4526479393393332296
Prelude System.Random> randomIO :: IO Integer
-5633929731840363662
Prelude System.Random> randomIO :: IO Integer
2656556158495842224
Prelude System.Random> randomIO :: IO Integer
7653980111264477478
Prelude System.Random> randomRIO (1, 6)
6
Prelude System.Random> randomRIO (1, 6)
3
Prelude System.Random> randomRIO (1, 6)
3
Prelude System.Random> randomRIO (0.0, 1.0)
0.6041613705219417
Prelude System.Random> randomRIO (0.0, 1.0)
0.9553834417651247
Prelude System.Random> randomRIO (0.0, 1.0)
0.16446221133916183
●モンテカルロ法
それでは、乱数を使った簡単な応用例を示しましょう。乱数を使って数学や物理などの問題を解くアルゴリズムを「モンテカルロ法 (Monte Carlo methods)」といいます。簡単な例として、円周率 \(\pi\) をモンテカルロ法で求めてみます。
正方形の領域 (\(0 \leq x \leq 1, \ 0 \leq y \leq 1\)) に乱数で点を打ちます。乱数であれば点は領域内に一様に分布するので、\(x^2 + y^2 \leq 1\) の円内に入る確率は \(\frac{\pi}{4}\) になります。つまり、(円内の点の個数 / 点の総数) の値は 0.7853... になるはずです。たくさん点を打つほど値は \(\frac{\pi}{4}\) に近づくはずですが、コンピュータの乱数は疑似乱数なので規則性が生じてしまい、値の精度にはどうしても限界があります。それでも、性能の良い疑似乱数ほど \(\frac{\pi}{4}\) に近い値になるでしょう。プログラムは次のようになります。
リスト : モンテカルロ法（πを求める）
import System.Random
test :: Int -> IO Double
test n = iter n 0.0
where iter 0 a = return (a * 4 / (fromIntegral n))
iter m a = do
x <- randomRIO (0.0, 1.0) :: IO Double
y <- randomRIO (0.0, 1.0)
if x * x + y * y <= 1.0
then iter (m - 1) (a + 1)
else iter (m - 1) a
簡単なプログラムなので説明は省略します。リストを読んでくださいね。それでは実行結果を示します。
*Main> test 1000000
3.140892
*Main> test 1000000
3.141076
*Main> test 1000000
3.143876
*Main> test 1000000
3.140868
*Main> test 1000000
3.1409
このように、乱数を使って \(\pi\)の近似値を求めることができます。
●乱数とクイックソート
次は、乱数を使って最悪のケースを回避する方法を紹介します。拙作のページ 配列のソート (2) で説明したクイックソートは、枢軸 (pivot) を基準にして大きなデータと小さなデータの 2 つの区間に分けてソートを行います。このとき、pivot の選び方でクイックソートの効率は大きく左右されます。最悪のケースが最小値または最大値を選ぶ場合で、クイックソートの実行時間はデータ数の 2 乗に比例する遅いソートになってしまいます。
このため、ソートする区間から数個のデータを選び、その中から中央の値を枢軸に選ぶ方法があります。この方法は 配列のソート (2) で詳しく説明しました。このほかに、もう一つ方法があります。それが乱数で枢軸を選ぶ方法です。次のリストを見てください。
リスト : クイックソート (乱択アルゴリズム)
search1 :: IOUArray Int Int -> Int -> 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 :: IOUArray Int Int -> Int -> 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 :: IOUArray Int Int -> Int -> 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)
quickSortR :: IOUArray Int Int -> IO ()
quickSortR buff = do
(low, high) <- getBounds buff
qsort low high
where
selectPv low high = do
a <- randomRIO (low, high)
x <- readArray buff a
return x
qsort low high = do
pv <- selectPv low high
(i, j) <- quickPartition pv low high
when (low < i - 1) $ qsort low (i - 1)
when (high > j + 1) $ qsort (j + 1) high
区間 low, high の中から枢軸 pivot を選ぶ関数 selectPv の中で randomRIO を使っています。区間の中央を枢軸にすると、山型データが最悪のケースになり、要素数の二乗に比例する遅いソートになってしまいます。この場合、枢軸をランダムで選ぶことにより、最悪のケースを回避することができます。もちろん、完璧に回避できるわけではありませんが、最悪のケースを何回も続けて選ぶ確率は相当に低くなると思われます。山型データでも高速にソートできるでしょう。
実行結果は次のようになりました。
表 : quickSortR の結果 (単位 : 秒)
[IOUArray]
個数 乱数 昇順 逆順 山型
-----------------------------------
80000 : 0.033 0.018 0.017 0.030
160000 : 0.068 0.035 0.046 0.061
320000 : 0.137 0.078 0.081 0.124
640000 : 0.287 0.153 0.156 0.243
実行環境 : Ubunts 18.04 (WSL), Intel Core i5-6200U 2.30GHz
乱数を生成する分だけ、実行速度はどうしても遅くなってしまいます。ですが、山型データでも高速にソートすることができました。このように、乱数を使うことにより、簡単な方法で最悪のケースを回避することができます。
●線形合同法
最後に「線形合同法」という簡単な方法を使って乱数を生成するプログラムを作ってみましょう。線形合同法は適当な初期値からスタートして、次の漸化式で乱数を生成します。
\(
X_i \equiv (A \times X_{i-1} + C) \pmod M
\)
これで 0 以上 M 未満の値を生成することができます。A, C, M の選び方は 参考文献 1 より引用します。
M が 2 の累乗なら A mod 8 を 5 または 3 とし（5 の方が安全）、定数項 C を奇数とする。このとき、周期はちょうど M となり、その 1 周期分には 0 から M - 1 までの整数が 1 個ずつ現れる。
高速にするには C = 0 とし、初期値 \(X_0\) を奇数にする。ただし、周期は M / 4 になる。
C = 0 の場合を「乗算合同法」といいます。初期値のことを「種 (seed)」と呼ぶ場合があります。
それでは実際にプログラムを作って確かめてみましょう。簡単な例として、M を 16 にして 0 - 15 の乱数を生成します。
Prelude> rand g = (g, (13 * g + 1) `mod` 16)
Prelude> rands g = let (n, g') = rand g in n : rands g'
Prelude> xs = rands 1
Prelude> take 16 xs
[1,14,7,12,13,10,3,8,9,6,15,4,5,2,11,0]
Prelude> take 32 xs
[1,14,7,12,13,10,3,8,9,6,15,4,5,2,11,0,1,14,7,12,13,10,3,8,9,6,15,4,5,2,11,0]
関数 rand は 1 周期が 16 で、0 から 15 までの値を 1 個ずつ生成していることがわかります。
実際には M の値を \(2^{32}\) のように大きくして周期を長くします。次のリストを見てください。
リスト : 線形合同法
module Rand where
import Data.Word
-- 乱数生成器の定義
newtype RandGen = RandGen {seed :: Word32} deriving Show
mkRandGen :: Int -> RandGen
mkRandGen n = RandGen (fromIntegral n)
-- 整数の一様乱数
rand :: RandGen -> (Integer, RandGen)
rand g = (fromIntegral (seed g), RandGen (69069 * seed g + 1))
rands :: RandGen -> [Integer]
rands g =
let (n, g') = rand g
in n : rands g'
-- 実数の一様乱数
random :: RandGen -> (Double, RandGen)
random g =
let (n, g') = rand g
in (1.0 / 4294967296.0 * (0.5 + fromIntegral n), g')
randoms :: RandGen -> [Double]
randoms g =
let (n, g') = random g
in n : randoms g'
newtype で乱数生成器 RandGen を定義します。RandGen には無符号 32 ビット整数 Word32 を格納します。これで\(\mod 2^{32}\) の計算を省くことができます。関数 mkRandGen は整数から RandGen を作ります。fromIntegral で Int を Word32 に変換するだけです。
関数 rand は RandGen に格納されている値を乱数として返します。生成される乱数の範囲は 0 以上 \(2^{32}\) (4294967296) 未満になります。乱数のデータ型は Integer としました。このとき、次の乱数を計算して新しい RandGen を生成し、乱数と RandGen をタプルに格納して返します。関数 rands は乱数をリストに格納して返します。遅延リスト (無限リスト) になることに注意してください。
整数の次は実数です。参考文献 1 によると、0 以上 1 未満の範囲の実数の一様乱数 U は、整数の一様乱数を生成する関数 rand を使って次式のように作ることができるそうです。
(1) U = (1.0 / (RAND_MAX + 1.0)) * (rand() + 0.5)
(2) U = (1.0 / (RAND_MAX + 1.0)) * rand()
RAND_MAX は関数 rand() が生成する乱数の最大値を表します。式 (1) は 『正確に 0 や 1 になることがないので、対数 log( U ) を求める際にも安心である。』 とのことですが、通常は式 (2) で十分だそうです。関数 random は式 (1) で 0 以上 1 未満の実数を生成します。関数 randoms は実数の乱数をリストに格納して返します。
実行結果は次のようになります。
*Rand> xs = rands (mkRandGen 1)
*Rand> take 10 xs
[1,69070,475628535,3277404108,772999773,3877832058,3821835443,1662200408,2044158073,3788989926]
*Rand> ys = rands (mkRandGen 123456789)
*Rand> take 10 ys
[123456789,1526876882,1232376875,1376507248,583047857,919067838,3774835239,2400386108,2235500557,4008647530]
*Rand> xs1 = randoms (mkRandGen 1)
*Rand> take 10 xs1
[3.4924596548080444e-10,1.6081728972494602e-5,0.11074089806061238,0.7630801080958918,
0.17997803480830044,0.9028781341621652,0.8898404062492773,0.38701119099278003,
0.47594263998325914,0.8821929633850232]
*Rand> ys1 = randoms (mkRandGen 123456789)
*Rand> take 10 ys1
[2.8744523762725294e-2,0.35550372733268887,0.28693510114680976,0.32049306866247207,
0.13575140794273466,0.2139871563995257,0.8788973185000941,0.5588834426598623,
0.5204930336913094,0.93333598470781]
●線形合同法の欠点
このように、線形合同法はとても簡単なアルゴリズムなのですが、使用上の注意点がいくつかあります。参考文献 1 によると、線形合同法には一般に次の欠点があるそうです。
- 合同法乱数は、上位の桁はランダムだが、下位の桁はランダムでない
- 合同法乱数は、1 個ずつ使えばランダムだが、いくつか組にして使えばランダムでない
まず、1 から説明します。線形合同法の場合、下位 k bit だけを見ると周期は高々 2 k になります。たとえば、線形合同法で生成した乱数 N から、0 以上 8 未満の乱数を作ってみましょう。この場合、N mod 8 として乱数を生成すると周期が 8 にしかなりません。下位 3 bit ではなく、上位 3 bit を使った方が良いようです。
それでは、実際に確かめてみましょう。線形合同法のプログラムで、上位 3 bit と下位 3 bit の値を出力します。結果は次のようになりました。
rand L H (L:下位 3 bit, H:上位 3 bit)
------------------
1, 1, 0
69070, 6, 0
475628535, 7, 0
3277404108, 4, 6
772999773, 5, 1
3877832058, 2, 7
3821835443, 3, 7
1662200408, 0, 3
2044158073, 1, 3
3788989926, 6, 7
797919023, 7, 1
2743624612, 4, 5
1156259413, 5, 2
1059494674, 2, 1
584849259, 3, 1
786050992, 0, 1
3369345009, 1, 6
下位 3 bit で乱数を生成すると、周期が 8 にしかならないことがわかります。また、rands の結果を見ると、偶数と奇数が交互に出現しています。これも線形合同法の特徴 (欠点) です。このように、線形合同法の下位の桁には規則性があるのです。
次は、もう一つの欠点 2 について説明します。たとえば 0 から 15 までの乱数で、(x, y) の座標を生成することを考えます。ここで周期 16 の線形合同法で乱数を生成してみましょう。すると、生成される (x, y) の座標は 16 / 2 = 8 通りにしかなりません。(x, y) の座標は 256 個もあるのに、8 個の座標しか生成できないのです。(x, y, z) の場合は 4096 個の座標があるのに、16 個の座標しか生成できません。極端な例ですが、これが線形合同法の欠点になります。
線形合同法の欠点の詳細な説明と改良方法については、拙作のページ Algorithms with Python 乱数 [1] をお読みください。
●参考文献
- 奥村晴彦, 『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
初版 2013 年 6 月 2 日
改訂 2021 年 8 月 8 日
選択アルゴリズム
N 個のデータのうち小さい方から数えて k 番目の値を求める方法を「選択アルゴリズム (Selection Algorithm)」と呼びます。簡単な方法はクイックソートなどでソートすることですが、データ数 N が多くなるとクイックソートでも N * log2 N に比例する時間が必要になります。ところが、クイックソートを変形することで、データ数 N に比例する程度の時間で k 番目のデータを求めることができます。これを「クイックセレクト (quick select)」と呼びます。
●クイックセレクト
クイックソートはある値（枢軸）を基準にして、要素をそれより大きいものと小さいものの 2 つに分割していくことでソートを行います。2 つに分けた各々のグループを同様に分割して、さらに 2 つのグループに分けます。分割を続けていくと最後は要素がひとつになりソートが完了します。
クイックセレクトも枢軸を基準に二分割していきますが、枢軸の決め方と二分割したグループの扱い方が異なるのです。
上図を見てください。配列 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 上で実行したところ、実行結果は次のようになりました。
*Main> test 5000
27909750426213363
0.0418992s
27909750426213363
0.1871055s
*Main> test 10000
34734837088680882
0.0872716s
34734837088680882
0.4090038s
*Main> test 20000
-41460865972116048
0.0728619s
-41460865972116048
0.8497877s
*Main> test 40000
3380968318790221
0.279833s
3380968318790221
1.9006622s
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)
初版 2013 年 6 月 2 日
改訂 2021 年 8 月 8 日
