線形合同法数列が最大周期になる条件(ただし法が2のべきの場合)

はじめに

一部界隈の皆さんはよくご存知の線形合同法数列
   a_0 = F \bmod 2^{32} (Fは初期シード)
   a_{n+1} = ({\rm 0x41c64e6d} \times a_n + {\rm 0x6073}) \bmod 2^{32}

は最大周期となることが知られています。
すなわちa_1, a_2, ... と進めていって a_{2^{32}}で初めてa_0と等しくなります。
なぜそんなことがいえるのでしょうか。またどういう条件のときに最大周期になるのでしょうか。

問題

線形合同法数列
   a_0 = F \bmod M
   a_{n+1} = (A a_n + B) \bmod M
が最大周期となる条件はなにか。

この記事ではM4以上の2の累乗のときを扱いましょう。つまりs2以上の整数としてM = 2^sと表せるときです。
(M = 2のときも含めると後に場合分けが必要になるので省きました。ご都合主義です^^;)

具体例

まずは実際にA, B, M, Fに具体的な値を入れて数列を眺めてみましょう。


  • a_0 = 0, a_{n+1}=(5 a_n + 1) \bmod 8

    n 0 1 2 3 4 5 6 7 8 9
    a_n 0 1 6 7 4 5 2 3 0 1

    …最大周期となる


  • a_0 = 0, a_{n+1}=(5 a_n + 2) \bmod 8

    n 0 1 2 3 4 5 6 7 8 9
    a_n 0 2 4 6 0 2 4 6 0 2

    …最大周期とならない


  • a_0 = 1, a_{n+1}=(5 a_n + 2) \bmod 8

    n 0 1 2 3 4 5 6 7 8 9
    a_n 1 7 5 3 1 7 5 3 1 7

    …最大周期とならない


  • a_0 = 0, a_{n+1}=(4 a_n + 1) \bmod 8

    n 0 1 2 3 4 5 6 7 8 9
    a_n 0 1 5 5 5 5 5 5 5 5

    …最大周期とならない

F=0としても一般性は失われない

実はA, B, M3つが決まれば、初期シードFによらず最大周期になるかどうかは決まります。
これを証明すれば、後の議論においてF = 0の場合だけ考えればよいことになり楽できます。

それを証明するための準備として周期に関する性質をまとめておきます。

周期に関する性質

具体的なところから、線形合同法数列\< a_n\>a_4で初めてa_0の値に戻ってくるものとします。数列の項の具体的な値は今は重要ではないのでa, b, c, dとしましょう。

ここで重要になってくるのは線形合同法数列がもつ、a_i = a_jならばa_{i+1} = a_{j+1}という性質です。

数列の項の変化を図に表すと次のようになりますね。

ここでa_0, a_1, a_2, a_3の値にダブりはあるでしょうか? ありません。
まず、a_1, a_2, a_3のいずれかがa_0に等しいとすると「a_4で初めてa_0の値に戻る」ということに反するのでありえません。
次にa_1, a_2, a_3の中でダブりがあるとするとそこでループができてしまって二度とa_0の値には戻ってきません。
たとえばa_1 = a_3のとき次のようになります。

以上の議論は周期が4以外のときも成り立ちます。これは直感的理解のための説明です。ちゃんとした証明も難しくはありませんが、この記事では省略します。

上記の説明の通り、線形合同法数列a_na_Mで初めてa_0の値に戻ってくるとき、a_0, a_1, a_2, \cdots, a_{M-1}は互いに異なります。互いに異なるのだからこれらはM個の値からなります。
ここで数列a_nの項はすべて0以上M未満であるM個の整数の中からとるのですから、a_0, a_1, a_2, \cdots, a_{M-1}には0以上M未満の整数が一つずつ現れるといえます。

F=0としても一般性は失われないことの証明

任意のf_1, f_2について、F = f_1で最大周期となるならばF = f_2で最大周期となることを証明すればよいです。

なぜならこれが正しいとすると、

  • f_1=0, f_2=fをあてはめて、F = 0で最大周期ならばF = fで最大周期であること
  • f_1=f, f_2=0をあてはめて、F = fで最大周期ならばF = 0で最大周期であること

