数学アルゴリズム研究ノート『非線形合同法による擬似乱数』

数学アルゴリズム研究ノート（１）

１．非線形合同法による擬似乱数

2003年７月　田沼英樹

１－１．線形合同法について

　擬似乱数やハッシュ関数のアルゴリズムの基本として，まず何らかの整数から整数への変換（数学的には全単射）を与える関数が必要になります．

　ここでいう整数とは，数学的な意味での無限集合 $Z$ のことではなく，例えばC言語の32ビット整数に相当する $Z/2^{32}Z$ の様な $2$ の冪を法とする剰余環あるいはその元のことを指します．（以後，暗黙に仮定）

　このような関数の自明な例として，以下のものが挙げられます．

　　C言語表記数式表記

　　 x += 11 $f(x) = x+11$

　　 x *= 13 $f(x) = 13x$

　　 x ^= 15 $f(x) = x$ xor $15$ 　（正式な数式表記ではない）

　　 x = ~x $f(x) = -x-1$

　整数から整数への全単射を与える関数のうち，さらにいくつかの望ましい性質をもったものが，擬似乱数として有用なものとなり得ます．

　これらの性質をまとめると以下のようになりますが，一般に下に行くほど条件として厳しく，また検証や証明が難しくなります．

　　(0) 計算が速い．

　　(1) 全単射である．

　　(2) 周期が長い．

　　(3) 値の偏りが（実用上）十分少ない．

　このような擬似乱数として特に一次式 $f(x)=ax+c$ による変換を利用したものは，（混合型）線形合同法と呼ばれて昔から広く使われおり，また様々な研究がなされています．

　例えば多くのC言語処理系の32ビット整数上では，加減算や乗算の結果自動的に $2^{32}$ を法とする剰余が得られるために，この変換は非常に高速に行われます．ここで使われる乗算は，加減算や論理演算に比べると負荷の大きい演算ですが，最新のCPUではパイプライン処理による高速化が可能なのでそれほど低速にはなりません．

　一般に一次式 $f(x)=ax+c$ が全単射であるための必要十分条件は， $a$ が奇数であることが容易にわかります．また，さらに詳細な解析により，周期軌道が整数全体にわたるための条件が $a\equiv 1{\rm mod}4, c\equiv 1{\rm mod}2$ であることもわかります．

　しかし，これが乱数として望ましいものであるかどうかを判定するのには非常に困難が伴います．Knuth の本 [1] には多種多様な乱数検定法の紹介があり，特に線形合同法については多次元分布を調べるための強力な理論的手法として，スペクトル検定というものがあります．例えば Numerical Recipes [2] によると，32ビット整数の線形合同法による擬似乱数としては

$f(x) = 1664525x+1013904223$

で表されるものが最上のものであるとされていますが，たしかに Knuth [1] を見ると，この係数については６次元までのスペクトル検定で比較的良好な結果が得られていることがわかります．

　ただ，この種の乱数は基本的に周期がワード長に制限されてしまうこと，そして，下位ビットの周期が短い，多次元の分布があまり良くないといった弱点があることから，大規模なシミュレーション等に利用するには不十分である可能性があります．

　そこで，これらの弱点を改善するための様々な拡張について，これから考察していきたいと思います．

１－２．整数をかき混ぜる二次式

　それではまず，線形合同法の単純な拡張として，二次式 $f(x)=ax^{2}+bx+c$ による変換を利用することを考えます．

　二次式が全単射であるかどうかの判定は比較的容易にできます．

【定理１－２】

　整数係数の二次式 $\mathbf{f}(\mathbf{x})=\mathbf{a}\mathbf{x}^{\mathbf{2}}+\mathbf{b}\mathbf{x}+\mathbf{c}$ が $\mathbf{Z}/\mathbf{2}^{\mathbf{n}}\mathbf{Z}, \mathbf{n}\geq\mathbf{2}$ 上全単射であるための必要十分条件は $\mathbf{a}\equiv\mathbf{0}{\rm mod}\mathbf{2}, \mathbf{b}\equiv\mathbf{1}{\rm mod}\mathbf{2}$ ．

［証明］

　単射性を示す．（有限集合上の自己写像では単射と全単射は同値）

　以下等号は $Z/2^{n}Z$ 上の同値関係とする．

　必要性：

　　もし $b\equiv 0{\rm mod}2$ ならば $f(2^{n-1})=2^{2n-2}a+2^{n-1}b+c=c=f(0)$ より $f$ は単射でない．したがって $b\equiv 1{\rm mod}2$ ．

　　もし $a\equiv 1{\rm mod}2$ ならば $f(-a^{-1}b)=a^{-1}b^{2}-a^{-1}b^{2}+c=c=f(0)$ より $f$ は単射でない．したがって $a\equiv 0{\rm mod}2$ ．

　十分性：

　　 $f(x)=f(y)$ $\Leftrightarrow$ $ax^{2}+bx+c=ay^{2}+by+c$ $\Leftrightarrow$ $a(x^{2}-y^{2})+b(x-y)=0$ $\Leftrightarrow$ $\{a(x+y)+b\}(x-y)=0$ ．

　　 $a\equiv 0{\rm mod}2, b\equiv 1{\rm mod}2$ のとき $\{a(x+y)+b\}$ は常に奇数となるので $f(x)=f(y)$ $\Leftrightarrow$ $(x-y)=0$ ．

　　よって， $f$ は単射．■

　この定理により，例えば次の二次式は整数の全単射であることがわかります．

C言語表記： x *= ~(x<<1)

