\[ %汎用 \newcommand{\ctext}[1]{\raise0.2ex\hbox{\textcircled{\scriptsize{#1}}}} %数学 %汎用 \newcommand{\as}{{\quad\textrm{as}\quad}} \newcommand{\st}{{\textrm{ s.t. }}} \newcommand{\set}[2]{\left\{\left.#1\;\right|\;#2\right\}} \newcommand{\naturalNumbers}{\mathbb{N}} \newcommand{\integers}{\mathbb{Z}} \newcommand{\rationalNumbers}{\mathbb{Q}} \newcommand{\realNumbers}{\mathbb{R}} \newcommand{\complexNumbers}{\mathbb{C}} \newcommand{\field}{\mathbb{F}} \newcommand{\func}[2]{{#1}\left({#2}\right)} \newcommand{\argmax}{\mathop{\textrm{arg~max}}} \newcommand{\argmin}{\mathop{\textrm{arg~min}}} %集合論 \newcommand{\range}[2]{\{#1,\dotsc,#2\}} \renewcommand{\complement}{\mathrm{c}} \newcommand{\ind}[2]{\mathbbm{1}_{#1}\left(#2\right)} \newcommand{\indII}[1]{\mathbbm{1}\left\{#1\right\}} %数論 \newcommand{\abs}[1]{\left|#1\right|} \newcommand{\combi}[2]{{_{#1}\mathrm{C}_{#2}}} \newcommand{\perm}[2]{{_{#1}\mathrm{P}_{#2}}} \newcommand{\GaloisField}[1]{\mathrm{GF}\left(#1\right)} %解析学 \newcommand{\sgn}[1]{\operatorname{sgn}\left(#1\right)} \newcommand{\cl}[1]{\operatorname{cl}#1} \newcommand{\Img}[1]{\operatorname{Img}\left(#1\right)} \newcommand{\dom}[1]{\operatorname{dom}\left(#1\right)} \newcommand{\norm}[1]{\left\|#1\right\|} \newcommand{\floor}[1]{\left\lfloor#1\right\rfloor} \newcommand{\ceil}[1]{\left\lceil#1\right\rceil} \newcommand{\expo}[1]{\exp\left(#1\right)} \newcommand{\sinc}{\mathop{\textrm{sinc}}} \newcommand{\GammaFunc}[1]{\Gamma\left(#1\right)} %逆三角関数 \newcommand{\asin}[1]{\operatorname{Sin}^{-1}{#1}} \newcommand{\acos}[1]{\operatorname{Cos}^{-1}{#1}} \newcommand{\atan}[1]{\operatorname{{Tan}^{-1}}{#1}} \newcommand{\atanEx}[2]{\atan{\left(#1,#2\right)}} %微分 \newcommand{\deriv}[3]{\frac{\operatorname{d}^{#3}#1}{\operatorname{d}{#2}^{#3}}} \newcommand{\derivLong}[3]{\frac{\operatorname{d}^{#3}}{\operatorname{d}{#2}^{#3}}#1} \newcommand{\partDeriv}[3]{\frac{\operatorname{\partial}^{#3}#1}{\operatorname{\partial}{#2}^{#3}}} \newcommand{\partDerivLong}[3]{\frac{\operatorname{\partial}^{#3}}{\operatorname{\partial}{#2}^{#3}}#1} \newcommand{\partDerivIIHetero}[3]{\frac{\operatorname{\partial}^2#1}{\partial#2\operatorname{\partial}#3}} \newcommand{\partDerivIIHeteroLong}[3]{{\frac{\operatorname{\partial}^2}{\partial#2\operatorname{\partial}#3}#1}} %積分 \newcommand{\integrate}[5]{\int_{#1}^{#2}{#3}{\mathrm{d}^{#4}}#5} \newcommand{\LebInteg}[4]{\int_{#1} {#2} {#3}\left(\mathrm{d}#4\right)} %複素解析 \newcommand{\conj}[1]{\overline{#1}} \renewcommand{\Re}[1]{{\operatorname{Re}{\left(#1\right)}}} \renewcommand{\Im}[1]{{\operatorname{Im}{\left(#1\right)}}} \newcommand{\Arg}[1]{\operatorname{Arg}{\left({#1}\right)}} \newcommand{\Log}[1]{\operatorname{Log}{#1}} %ラプラス変換 \newcommand{\LPLC}[1]{\operatorname{\mathcal{L}}\left[#1\right]} \newcommand{\ILPLC}[1]{\operatorname{\mathcal{L}}^{-1}\left[#1\right]} %線形代数 \newcommand{\bm}[1]{{\boldsymbol{#1}}} \newcommand{\Span}[1]{\operatorname{span}\left[#1\right]} \newcommand{\Ker}[1]{\operatorname{Ker}\left(#1\right)} \newcommand{\rank}[1]{\operatorname{rank}\left(#1\right)} \newcommand{\inprod}[2]{\left\langle#1,#2\right\rangle} \newcommand{\HadamardProd}{\odot} \newcommand{\HadamardDiv}{\oslash} \newcommand{\matEntry}[3]{#1\left[#2\right]\left[#3\right]} \newcommand{\matPart}[5]{\matEntry{#1}{#2:#3}{#4:#5}} \newcommand{\diag}[1]{\operatorname{diag}\left(#1\right)} \newcommand{\tr}[1]{\operatorname{tr}{\left(#1\right)}} %ベクトル %単位ベクトル \newcommand{\vix}{\bm{i}_x} \newcommand{\viy}{\bm{i}_y} \newcommand{\viz}{\bm{i}_z} %確率論 \newcommand{\PDF}[2]{\operatorname{PDF}\left[#1,\;#2\right]} \newcommand{\Ber}[1]{\operatorname{Ber}\left(#1\right)} \newcommand{\Beta}[2]{\operatorname{Beta}\left(#1,#2\right)} \newcommand{\ExpDist}[1]{\operatorname{ExpDist}\left(#1\right)} \newcommand{\ErlangDist}[2]{\operatorname{ErlangDist}\left(#1,#2\right)} \newcommand{\PoissonDist}[1]{\operatorname{PoissonDist}\left(#1\right)} \newcommand{\GammaDist}[2]{\operatorname{Gamma}\left(#1,#2\right)} \newcommand{\cind}[2]{\ind{#1\left| #2\right.}} %条件付き指示関数 \renewcommand{\Pr}[1]{\operatorname{Pr}\left[#1\right]} \newcommand{\cPr}[2]{\Pr{#1\left| #2\right.}} \newcommand{\E}[2]{\operatorname{E}_{#1}\left[#2\right]} \newcommand{\cE}[3]{\E{#1}{\left.#2\right|#3}} \newcommand{\Var}[2]{\operatorname{Var}_{#1}\left[#2\right]} \newcommand{\Cov}[2]{\operatorname{Cov}\left[#1,#2\right]} \newcommand{\CovMat}[1]{\operatorname{Cov}\left[#1\right]} %グラフ理論 \newcommand{\neighborhood}{\mathcal{N}} %プログラミング \newcommand{\plpl}{\mathrel{++}} \newcommand{\pleq}{\mathrel{+}=} \newcommand{\asteq}{\mathrel{*}=} \]

LFSRによる擬似ランダムビット列生成

経緯

 筆者はソフト屋だが、仕事でディジタル無線のベースバンド信号処理に関わっている。 この世界では通信路,送受信機の性能の測定に疑似ランダムビット列を使うのは普通のことらしい。 ビット列の生成に Linear Feedback Shift Register (線形帰還シフトレジスタ, LFSR)が使われることを知った。 機構が単純なのでソフトウェアに頼らず論理回路だけで極めて高速にビット列を生成できる。

 いつものことながら背後にある理論が気になってしまうわけで、生成多項式が原始多項式であれば周期が最長となる理由について随分と考えてしまった。自分なりに考えたことを書き残しておく。

 こちらのページで紹介されている周期7のLFSRを例にとる。

周期7のLFSR
図1: 周期7のLFSR
上の図に於いて$y(n) \in \GaloisField{2}$(位数2の有限体)であり、Dは遅延素子、$\bigoplus$はXORを表す。

出力の差分方程式

 $\bigoplus$の可換性に注意し、このシステムの出力は次の差分方程式で表される。 \[ y(n) = y(n-2) + y(n-3) \]

z変換と生成多項式

 初期条件を$y(-3)=1,\;y(m)=0\;(m<0,m \neq -3)$とする(こうすると後で式が綺麗になる)。 これは遅延素子(レジスタ)の初期状態を左から0,0,1に定め、初回のレジスタシフト処理の出力を$y(0)$とすることで達成できる。 この条件の下で上式をz変換すると \begin{align*} Y(z) &= z^{-2}Y(z) + z^{-3}Y(z) + 1 \\ (1 + z^{-2} + z^{-3})Y(z) &= 1 \\ Y(z) &= (1 + z^{-2} + z^{-3})^{-1} \end{align*} 上式の最後に現れる$(1 + z^{-2} + z^{-3})^{-1}$は$z^{-1}$に関する$\GaloisField{2}$上の形式的冪級数$\GaloisField{2}[[z^{-1}]]$に於ける$1 + z^{-2} + z^{-3}$の乗法逆元である。 この$1 + z^{-2} + z^{-3}$に於ける$z^{-1}$を$x$(記号は重要ではない。負冪を消せれば何でも良い)に置き換えたものをLFSRの生成多項式と呼ぶ。 図1の例では$x^3+x^2+1$である。 生成多項式を機械的に求めるには、$y(n)$からの入力を受けている遅延素子(図1の例では右から3番目と2番目)を出力から遠い順に$x$の高次の項としながら列挙し、最後に定数項1を足せばよい。

出力の計算と周期性の確認

 $Y(z)$を$z^{-1}$に関する形式的冪級数で表示すれば$y(n)$が求まる。 具体的には \[ Y(z) = \sum_{i=0}^\infty (z^{-2} + z^{-3})^i \] を計算すれば良い(これは理論考察の為の式であり、実際に計算するならこれよりも筆算の方が良い。このとき除数の「最高次」の項は「1」(次数0)を意味することに注意)。 最初のいくつかの項をMathematicaで計算すると次のようになる。 $\newcommand{\redtext}[1]{\textcolor{red}{#1}}$ \begin{align*} Y(z) &=\phantom{+} \redtext{1} + 0z^{-1} + \redtext{1}z^{-2} + \redtext{1}z^{-3} + \redtext{1}z^{-4} + 0z^{-5} + 0z^{-6} \\ &\phantom{=}+ \redtext{1}z^{-7} + 0z^{-8} + \redtext{1}z^{-9} + \redtext{1}z^{-10} + \redtext{1}z^{-11} + 0z^{-12} + 0z^{-13} \\ &\phantom{=}+ \redtext{1}z^{-14} + 0z^{-15} + \redtext{1}z^{-16} + \redtext{1}z^{-17} + \redtext{1}z^{-18} + 0z^{-19} + 0z^{-20} \\ &\phantom{=}+ \redtext{1}z^{-21} + 0z^{-22} + \cdots \end{align*} 係数列が1011100という周期7$(=2^3-1)$の繰り返しになっていることがわかる。

最長LFSR

 レジスタ長(遅延素子の個数)を$n$とし、レジスタの状態を$n$桁の2進数で表現することにする。 但し出力から最も遠い遅延素子をMSBとする。 レジスタの状態は高々$2^n$通りしかない。 そして一度$00 \dots 0$の状態になると永久にそこに留まることは簡単にわかる。 しかし先述のz変換から明らかなように、初期状態を$00 \dots 0$にしない限りこの状態にはならない。 また、ある時刻$n$のレジスタの状態は直前の時刻$n-1$のレジスタの状態のみに依存して決まる。 以上より、LFSRの出力すなわち$Y(z)$の係数が周期的であること、およびその周期が高々$2^n-1$であることがわかる。

 工学の立場からは、レジスタ長$n$のLFSRを設計するならば、その周期が最長($2^n-1$)となるようにしたい。 これは特に「最長LFSR」と呼ばれる。 実は、生成多項式が原始多項式であれば周期が最長($2^n-1$)になることが知られている。 詳細は次で述べる。

生成多項式に原始多項式を用いると周期が最長になる理由

 次の主張に於いて$p=2$とし、$x$を$z^{-1}$で置き換えることで先述のz変換に於ける$Y(z)$の係数列の周期が$2^n-1$であることがわかる。

 主張: 「$p$を素数とする。 $n$次の原始多項式$f(x) \in \GaloisField{p}[[x]]$の逆数$f(x)^{-1}$の係数列は周期的であり、その周期は$p^n-1$である。」

$\textit{Proof}$
 まず係数列が周期的であることを示す。 $f(x)$が原始多項式だから、適当な$p^n-1-n$次多項式$g(x) \in \GaloisField{p}[[x]]$が存在して次式が成り立つ。 \begin{align*} x^{p^n-1} - 1 &= g(x)f(x) \\ \therefore f(x)^{-1} &= -g(x)\left(1 - x^{p^n-1}\right)^{-1} \\ &= -g(x) \sum_{i=0}^\infty x^{i(p^n-1)} \end{align*} $\deg{g(x)} < p^n - 1$だから右辺の係数列は周期的であり、その周期は高々$p^n-1$である。 「高々」と言ったのは、$g(x)$の係数列のパターンに依ってはより短い周期になり得るからである。 この周期が真に$p^n-1$であることを背理法で示す。

 仮に周期が$p^n-2$以下であるとする、つまりある$l\in\{1,2,\cdots,p^n-2\}$と適当な$l-1$次以下の多項式$h(x) \in \GaloisField{p}[[x]]$が存在して \[ f(x)^{-1} = -h(x) \sum_{i=0}^\infty x^{il} \] が成り立つと仮定すると次式を得る。 \begin{align*} f(x)^{-1} &= -h(x) \sum_{i=0}^\infty x^{il} = h(x)(x^l - 1)^{-1} \\ \therefore x^l - 1 &= h(x)f(x) \equiv 0 \mod f(x) \end{align*} これは$f(x)$が原始多項式であることに矛盾する。

$\square$

もう一つの生成多項式

 次数$n$の原始多項式$f_1(x)$が一つ与えられたとき(すなわちそれに対応する最長LFSRが一つ与えられたとき)、同じ次数のもう一つの原始多項式$f_2(x)$が(すなわち遅延素子の個数が等しいもう一つの最長LFSRが)直ちに得られる。 $f_1(x)$の係数の順番を逆にしたもの(或いは、同じことだが、$f_1(x)$の$x^i\;(i=0,1,\cdots,n)$の項を$x^{n-i}$に置き換える)を$f_2(x)$とすると、これが原始多項式になる。 証明は簡単なのでここで述べておく。 次に述べる補題2で$p=2$とすれば良い。

補題1

 主張: 「体$F$上の$n$次多項式$f_1(x) := \sum_{i=0}^n a_ix^i$が$F$上で既約であれば、係数の順番を逆にした多項式$f_2(x) := \sum_{i=0}^n a_{n-i}x^i$もまた既約である。」

$\textit{Proof}$
 背理法で示す。 $f_2(x)$が既約でないと仮定すると、適当な多項式$g_{2,1}(x),g_{2,2}(x) \in F[x],\;d_{2,1} + d_{2,2} := \deg(g_{2,1}(x)) + \deg(g_{2,1}(x)) = n$が存在して次式が成り立つ。 \[ f_2(x) = g_{2,1}(x)g_{2,2}(x) \] 上式をLaurent多項式に拡張して両辺の$x$を$x^{-1}$で置き換えると \[ f_2(x^{-1}) = g_{2,1}(x^{-1})g_{2,2}(x^{-1}) \] 両辺に$x^n$を掛けると \begin{align*} f_1(x) &= x^n f_2(x^{-1}) \\ &= x^{d_{2,1}}g_{2,1}(x^{-1})x^{d_{2,2}}g_{2,2}(x^{-1}) \\ &=: g_{1,1}(x)g_{1,2}(x) \end{align*} ここに$g_{1,1}(x), g_{1,2}(x)$は$g_{2,1}(x), g_{2,2}(x)$の係数の順番を逆にしたものである。 これは$f_1(x)$が既約であることに矛盾する。

$\square$

補題2

 主張: 「$p$を素数とする。 $n$次多項式$f_1(x) := \sum_{i=0}^n a_ix^i \in \GaloisField{p}[x]$が原始多項式であれば、係数の順番を逆にした多項式$f_2(x) := \sum_{i=0}^n a_{n-i}x^i$もまた原始多項式である。」

$\textit{Proof}$
 補題1より$f_2(x)$は既約である。 あとは$f_2(x) \mid x^l - 1$となる最小の自然数$l$が$p^m-1$であることを示せば良い。

 まず$f_2(x) \mid x^{p^n-1}-1$を示す。 $f_1(x)$が原始多項式だから適当な$p^n-1-n$次多項式$g_1(x) \in \GaloisField{p}[[x]]$が存在して次式が成り立つ。 \[ x^{p^n-1} - 1 = g_1(x)f_1(x) \] 上式をLaurent多項式に拡張して両辺の$x$を$x^{-1}$で置き換えると \[ \frac{1}{x^{p^n-1}} - 1 = g_1(x^{-1})f_1(x^{-1}) \] 両辺に$x^{p^n-1}$を掛けると \begin{align*} 1 - x^{p^n-1} &= x^{p^n-1-n}g_1(x^{-1}) x^n f_1(x^{-1}) \\ &=: g_2(x)f_2(x) \\ \therefore x^{p^n-1} -1 &\equiv 0 \mod f_2(x) \end{align*} ここに$g_2(x)$は$g_1(x)$の係数の順番を逆にしたものである。

 次に$l=1,2,\cdots,p^n-2$について$f_2(x) \nmid x^l-1$を背理法で示す。 仮にある$l$に対して$l-n$次多項式$h_1(x) \in \GaloisField{p}[x]$が存在して \[ x^l - 1 = h_1(x)f_2(x) \] が成り立つとする。 上式をLaurent多項式に拡張して両辺の$x$を$x^{-1}$で置き換えると \[ \frac{1}{x^l} - 1 = h_1(x^{-1})f_2(x^{-1}) \] 両辺に$x^l$を掛けると \begin{align*} 1 - x^l &= x^{l-n}h_1(x^{-1}) x^n f_2(x^{-1}) \\ &=: h_2(x)f_1(x) \\ \therefore x^l - 1 &\equiv 0 \mod f_1(x) \end{align*} ここに$h_2(x)$は$h_1(x)$の係数の順番を逆にしたものである。 これは$f_1(x)$が原始多項式であることに矛盾する。

$\square$

出力の順番が逆転すること

 この方法で得られる生成多項式に対応するLFSRの出力に現れる周期的なパターンは元のLFSRのちょうど逆順になる。 図1の例の生成多項式$x^3+x^1+1$の係数の順番を反転した$x^3+x+1$に対応するLFSRの出力を$y_2(n)$とする。 そのz変換を計算すると次式になる。 \begin{align*} Y_2(z) &=\phantom{+} (1 + z^{-1} + z^{-3})^{-1} \\ &=\phantom{+} \redtext{1} + \redtext{1}z^{-1} + \redtext{1}z^{-2} + 0z^{-3} + \redtext{1}z^{-4} + 0z^{-5} + 0z^{-6} \\ &\phantom{=}+ \redtext{1}z^{-7} + \redtext{1}z^{-8} + \redtext{1}z^{-9} + 0z^{-10} + \redtext{1}z^{-11} + 0z^{-12} + 0z^{-13} \\ &\phantom{=}+ \redtext{1}z^{-14} + \redtext{1}z^{-15} + \redtext{1}z^{-16} + 0z^{-17} + \redtext{1}z^{-18} + 0z^{-19} + 0z^{-20} \\ &\phantom{=}+ \redtext{1}z^{-21} + \redtext{1}z^{-22} + \cdots \end{align*} 任意の最長LFSRについてこれが成り立つことは、次に述べる補題で$p=2$とし、$x$を$z^{-1}$で置き換えればわかる。

補題3

 主張: 「$p$を素数とする。 $n$次の原始多項式$f_1(x) \in \GaloisField{p}[[x]]$の係数の順番を逆転させ$f_2(x)$とする。 $f_1(x)^{-1} =: a_0 + a_1x + \cdots$とする。 このとき次式が成り立つ。 \[ f_2(x)^{-1} = -\sum_{j=0}^{p^n-1-n} a_{p^n-1-n-j}x^j \sum_{i=0}^\infty x^{i(p^n-1)} \] 視覚的に言えば、$f_2(x)^{-1}$の係数列は周期的であり、$f_1(x)^{-1}$の係数列に現れるパターンを逆転させ、符号を逆にしたものになる。」

$\textit{Proof}$
 既に述べたように$f_1(x)^{-1}$の係数列は周期$p^n-1$であるから、適当な$p^n-1$次未満の多項式$g_1(x) \in \GaloisField{p}[[x]]$が存在して次式が成り立つ。 \[ f_1(x)^{-1} = g_1(x) \sum_{i=0}^\infty x^{i(p^n-1)} \] この式をLaurent多項式に拡張して変形してゆく。 \begin{align*} f_1(x)^{-1} &= g_1(x) \sum_{i=0}^\infty x^{i(p^n-1)} \\ &= -g_1(x) \left(x^{p^n-1} - 1\right)^{-1} \\ \therefore x^{p^n-1} - 1 &= -g_1(x)f_1(x) \\ &= -g_1(x) x^n f_2(x^{-1}) \end{align*} ここに$f_2(x)$は$f_1(x)$の係数の順番を逆転させたものである。 初めに$\deg g_1(x) < p^n-1$としていたが、上式より$\deg g_1(x) = p^n-1-n$が確定する。 上式の$x$を$x^{-1}$で置き換えると次式を得る。 \[ x^{-(p^n-1)} - 1 = -g_1(x^{-1}) x^{-n} f_2(x) \] 両辺に$x^{p^n-1}$を掛けると \[ 1 - x^{p^n-1} = -x^{p^n-1-n}g_1(x^{-1}) f_2(x) =: -g_2(x)f_2(x) \] ここに$g_2(x)$は$g_1(x)$の係数の順番を逆転させたものである。 これより次式を得る。 \[ f_2(x)^{-1} = -g_2(x) \left(1 - x^{p^n-1}\right)^{-1} = -g_2(x) \sum_{i=0}^\infty x^{i(p^n-1)} \]

$\square$

実装

 最近私が関心を寄せているプログラミング言語Rustを使って実装してみた。

周期7 ($=2^3-1$)

 次に示すのは図1の例を実装したものである。

 先述の通り上の例の生成多項式の係数の順番を逆転させた$x^3+x+1$に対応するLFSRも周期7をもつ。 これは上のコードのconst XOR_MASK: u8 = 0b110const XOR_MASK: u8 = 0b101に書き換えることで実現される。

周期511 ($=2^9-1$)

 次に示すのは生成多項式が$x^9+x^5+1$である最長LFSRを実装したものである。

 周期7の例と同様、生成多項式の係数の順番を逆転させた$x^9+x^4+1$に対応するLFSRも周期511をもつ。 これは上のコードのconst XOR_MASK: u16 = 0b100_010_000const XOR_MASK: u16 = 0b100_001_000に書き換えることで実現される。

スポンサーサイト



コメント

非公開コメント

プロフィール

motchy

Author:motchy
・信号処理
・無線
・組込

GitHub/motchy869

検索フォーム