がいえますから、任意のfについて、F = 0で最大周期であることとF = fで最大周期であることが同値になります。

さて、任意のf_1, f_2について、F = f_1で最大周期となるならばF = f_2で最大周期となることを証明しましょう。
F=f_1のときの数列を\< a_n\>, F=f_2のときの数列を\< b_n\>とします。

\< a_n\>は最大周期となるのだから、「周期に関する性質」よりa_i = f_2 \bmod M, 0 \le i < Mとなる整数iが存在しますね。
a_i = b_0 であり、二つの数列の漸化式は同じなのだから任意のkについてa_{i+k} = b_kといえます。

えーここでもちゃんとした証明がめんどくさくなったので^^;、図を使った直感的な説明に逃げます。

これを見るとa_iからa_{i+M-1}までにダブりがない、すなわちb_0からb_{M-1}までにダブりがない。また、a_i = a_{i+M}すなわちb_0 = b_Mだと直感的理解ができますね。

よって数列\< b_n\>は最大周期になります。
したがって、任意のf_1, f_2について、F = f_1で最大周期となるならばF = f_2で最大周期となるといえます。

F=0としても一般性は失われないことが示されましたので、以降からはF=0の場合だけ考えましょう。

一般項

まずa_nは必ず0以上M未満の整数となりますから a_n = (a_n \bmod M)が常に成り立ちます。
ここで数列a_nの一般項を求めてみましょう。

実際に最初の数項を観察してみると

いずれもMを法として
   a_0 \equiv 0
   a_1 \equiv A \times a_0 + B \equiv B
   a_2 \equiv A \times a_1 + B \equiv A \times B + B = (A + 1) B
   a_3 \equiv A \times a_2 + B \equiv A \times (A + 1) B + B = (A^2 + A + 1) B

より
   a_n \equiv (1 + A + \cdots + A^{n-1}) B \pmod{M}
と予想できます。この予想が正しければ、

   a_n \bmod M = ( (1 + A + \cdots + A^{n-1}) B) \bmod M
   よって a_n = ( (1 + A + \cdots + A^{n-1}) B) \bmod Ma_n = a_n \bmod M より)

となります。
つまり一般項は a_n = ( (1 + A + \cdots + A^{n-1}) B) \bmod M と予想できますね。
この予想は数学的帰納法で簡単に証明できます。(省略)

補足

ここではさらっと次の定理を使いました。
a \equiv b \pmod{m}ならばa + c \equiv b + c \pmod{m}ac \equiv bc \pmod{m}
この定理はm=10のときで考えてみるととっても当たり前です。
足し算や掛け算した結果の一の位は、足す数や掛ける数の一の位だけで決まりますからね。

Bが偶数のとき?

まずは簡単なところから条件を絞っていきます。Bが偶数であるとするとどうなるでしょう。
最初の具体例の一つにもう一度ご登場願いましょう。

a_0 = 0, a_{n+1}=(5 a_n + 2) \bmod 8

n 0 1 2 3 4 5 6 7 8 9
a_n 0 2 4 6 0 2 4 6 0 2

Mも偶数なのでa_nは常に偶数となってしまいます。(a_nが(偶数)をMで割った余りでありa_n = (偶数) - M \times (商)と表せるから)
しかし、最大周期のときa_0, a_1, \cdots, a_{M-1}0からM-1までのすべての整数をとるはずだったので、Bが偶数だとこれに反しますね。
よって、Bが偶数のときは最大周期になりません。

Aが偶数のとき?

では、Aが偶数のときはどうでしょう。

Aが偶数かつ最大周期だと仮定します。
ここで周期に関する性質よりa_i = M/2, 0 \le i < Mを満たす整数iがあります。M/2\ne 0だからi=0はないですね。
A2の倍数だからA a_iMの倍数となります。
Mの倍数の違いはMで割った余りには影響しませんから、
   \begin{eqnarray}\\ a_{i+1} &=& (A a_i + B) \bmod M\\ &=& B \bmod M\\ \end{eqnarray}
よってa_{i+1} = a_1 です。a_1, a_2, ..., a_Mは互いに異なりますから矛盾しました。(これはa_0, a_1, ..., a_{M-1}が互いに異なることとa_M = a_0からすぐ分かります)


