今回は「乱数 (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 を調べると次のように表示されます。
ghci> :m + System.Random
ghci> :i Random
type Random :: * -> Constraint
class Random a where
randomR :: RandomGen g => (a, a) -> g -> (a, g)
default randomR :: (RandomGen g, UniformRange a) =>
(a, a) -> g -> (a, g)
random :: RandomGen g => g -> (a, g)
default random :: (RandomGen g, Uniform a) => g -> (a, g)
randomRs :: RandomGen g => (a, a) -> g -> [a]
randoms :: RandomGen g => g -> [a]
-- Defined in ‘System.Random’
instance (Random a, Random b) => Random (a, b)
-- Defined in ‘System.Random’
instance (Random a, Random b, Random c) => Random (a, b, c)
-- Defined in ‘System.Random’
instance (Random a, Random b, Random c, Random d) =>
Random (a, b, c, d)
-- Defined in ‘System.Random’
instance (Random a, Random b, Random c, Random d, Random e) =>
Random (a, b, c, d, e)
-- Defined in ‘System.Random’
instance (Random a, Random b, Random c, Random d, Random e,
Random f) =>
Random (a, b, c, d, e, f)
-- Defined in ‘System.Random’
instance (Random a, Random b, Random c, Random d, Random e,
Random f, Random g) =>
Random (a, b, c, d, e, f, g)
-- Defined in ‘System.Random’
instance Random Bool -- Defined in ‘System.Random’
instance Random Char -- Defined in ‘System.Random’
instance Random Double -- Defined in ‘System.Random’
instance Random Float -- Defined in ‘System.Random’
instance Random Int -- Defined in ‘System.Random’
instance Random Integer -- Defined in ‘System.Random’
instance Random Word -- Defined in ‘System.Random’
RandomGen は「乱数生成器」を表す型クラスです。
ghci> :i RandomGen
type RandomGen :: * -> Constraint
class RandomGen g where
next :: g -> (Int, g)
genWord8 :: g -> (GHC.Word.Word8, g)
genWord16 :: g -> (GHC.Word.Word16, g)
genWord32 :: g -> (GHC.Word.Word32, g)
genWord64 :: g -> (GHC.Word.Word64, g)
genWord32R :: GHC.Word.Word32 -> g -> (GHC.Word.Word32, g)
genWord64R :: GHC.Word.Word64 -> g -> (GHC.Word.Word64, g)
genShortByteString :: Int
-> g -> (Data.ByteString.Short.Internal.ShortByteString, g)
genRange :: g -> (Int, Int)
split :: g -> (g, g)
{-# MINIMAL split, (genWord32 | genWord64 | next, genRange) #-}
-- Defined in ‘System.Random.Internal’
instance RandomGen StdGen -- Defined in ‘System.Random.Internal’
RandomGen は整数の一様乱数を生成します。next と genRange は random-1.2 で非推奨 (Deprecated) になりました。そのかわり、無符号整数 (Word) の一様乱数を生成する関数 genWord* が用意されています。今までのように、型クラス Random に定義されている関数を使うこともできます。
Haskell で乱数を生成する場合、何らかの乱数生成器が必要になります。System.Random に用意されている乱数生成器が StdGen です。
ghci> :i StdGen
type StdGen :: *
newtype StdGen
= System.Random.Internal.StdGen {System.Random.Internal.unStdGen :: System.Random.SplitMix.SMGen}
-- Defined in ‘System.Random.Internal’
instance RandomGen StdGen -- Defined in ‘System.Random.Internal’
instance Eq StdGen -- Defined in ‘System.Random.Internal’
instance Show StdGen -- Defined in ‘System.Random.Internal’
StdGen を生成する関数が mkStdGen です。
mkStdGen :: Int -> StdGen
mkStdGen は整数 (Int) から StdGen を生成します。また、Haskell は起動時にシステム内の適当な値 (たとえば時刻など) を使って StdGen をひとつ生成します。この値は大域変数 getStdGen に格納されています。
getStdGen :: IO StdGen
それでは実際に試してみましょう。
ghci> g = mkStdGen 11
ghci> genWord32 g
(4020123995,StdGen {unStdGen = SMGen 4664641791676752737 5833679380957638813})
ghci> (n1, g1) = genWord32 g
ghci> n1
4020123995
ghci> g1
StdGen {unStdGen = SMGen 4664641791676752737 5833679380957638813}
ghci> (n2, g2) = genWord32 g1
ghci> n2
2965546545
ghci> g2
StdGen {unStdGen = SMGen 10498321172634391550 5833679380957638813}
ghci> g3 = mkStdGen 256
ghci> g3
StdGen {unStdGen = SMGen 15401116602503760187 7350668447369684287}
ghci> (n4, g4) = genWord32 g3
ghci> n4
1165258315
ghci> g4
StdGen {unStdGen = SMGen 4305040976163892858 7350668447369684287}
ghci> (n5, g5) = genWord32 g4
ghci> n5
3516030425
ghci> g5
StdGen {unStdGen = SMGen 11655709423533577145 7350668447369684287}
ghci> a <- getStdGen
ghci> a
StdGen {unStdGen = SMGen 7146799629592057685 9833312923008228965}
ghci> (m1, a1) = genWord32 a
ghci> m1
1503826549
ghci> a1
StdGen {unStdGen = SMGen 16980112552600286650 9833312923008228965}
ghci> (m2, a2) = genWord32 a1
ghci> m2
2905235083
ghci> a2
StdGen {unStdGen = SMGen 8366681401898963999 9833312923008228965}
0 以上 upper 以下の乱数を生成する場合、関数 genRand32R, genRand64R を使うと便利です。
関数が randomR は lower 以上 upper 以下の乱数を生成します。
randomR :: (RandomGen g, Random a) => (a, a) -> g -> (a, g)
randomR は生成する乱数の範囲と乱数生成器を受け取り、乱数と新しい乱数生成器を返します。同じ乱数生成器を渡すと同じ乱数しか生成しません。返された乱数生成器を使って新しい乱数を生成します。ご注意くださいませ。
簡単な実行例を示します。
ghci> (n1, g1) = randomR (1, 6) (mkStdGen 11)
ghci> n1
4
ghci> g1
StdGen {unStdGen = SMGen 4664641791676752737 5833679380957638813}
ghci> (n2, g2) = randomR (1, 6) g1
ghci> n2
2
ghci> g2
StdGen {unStdGen = SMGen 10498321172634391550 5833679380957638813}
ghci> (n3, g3) = randomR (1, 6) g2
ghci> n3
1
ghci> g3
StdGen {unStdGen = SMGen 3718935860840117560 5833679380957638813}
関数 randomRs は乱数をリストに格納して返します。返り値は遅延リスト (無限リスト) になることに注意してください。
randomRs :: (RandomGen g, Random a) => (a, a) -> g -> [a]
簡単な実行例を示します。
ghci> xs = randomRs (1, 6) (mkStdGen 11)
ghci> take 16 xs
[4,2,1,5,6,3,3,4,5,6,5,3,4,5,1,2]
ghci> ys = randomRs ('a', 'z') (mkStdGen 11)
ghci> take 26 ys
"yrbjxzjklmtfisiyptroofihvv"
関数 random, randoms は randomR, randomRs と似ていますが、生成する乱数の範囲の指定がありません。
random :: (RandomGen g, Random a) => g -> (a, g) randoms :: (RandomGen g, Random a) => g -> [a]
乱数の範囲は指定したデータ型の (最小値, 最大値) となります。Integer の場合は Int と同じになります。実数の場合は 0 以上 1.0 未満になります。
簡単な実行例を示します。
ghci> (n1, g1) = random (mkStdGen 11) :: (Integer, StdGen)
ghci> n1
-8691282579076269733
ghci> g1
StdGen {unStdGen = SMGen 4664641791676752737 5833679380957638813}
ghci> (n2, g2) = random g1 :: (Integer, StdGen)
ghci> n2
-6401397088113744335
ghci> g2
StdGen {unStdGen = SMGen 10498321172634391550 5833679380957638813}
ghci> (n3, g3) = random g2 :: (Integer, StdGen)
ghci> n3
-7708536135073118066
ghci> g3
StdGen {unStdGen = SMGen 16332000553592030363 5833679380957638813}
ghci> xs = randoms (mkStdGen 11) :: [Integer]
ghci> take 10 xs
[-8691282579076269733,-6401397088113744335,-7708536135073118066,1428399377092000640,
-6279174627358091186,-203332483515836804,-9116495703984563675,-976259336053638646,
3305629617079671367,2460966432541861862]
ghci> ys = randoms (mkStdGen 11) :: [Double]
ghci> take 5 ys
[0.4711553726959955,0.34702043149376516,0.41788058121646443,0.9225663146089955,
0.340394738619875]
大域変数 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 アクションになります。
簡単な実行例を示します。
ghci> rand = getStdRandom random :: IO Integer ghci> rand 8359930581553067021 ghci> rand 4962155463747272277 ghci> rand 910859405987817638 ghci> rand 4249964703329071031 ghci> dice = getStdRandom (randomR (1,6)) ghci> dice 4 ghci> dice 4 ghci> dice 3 ghci> dice 3 ghci> dice 1 ghci> dice 5 ghci> dice 3 ghci> randomIO :: IO Integer 3439421987660092452 ghci> randomIO :: IO Integer 5035637377291934312 ghci> randomIO :: IO Integer 4526479393393332296 ghci> randomIO :: IO Integer -5633929731840363662 ghci> randomIO :: IO Integer 2656556158495842224 ghci> randomIO :: IO Integer 7653980111264477478 ghci> randomRIO (1, 6) 6 ghci> randomRIO (1, 6) 3 ghci> randomRIO (1, 6) 3 ghci> randomRIO (0.0, 1.0) 0.6041613705219417 ghci> randomRIO (0.0, 1.0) 0.9553834417651247 ghci> 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
簡単なプログラムなので説明は省略します。リストを読んでくださいね。それでは実行結果を示します。
ghci> test 1000000 3.142244 ghci> test 1000000 3.142428 ghci> test 1000000 3.142156 ghci> test 1000000 3.14072 ghci> test 1000000 3.143336
このように、乱数を使って \(\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 buff 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.025 0.005 0.004 0.008 160000 : 0.031 0.008 0.007 0.017 320000 : 0.061 0.015 0.017 0.035 640000 : 0.127 0.031 0.033 0.073 実行環境 : Ubunts 22.04 (WSL2), Intel Core i5-6200U 2.30GHz
乱数を生成する分だけ、実行速度はどうしても遅くなってしまいます。ですが、山型データでも高速にソートすることができました。このように、乱数を使うことにより、簡単な方法で最悪のケースを回避することができます。
最後に「線形合同法」という簡単な方法を使って乱数を生成するプログラムを作ってみましょう。線形合同法は適当な初期値からスタートして、次の漸化式で乱数を生成します。
これで 0 以上 M 未満の値を生成することができます。A, C, M の選び方は参考文献『Ｃ言語による最新アルゴリズム事典』より引用します。
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 の乱数を生成します。
ghci> rand g = (g, (13 * g + 1) `mod` 16) ghci> rands g = let (n, g') = rand g in n : rands g' ghci> xs = rands 1 ghci> take 16 xs [1,14,7,12,13,10,3,8,9,6,15,4,5,2,11,0] ghci> 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 は乱数をリストに格納して返します。遅延リスト (無限リスト) になることに注意してください。
整数の次は実数です。参考文献『Ｃ言語による最新アルゴリズム事典』によると、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 は実数の乱数をリストに格納して返します。
実行結果は次のようになります。
ghci> xs = rands (mkRandGen 1) ghci> take 10 xs [1,69070,475628535,3277404108,772999773,3877832058,3821835443,1662200408,2044158073, 3788989926] ghci> ys = rands (mkRandGen 123456789) ghci> take 10 ys [123456789,1526876882,1232376875,1376507248,583047857,919067838,3774835239,2400386108, 2235500557,4008647530] ghci> xs1 = randoms (mkRandGen 1) ghci> 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] ghci> ys1 = randoms (mkRandGen 123456789) ghci> 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 から説明します。線形合同法の場合、下位 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] をお読みください。