　数式表記： $f(x) = -2x^{2}-x$

　ただし，この変換は軌道周期があまり長くなく，例えば32ビット整数上で0x55555555を初期値とすると，周期 536870912（=0x20000000）で元の初期値に戻ります．また，初期値を 0 とした場合，明らかにその後の軌道も 0 が続きます．

　一方，次の場合は初期値によらず，周期 4294967296（=0x100000000）で軌道が32ビット整数を全て埋め尽くします．

C言語表記： x = ((x<<1)+1) * (x+1)

　数式表記： $f(x) = 2x^{2}+3x+1$

　このように，全単射である二次式が整数を埋め尽くす周期軌道を持つためには，さらに何らかの条件が必要となります．

１－３．非線形合同法

　一次式による変換の場合は漸化式の一般項を具体的に求めることが可能なことから，周期に関する条件は比較的容易に得ることができます．（詳細は Knuth [1] を参照）

　しかし，二次式の場合は一般項の式を直接求めることはおそらく不可能なので，何か別の方策が必要となります．そこで，求める条件を周期軌道が整数を埋め尽くすかどうかのみに絞ると，ちょっとした工夫で次の条件が得られます．

【定理１－３】

　整数係数の二次式 $\mathbf{f}(\mathbf{x})=\mathbf{a}\mathbf{x}^{\mathbf{2}}+\mathbf{b}\mathbf{x}+\mathbf{c}$ が任意の $\mathbf{n}\in\mathbf{N}$ について $\mathbf{Z}/\mathbf{2}^{\mathbf{n}}\mathbf{Z}$ 上周期 $\mathbf{2}^{\mathbf{n}}$ をもつための必要十分条件は $\mathbf{b}\equiv\mathbf{c}\equiv\mathbf{1}{\rm mod}\mathbf{2}, \mathbf{a}+\mathbf{b}\equiv\mathbf{1}{\rm mod}\mathbf{4}$ ．

［証明の概略］

　下位 $n$ ビットの周期が $2^{n}$ であると仮定した上で，延長して $n+1$ ビットにしたときの周期が $2^{n+1}$ となる条件を求める．

　詳しい証明はこちら．

　この定理の条件は， $a=0$ としたときの線形合同法の条件を含んでいます．また，実用性はともかく原理的には，この証明の手法で三次式以上にも拡張が可能です．

　実はこれをより一般化した形の定理が Knuth [1] にあり，この二次式による擬似乱数自体に新規性はありません．では，なぜこの二次合同法が線形合同法ほどよく使われていないのかというと，おそらく線形合同法の場合のスペクトル検定のような強力な理論的裏付けがないために，乱数として良いという確たる証拠がないのと，特に周期が長くなるわけでもないのに演算の負荷が増えてしまうからではないかと思われます．

　では，果たして二次合同法は擬似乱数として検討するに値しないものなのでしょうか？

　例えば，前に出た二次式 $f(x) = 2x^{2}+3x+1$ の例では乗算１回，シフト１回，加算２回と演算負荷は比較的軽く，また最新のスーパースカラー式のCPUでは ((x<<1)+1) と (x+1) の演算は並列化が期待できます．より一般の場合でも，係数によっては乗算を１回のみに抑えることが可能になります．

　そして，下位ビットの様子を視覚化すると，さらに違いがはっきり現れます．下図は，下位10ビットについて値の組 $\left(x,f(x)\right)$ をプロットしたものです．

　　　左：線形合同法　 $f(x) = 1664525x+1013904223$

　　　中：二次合同法　 $f(x) = 2x^{2}+3x+1$ 　　　右：二次合同法　 $f(x) = 4x^{2}+5x+1$

　図の左端は代表的な線形合同法のプロットで，明らかに格子模様が出現しています．乱数としてはかなり不自然な分布ですが，実際にはビット数を増やせば点の総数も増えるので，周期の大部分を埋め尽くさない限り格子模様は現れません．とはいえ，10ビット程度では十分なランダム性が得られないことをこの図は示しています．

　図の中央は最も簡単な二次合同法で，一見線形合同法よりもランダムに見えますが，斜めに不自然な線状の偏りがあります．これは， $f(16x+7)+16x+7$ $=$ $512x^{2}+512x+127$ $\equiv$ $127{\rm mod}1024$ という関係式によるものです．この偏りを除くと，線形合同法と比べて比較的少ないビット数でランダム性が得られると見ることもできます．

　図の右端は係数を少し変えて，C言語表記で x = ((x<<2)+1) * (x+1) としたものです．若干ひし形の周期が見えますが，偏りが分散している分こちらの方が良いかもしれません．

　これだけで二次合同法が乱数として優れているという証拠になり得るわけではありませんが，ここで挙げた二次式はあくまでも最も簡単な例であり，係数を改良することによってさらに良質の乱数が得られる可能性を示唆しているようにも見えます．

　次回はこの基本的なアイデアをもとに，さらに拡張した形の擬似乱数について考察します．

参考文献

[1] Knuth, D.E., Seminumerical Algorithms, 2nd ed., The Art of Computer Programming Volume II (Addison−Wesley, 1981).

　　渋谷政昭訳　『準数値算法／乱数』（サイエンス社，1981）．

[2] Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., Numerical Recipes in C: The Art of Scientific Computing (Cambridge University Press, 1988).

　　丹慶勝市・奥村晴彦・佐藤俊郎・小林誠訳　『Numerical Recipes in C ［日本語版］Ｃ言語による数値計算のレシピ』（技術評論社，1993）．

次に進む

戻る