補足
二進数で考えてみるとAが偶数だったら現在の項の最上位ビットが次の項にまったく影響しません。すると最上位ビットが異なりほかのビットは一致する2つのシード値において、それぞれの次のシード値が一致することになります。このとき直感的に明らかに最大周期にはなりませんね。このことをちゃんと証明したのが上の内容です。

さらに補足
a_0 = 0, a_{n+1}=(4 a_n + 1) \bmod 8

n 0 1 2 3 4 5 6 7 8 9
a_n 0 1 5 5 5 5 5 5 5 5

この具体例のようにAが偶数だった場合は必ず一定値に収束してしまいます。擬似乱数としてはまったく使い物になりませんね。
余談ですが、あるゲームのあるデータの難読化用XORマスクとしてAが偶数である線形合同法数列が使われていました^^;この例では別に擬似乱数としてぶっ壊れていても全然差し支えはないですが、もしゲーム中のランダム要素でこうなっていたら大変なことになりますね

さて、ABのどちらか少なくとも一方が偶数のときは最大周期にならないことは分かったので、以降ではABともに奇数である場合を調べましょう。

a_n = 0という条件

最大周期となるというのは、a_Mで初めてa_0の値に戻ってくるということでした。
つまりa_M = 0かつ1 \le n < M であるすべてのnについてa_n \ne 0であるということです。

ということはnに対してa_n0になるかどうかが重要になります。a_n = 0であるかどうかという条件を変形していきましょう。

   a_n = 0
   \Leftrightarrow ( (1 + A + \cdots + A^{n-1}) B) \bmod M = 0
   \Leftrightarrow (1 + A + \cdots + A^{n-1}) BMの倍数
   \Leftrightarrow (1 + A + \cdots + A^{n-1}) B2^sの倍数
   \Leftrightarrow (1 + A + \cdots + A^{n-1})2^sの倍数 (Bは奇数より)

よって (1 + A + \cdots + A^{n-1})2の何乗で割れるかに注目しましょう。

正整数x2^kで割り切れるような最大のkv_2(x) で表すことにします。

D_A(n) = (1 + A + \cdots + A^{n-1}) とおきます。

すると
   a_n = 0 \Leftrightarrow v_2(D_A(n)) \ge s
となります。(ただしn>0とする。n=0のときD_A(n)=0v_2(0)は求まらないから)

v_2(D_A(n))を観察する

ここで v_2(D_A(n)) が実際どうなるのか具体的に見てみましょう。

A=3のとき

n D_A(n) v_2(D_A(n))
1 1=1 0
2 1+3=4 2
3 1+3+3^2=13 0
4 1+3+3^2+3^3=40 3
5 1+3+3^2+3^3+3^4=121 0
6 1+3+3^2+3^3+3^4+3^5=364 2
7 1+3+3^2+3^3+3^4+3^5+3^6=1093 0
8 1+3+3^2+3^3+3^4+3^5+3^6+3^7=3280 4
9 1+3+3^2+3^3+3^4+3^5+3^6+3^7+3^8=9841 0
10 1+3+3^2+3^3+3^4+3^5+3^6+3^7+3^8+3^9=29524 2

v_2(D_A(n))になにか規則性はないか睨んでみると、v_2(D_A(n))v_2(n)と関係がありそうだとわかります。

n v_2(n) v_2(D_A(n))
1 0 0
2 1 2
3 0 0
4 2 3
5 0 0
6 1 2
7 0 0
8 3 4
9 0 0
10 1 2

nが偶数のときはv_2(D_A(n)) = v_2(n) + 1nが奇数のときはv_2(D_A(n)) = 0となっていそうです。

nが奇数のとき0になるのは当たり前ですね。奇数を奇数個集めた和は奇数になりますから。
よってnが偶数の場合のみを追ってA=3以外も試してみます。さすがにAが大きくなると計算が大変になるのでコンピュータに計算させました^^;
v_2(D_A(n))の表

An 2 4 6 8 10
1 1 2 1 3 1
3 2 3 2 4 2
5 1 2 1 3 1
7 3 4 3 5 3
9 1 2 1 3 1
11 2 3 2 4 2
13 1 2 1 3 1
15 4 5 4 6 4
17 1 2 1 3 1
19 2 3 2 4 2
21 1 2 1 3 1
23 3 4 3 5 3
25 1 2 1 3 1
27 2 3 2 4 2
29 1 2 1 3 1
31 5 6 5 7 5

