maximal length sequence - Rokiのチラ裏

M系列に関する学習メモ*1。

import Test.HUnit
import System.IO

mulcon :: Int -> Int
mulcon 0 = 1
mulcon n = (a * mulcon(n - 1) + b) `mod` m
    where
        a = 3
        b = 0
        m = 7

mc :: [Int] -> [Int]
mc = map mulcon

main :: IO (Counts, Int)
main = runTestText (putTextToHandle stderr False) tests
    where
        x = mc [0 .. 11]
        tests = TestList [ "mc test" ~: x ~?= [1, 3, 2, 6, 4, 5, 1, 3, 2, 6, 4, 5] ]

これは線形合同法 *2(漸化式 $\displaystyle X_{n} = ax_{n - 1} + b \pmod{m}$ )の単純な実装であるが、結果からも分かるように、周期は $m$ に強く依存するため*3、大きな値 $m$ を用いることで長い周期の乱数を得る事ができる。しかしそれに応じて計算処理の膨大化の他、既知の様々な問題*4がある。M系列は $\bigoplus$ を Bitwise XOR とした時、次の高次線形漸化式で発生する 1 ビット数列を言い、線形合同法の欠点を克服する。

$\displaystyle \large \begin{equation} x_{n} = x_{n - p} \bigoplus x_{n - q} \ \ \ \ (p \gt q) \end{equation}$

各項の値は 1 ビット、すなわち 0 か 1 で事前に $\displaystyle x_{n-1},\ x_{n-2},\ \cdots,\ x_{n-p}$ に対して初期値を与えておく。次のコードは、 $x_{0} = 1,\ x_{1} = 0,\ x_{2} = 0,\ p=3,q=1$ とした場合(話を簡単にするため整数をビットをみなしている)のM系列の実装である。加えて、得られたM系列の先頭から $p$ ビットを 1 ビットずつずらしてリスト内にリスト化していく関数と、リスト内のリストがその自身を内包するリスト内で唯一のリストであるかを検査する関数を実装した。

import Test.HUnit
import System.IO
import Data.Bits

p :: Int
p = 3

mseq :: Int -> Int
mseq 0 = 1
mseq 1 = 0
mseq 2 = mseq 1
mseq n
  | p > q = (mseq (n - p) `xor` mseq (n - q)) :: Int
  where
    q = 1

invmseq :: [Int] -> [Int]
invmseq [] = []
invmseq (x:xs)
  | x >= 0 = mseq x : invmseq xs

isUniqueImpl :: [Int] -> [[Int]] -> Bool
isUniqueImpl x (y:z)
  | x == y = False
  | otherwise = isUniqueImpl x z
isUniqueImpl x y = [x] /= y

isUnique' :: [[Int]] -> Bool
isUnique' [] = True
isUnique' [_] = True
isUnique' (x:xs)
  | not (isUniqueImpl x xs) = False
  | otherwise = isUnique' xs

shiftpPattern :: [Int] -> [[Int]]
shiftpPattern [] = []
shiftpPattern xs
  | length xs >= p = take p xs : shiftpPattern (drop 1 xs)
  | otherwise = []


main :: IO (Counts, Int)
main = runTestText (putTextToHandle stderr False) tests
  where
    x = invmseq [0 .. 8]
    shiftp = shiftpPattern x
    isunique = isUnique' $shiftp
    tests = TestList
        [ 
            "mseq" ~: x ~?= [1, 0, 0, 1, 1, 1, 0, 1, 0],
            "shift pattern" ~: shiftp ~?= [[1, 0, 0], [0, 0, 1], [0, 1, 1], [1, 1, 1], [1, 1, 0], [1, 0, 1], [0, 1, 0]],
            "unique check" ~: isunique ~?= True
        ]

shiftpが示しているように、先頭から 1 ビットずつずらした $p$ ビットの数列は、それぞれが他の数列と一切異なる数列、つまり $p$ ビットで表現できるビットパターンを網羅している事が分かる。しかしそれは全てではなく、全てのビットパターンが 0 となる数列は発生されない*5。よって、M系列の周期 $T$ は $T=2^{p}-1$ で求まる事が言える。
この時選ばれる $p$ と $q$ の値は次の式が既約多項式となるようにして選ばれる。

$\large x^{p} + x^{q} + 1 \ \ \ \ (p \gt q)$

この時 $x$ の係数はガロア体*6 $GF(2)$ であり 0 か 1 を取る。
次にこれまでで用いた $p=3, q=1$ が、上記の条件を満たす事を背理的に証明する。次の方程式 $x^{3} + x + 1 = 0 \tag{1}$ が有理数解を持つと過程する。有理数解であるならば、互いに素な $p, q$ を用いて $\displaystyle x = \dfrac{p}{q}$ とする事ができるため、 $\displaystyle \dfrac{p^{3}}{q^{3}}+\dfrac{p}{q}+1=0$ 、 $p^{3}+pq^{2}+q^{3}=0$ 、 $p(p^{2}+q^{2})=-q^{3}$ であり、これは $q$ が $p$ の倍数である事を意味する。 $p, q$ は互いに素であるため $p = 1$ であり、 $\displaystyle \dfrac{1}{q^{3}}+\dfrac{1}{q+1}=0$ 、 $1+q^{2}+q^{3}=0$ 、 $q^{2}(q+1)=-1$ となるが、左辺は $q(q+1)$ を因数に持つため偶数、右辺は奇数となり矛盾する。よって、方程式 $(1)$ は有理数解を持たず、既約である事が言える。
[※ 追記 もっと単純化できた...
$f(x) = x^{3} + x + 1$ とする。 $f(x)$ は三次式であるため $f(x)$ が位数 2 のガロア体 $GF(2)$ 上可約であるとすると $f(x)＝(x－a)(x^{2}＋bx＋c)$ となる $a, b, c \in GF(2)$ が存在するはずである。すなわち $f(a) = 0$ である。しかし $GF(2) = {0, 1}$ で $f(0) \neq 0, f(1) \neq 0$ であるためそれに矛盾する。よって $f(x)$ は既約である。
-- 追記ここまで]

蛇足

有名なメルセンヌ・ツイスタ(MT)はこのM系列を発展させたものであり、Twisted General Feedback Shift Register(TGFSR) がはじめに発案された。これは、次の漸化式 $x_{n+p} = x_{n+q} + x_{n}A\ \ \ \ (p \gt q \gt 0)$ によって発生する。 $A$ は Twister と呼ばれる $GF(2)$ の元を要素とする $\omega \times \omega\ \ \ \ (x_{n} \in GF(2)^{\omega})$ の行列である。TGFSR はM系列と同様 $2^{p\omega} - 1$ の最大周期を持つ、すなわちM系列よりも周期が長く、かつ高速な演算となるような $A$ を選択できることから、M系列よりも優れる。MT はこの TGFSR を改良したものであり、次の漸化式 $x_{n+p}=x_{n+q}+x_{n+1}B+x_{n}C\ \ \ \ (p \gt q \gt 0) \tag{2}$ によって発生する。式 $(2)$ による最大周期が $2^{19937} - 1$ である時、特に MT19937 と呼ばれる*7。