\< v_2(2), v_2(4), v_2(6), v_2(8), v_2(10)\> = \< 1,2,1,3,1\>ですから、何かにv_2(n)を足した結果がv_2(D_A(n))になっていそうですね。
v_2(n)を引いた結果を表にしてみると

v_2(D_A(n)) - v_2(n)の表

An 2 4 6 8 10
1 0 0 0 0 0
3 1 1 1 1 1
5 0 0 0 0 0
7 2 2 2 2 2
9 0 0 0 0 0
11 1 1 1 1 1
13 0 0 0 0 0
15 3 3 3 3 3
17 0 0 0 0 0
19 1 1 1 1 1
21 0 0 0 0 0
23 2 2 2 2 2
25 0 0 0 0 0
27 1 1 1 1 1
29 0 0 0 0 0
31 4 4 4 4 4

15のとき331のとき4…、152^4-1312^5-1だから…と考えるとこの数値はv_2(A+1)-1に等しいと気づきます。つまりv_2(D_A(n)) - v_2(n) = v_2(A+1) - 1です。

まとめますと、

nが奇数のとき v_2(D_A(n)) = 0
nが偶数のとき v_2(D_A(n)) = v_2(n) + v_2(A+1) - 1

ということです。

しかしこれはAnが小さいときの結果からの予想に過ぎないので、実際にすべてのAnで成り立つことは証明しないといけません。

偶数n>0についてv_2(D_A(n)) = v_2(n) + v_2(A+1) - 1の証明

さあ、ここが一番やっかいそうなところです。

関数v_2(x)v_2(x \times y) = v_2(x) + v_2(y)という対数関数のような性質を持ちます。
i = v_2(x), j = v_2(y)とするとm, nを奇数としてx = 2^i m, y = 2^j nと表せますから、x \times y = (2^i m) \times (2^j n) = 2^{i+j} (mn)mnは奇数よりv_2(x \times y) = i + jとなるからです。

変形

n = 2k とおいて、証明すべき命題を変形しましょう。
   v_2(D_A(2k)) = v_2(2k) + v_2(A+1) - 1
v_2(2k) = v_2(2) + v_2(k) を使って
   \Leftrightarrow v_2(D_A(2k)) = (1 + v_2(k)) + v_2(A+1) - 1
等式の右辺を整理して
   \Leftrightarrow v_2(D_A(2k)) = v_2(k) + v_2(A+1)
両辺にv_2(A-1)を足して
   \Leftrightarrow v_2(D_A(2k)) + v_2(A-1) = v_2(k) + v_2(A+1) + v_2(A-1)
v_2(x) + v_2(y) = v_2(x \times y)を使って
   \Leftrightarrow v_2(D_A(2k) \times (A-1)) = v_2(k) + v_2( (A+1) \times (A-1))
ここで
   \begin{eqnarray}\\ D_A(2k) \times (A-1) &=& (1 + A + \cdots + A^{2k-1}) \times (A - 1)\\ &=& A^{2k} - 1\\ \end{eqnarray}
だから (等比数列の和の公式の分母を払った形ですね)
   \Leftrightarrow v_2(A^{2k} - 1) = v_2(k) + v_2(A^2-1)
さらにc=A^2とおくと
   \Leftrightarrow v_2(c^{k} - 1) = v_2(k) + v_2(c-1)
となります。

おっと、A=1のときはv_2(A-1)=v_2(0)は求まりませんから、A=1のときは別に議論をちゃちゃっと済ませちゃいましょう。
A=1のときD_A(n)=nだからv_2(D_A(n)) = v_2(n)。また、v_2(A+1) - 1 = v_2(2) - 1 = 0だからv_2(D_A(n)) = v_2(n) + v_2(A+1) - 1が成り立ちます。

具体的なkの値で考えてみる

具体的なkの値としてk=10で考えてみます。

c^{10} - 1因数分解します。
   \begin{eqnarray}\\ c^{10} - 1 &=& (c^5 - 1) \times (c^5 + 1)\\ &=& (c - 1) \times (1 + c + c^2 + c^3 + c^4) \times (c^5 + 1)\\ \end{eqnarray}
より
   v_2(c^{10} - 1) = v_2(c - 1) + v_2(1 + c + c^2 + c^3 + c^4) + v_2(c^5 + 1)
です。

すると、証明すべき等式
   v_2(c^{10} - 1) = v_2(c - 1) + v_2(10)
をいうためには、二つの式を見比べると
   v_2(1 + c + c^2 + c^3 + c^4) + v_2(c^5 + 1) = v_2(10)
つまり、v_2(10) = 1より
   v_2(1 + c + c^2 + c^3 + c^4) + v_2(c^5 + 1) = 1
がいえればいいですね。

v_2(1 + c + c^2 + c^3 + c^4)はすぐわかります。Aは奇数よりc = A^2は奇数だから奇数を奇数個集めた和は奇数より、これは0になります。
するとあとはv_2(c^5 + 1) = 1がいえればいいはずです。v_2の値が1というのは2で割り切れ4で割り切れないということですから4で割った余りを考えましょう。

c=A^2であり、Aは奇数です。
奇数の2乗は必ず4で割った余りが1になります。また4で割った余りが1である数同士の積は再び4で割った余りが1になります。(これは奇数を2n+1とおいたり、4で割った余りが1である二数を4n+14m+1とおいたりして実際に計算してみれば明らかです)
よってcの累乗は常に4で割った余りが1になります。

よってc^5 \bmod 4 = 1より(c^5 + 1) \bmod 4 = 2です。よってv_2(c^5 + 1) = 1がいえます。
これでk=10のとき証明すべきv_2(c^k - 1) = v_2(c - 1) + v_2(k)がいえました。

また、同じようにしてk=5のときの
   \begin{eqnarray}\\ v_2(c^5 - 1) &=& v_2(c - 1) + v_2(1 + c + c^2 + c^3 + c^4)\\ &=& v_2(c - 1) + 1\\ &=& v_2(c - 1) + v_2(5)\\ \end{eqnarray}
もいえます。


そうなるとk=10の場合はk=5を先に証明しておいてから、その結果を使えばいいですね。

    v_2(c^{10} - 1) = v_2(c^5 - 1) + v_2(c^5 + 1)
よりk=5の結果 v_2(c^5 - 1) = v_2(c - 1) + v_2(5)を使って
   \begin{eqnarray}\\ v_2(c^{10} - 1) &=& v_2(c - 1) + v_2(5) + v_2(c^5 + 1)\\ &=& v_2(c - 1) + v_2(5) + 1\\ &=& v_2(c - 1) + v_2(10)\\ \end{eqnarray}

同様にk=20の場合は、k=10の結果を使って

   \begin{eqnarray}\\ v_2(c^{20} - 1) &=& v_2(c^{10} - 1) + v_2(c^{10} + 1)\\ &=& v_2(c - 1) + v_2(10) + v_2(c^{10} + 1)\\ &=& v_2(c - 1) + v_2(10) + 1\\ &=& v_2(c - 1) + v_2(20)\\ \end{eqnarray}

k=40の場合は、k=20の結果を使って
   \begin{eqnarray}\\ v_2(c^{40} - 1) &=& v_2(c^{20} - 1) + v_2(c^{20} + 1)\\ &=& v_2(c - 1) + v_2(20) + v_2(c^{20} + 1)\\ &=& v_2(c - 1) + v_2(20) + 1\\ &=& v_2(c - 1) + v_2(40)\\ \end{eqnarray}

…とk52の累乗数をかけた数のときは帰納的に示せれそうです。
それどころか、同じ方法で5でなくてもkが奇数m2の累乗数をかけた数のとき、つまり任意の正整数kについていえるんじゃないでしょうか。

やってみましょう。

証明

i=v_2 kとして、k = 2^i m (mは奇数)と表します。
v_2(c^{2^i m} - 1) = v_2(c - 1) + i が成り立つことをiに関する帰納法で示します。

i = 0のとき

   \begin{eqnarray}\\ v_2(c^{2^i m} - 1) &=& v_2(c^m - 1)\\ &=& v_2(c - 1) + v_2(1 + c + \cdots c^{m-1})\\ \end{eqnarray}
であり、1 + c + \cdots c^{m-1}は奇数を奇数個集めた和だからv_2(1 + c + \cdots c^{m-1}) = 0となり、
   v_2(c^{2^i m} - 1) = v_c(c - 1)
より成り立ちます。

i = lのとき成り立つと仮定します。
v_2(c^{2^l m} - 1) = v_2(c - 1) + lです。
cの累乗は4で割った余りが1だから、v_2(c^{2^l m} + 1) = 1

   \begin{eqnarray}\\ v_2(c^{2^{l+1} m} - 1) &=& v_2(c^{2^l m} - 1) + v_2(c^{2^l m} + 1)\\ &=& v_2(c - 1) + l + 1\\ \end{eqnarray}
よりi = l + 1 でも成り立ちます。

よって数学的帰納法により任意の整数i \ge 0についてv_2(c^{2^i m} - 1) = v_2(c - 1) + iが成り立ちます。

これで任意の正整数kについて v_2(c^{k} - 1) = v_2(k) + v_2(c-1) が示されました。
したがって最初の同値変形から任意の偶数n>0についてv_2(D_A(n)) = v_2(n) + v_2(A+1) - 1が成り立ちます。

ラストスパート!

ここでわかっていることを整理しておきましょう。

  • 最大周期となるというのはa_M = 0かつ1 \le n < Mであるすべてのna_n \ne 0であること
  • 任意のn \ge 1についてa_n = 0となる条件はv_2(D_A(n)) \ge s
  • nが奇数ならv_2(D_A(n)) = 0
  • nが偶数ならv_2(D_A(n)) = v_2(n) + v_2(A+1) - 1 (n \ne 0)

v_2(A+1)に具体的な値を入れて考えてみます。Aは奇数よりA+1は偶数だからv_2(A+1) = 0はないですね。よってまずは、1のときです。

v_2(A+1) = 1 のとき

このとき、nが偶数ならv_2(D_A(n)) = v_2(n) + v_2(A+1) - 1 = v_2(n)nが奇数ならv_2(D_A(n)) = 0となります。これはv_2(D_A(n)) = v_2(n)と一つにまとめられますね。

1 \le n < M (= 2^s) であるすべてのnについてn2^sで割り切れないので、v_2(n) < s よって v_2(D_A(n)) < s より a_n \ne 0です。
n = Mのときv_2(n) = sよりa_n = 0です。

よって数列は最大周期となります!

v_2(A+1) = 2 のとき

このときはnが偶数ならv_2(D_A(n)) = v_2(n) + 1nが奇数ならv_2(D_A(n)) = 0です。
1 \le n < 2^sの中にv_2(n) + 1s以上となるような偶数nはあるか…。
ありますね、n = 2^{s-1}です。
よって、v_2(D_A(2^{s-1})) = sとなりますから、a_{2^{s-1}} = 0です。はい、最大周期になりません!

さあ次はv_2(A+1) = 3のとき…とその必要はありません。2のときと同様に3以上のときもv_2(D_A(2^{s-1})) \ge sとなりますから。最大周期になりません。

よってv_2(A+1) = 1ならば最大周期となって、そうでないときは最大周期となりません。
v_2(A+1) = 1は言い換えると(A+1) \bmod 4 = 2、すなわちA \bmod 4 = 1です。

はい、これでひと仕事おしまい。

補足
ここでs \ge 2としておいた布石がきいてます^^; s=1だと2^{s-1}は偶数になりませんからね。

結論

線形合同法数列
   a_0 = F \bmod M
   a_{n+1} = (A a_n + B) \bmod M
ただしM = 2^s (s2以上の整数)

が最大周期となる条件は
   「Bが奇数かつA \bmod 4=1
です。
Aが奇数であることはA \bmod 4=1からいえるので省いていいのです)

おわりに

この記事のほとんどは上記のページに負っています。感謝です。

数学をちゃんと学び始める前の頃の自分でも理解できることを目標として平坦になるよう心がけましたがどうなったでしょうね^^;
最大周期となる条件だけでひいひい言っていますが、線形合同法にはまだまだ面白い宝物が隠されているんじゃないかと思っています。

筆者: oupo (連絡先: oupo.nejiki@gmail.com)