# 擬似乱数生成器：線形合同法の周期・格子構造・並列ストリーム

> 擬似乱数生成器を有限状態の機械として定式化し、線形合同法が全周期をもつ条件（Hull–Dobell）、連続する組が少数の超平面に載ること（Marsaglia）、F_2 線形生成器の最大周期と原始多項式、飛び越しと並列ストリームの重なりの確率を証明する。
> https://rikai.mugen-giken.com/computer-science/stochastic-computation/pseudorandom-generators

## 0. この記事の要点

コンピュータが出す「乱数」は、実はサイコロを振っているのではなく、決まった計算を繰り返して数を並べているだけです。そのため、並べ方に癖があると、見た目はばらばらでも計算の結果が偏ってしまいます。この記事では、代表的な作り方で数の列が何回で元に戻るか、数の組を点として描くとどんな模様になるか、複数の計算機で別々の列を使って重ならないかを、証明つきで調べます。読み終えると、乱数生成器の「癖」を式で見積もれるようになります。

- 擬似乱数生成器は**有限の状態集合の上の決定的な写像**です。出力は必ず最終的に周期的になり、周期は状態数を超えません（<Ref to="prop-eventually-periodic" />）。
- 線形合同法 $x_{k+1} = ax_k + c \bmod m$ が全周期 $m$ をもつ条件は、$\gcd(c,m)=1$、$m$ の各素因数 $p$ で $a\equiv1 \pmod p$、$4\mid m$ なら $a\equiv 1\pmod 4$ の 3 つです（<Ref to="thm-hull-dobell" />）。法が $2^e$ のとき下位 $r$ ビットの周期は $2^r$ しかありません（<Ref to="prop-low-bits" />）。
- 全周期の線形合同法の連続する $d$ 個の組は**格子の平行移動**であり（<Ref to="thm-lcg-lattice" />）、高々 $(d!\,m)^{1/d}$ 枚の平行な超平面の上に載ります（<Ref to="thm-marsaglia-planes" />）。RANDU の 3 つ組は 15 枚の平面に載ります。
- $\mathbb{F}_2$ 上の線形生成器 $x_{k+1} = Ax_k$ の周期が最大値 $2^n-1$ になるのは、$A$ の特性多項式が原始多項式のとき、かつそのときに限ります（<Ref to="thm-f2-linear-period" />）。
- 並列計算では、$k$ 歩先の状態を $O(\log k)$ 回の演算で求める**飛び越し**で列を重なりなく分けられます。開始点を無作為に選ぶ場合、$s$ 本の長さ $L$ の列のどれかが重なる確率は $s(s-1)L/P$ 以下です（<Ref to="thm-stream-overlap" />）。

<Aside type="tip" title="この記事で使う事前知識">

- **合同式**: 整数を $m$ で割った余りだけに注目して等式のように扱う記法。$a \equiv b \pmod m$ は「$a - b$ が $m$ の倍数」の意味。→ [合同](/mathematics/number-theory/congruences-and-fermat#def-congruence)
- **フェルマーの小定理・単元**: 素数 $p$ で割り切れない $a$ について $a^{p-1} \equiv 1 \pmod p$ が成り立つ、という性質。→ [合同式とフェルマーの小定理](/mathematics/number-theory/congruences-and-fermat)
- **位数**: $a^k \equiv 1$ となる最小の正の整数 $k$ のこと。→ [位数](/mathematics/number-theory/quadratic-reciprocity#def-order)
- **素数ごとの割り切れる回数（付値）**: $n$ を素数 $p$ で何回割れるか。たとえば $12 = 2^2 \cdot 3$ なので 2 では 2 回。→ [整除と素因数分解](/mathematics/number-theory/elementary-number-theory)
- **既約多項式**: 次数の低い多項式の積に分解できない多項式。素数の多項式版。→ [多項式環と体の拡大](/mathematics/algebra/polynomial-rings-field-extensions)

</Aside>

## 1. 動機

Monte Carlo 法の定理はすべて「独立で、$[0,1)$ 上に一様に分布する確率変数の列 $U_1, U_2, \ldots$」を仮定します。ところが計算機が実際に返すのは、種（seed）から決定的な規則で作った数列です。同じ種からは同じ列が出るので、確率変数ではありません。それでも Monte Carlo 法が使えるのは、良い生成器の出力が、多くの計算にとって独立一様な列と区別がつかないからです。

問題は「区別がつかない」の範囲です。1960 年代に IBM の計算機で広く使われた RANDU という生成器は、連続する 3 つの出力を 3 次元の点として並べると、すべての点がわずか 15 枚の平面の上に載ります。3 次元の積分にこの生成器を使うと、立方体のほとんどの部分に点が 1 つも落ちません。この欠陥は、出力を 1 つずつ眺めても分かりません。

この章では、生成器の**構造から証明できる性質**を扱います。周期はどれだけか、連続する組はどう並ぶか、並列計算で複数の列に分けたとき列どうしが重ならないと言えるか、の 3 点です。これらは「独立一様な列」というモデルが擬似乱数に当てはまる範囲の下限を与えます。逆に、統計的な検定に合格するかどうかは経験的な事実で、定理にはなりません（<Ref to="rem-empirical" />）。

以下、$\mathbb{N} = \{1, 2, \ldots\}$ とし、整数 $x$ を $m$ で割った余りを $x \bmod m \in \{0, 1, \ldots, m-1\}$ と書きます。合同式の記号と性質は[合同式とフェルマーの小定理](/mathematics/number-theory/congruences-and-fermat)（<Ref to="mathematics/number-theory/congruences-and-fermat#def-congruence" />、<Ref to="mathematics/number-theory/congruences-and-fermat#prop-arithmetic" />）のものを使います。

## 2. 生成器と周期

<Definition id="def-prng" title="擬似乱数生成器">
**擬似乱数生成器**とは、有限集合 $S$（状態集合）、写像 $T : S \to S$（遷移）、有限集合 $O$ への写像 $g : S \to O$（出力関数）の組 $(S, T, g)$ のことである。初期状態（種）$s_0 \in S$ から、状態の列 $s_{k+1} = T(s_k)$ と出力の列 $y_k = g(s_k)$（$k = 0, 1, 2, \ldots$）が決まる。
</Definition>

出力集合 $O$ は、たとえば 64 ビットの整数の全体です。$[0,1)$ の実数がほしいときは、出力の整数を $2^{64}$ などで割ります（倍精度の一様乱数の作り方は[次章](/computer-science/stochastic-computation/sampling-from-distributions)で扱います）。

<Definition id="def-period" title="周期と前周期">
列 $(z_k)_{k \ge 0}$ が**最終的に周期的**であるとは、整数 $\mu \ge 0$ と $\lambda \ge 1$ があって、すべての $k \ge \mu$ で $z_{k+\lambda} = z_k$ となることをいう。そのような $\lambda$ の最小値を**周期**、周期 $\lambda$ に対して条件を満たす $\mu$ の最小値を**前周期**という。前周期が $0$ のとき**純周期的**という。
</Definition>

<Proposition id="prop-eventually-periodic" title="有限状態の生成器は最終的に周期的">
生成器 $(S, T, g)$ と任意の種 $s_0$ について次が成り立つ。

1. 状態の列 $(s_k)$ は最終的に周期的で、周期 $\lambda$ と前周期 $\mu$ は $\mu + \lambda \le |S|$ を満たす。
2. 出力の列 $(y_k)$ も最終的に周期的で、その周期は $\lambda$ の約数である。特に $|S|$ 以下である。
3. $T$ が全単射ならば、状態の列も出力の列も純周期的である。
</Proposition>

<Proof of="prop-eventually-periodic">
**1.** $|S| + 1$ 個の状態 $s_0, s_1, \ldots, s_{|S|}$ の中には、鳩の巣原理によって等しいものがあります。$s_i = s_j$（$i < j$）となる組のうち $j$ が最小のものを取り、$\mu = i$、$\lambda = j - i$ と置きます。$j$ の最小性から $s_0, \ldots, s_{j-1}$ は互いに異なるので、$j = \mu + \lambda \le |S|$ です。

$s_{\mu+\lambda} = s_\mu$ の両辺に $T$ を $k - \mu$ 回施すと、すべての $k \ge \mu$ で $s_{k+\lambda} = s_k$ となります。これが最小の周期であることを確かめます。ある $\lambda'$（$1 \le \lambda' < \lambda$）と $\mu'$ ですべての $k \ge \mu'$ について $s_{k+\lambda'} = s_k$ だとします。$\lambda$ の倍数 $q\lambda$ を $\mu + q\lambda \ge \mu'$ となるように取ると、$s_\mu = s_{\mu+q\lambda} = s_{\mu+q\lambda+\lambda'} = s_{\mu+\lambda'}$ です（最初と最後の等号は周期 $\lambda$ から）。$\mu < \mu + \lambda' < \mu + \lambda = j$ なので、$s_0, \ldots, s_{j-1}$ が互いに異なることに反します。前周期が $\mu$ であることも同じ理由によります（$\mu' < \mu$ で $s_{\mu'+\lambda} = s_{\mu'}$ なら、$\mu' + \lambda < j$ で等しい状態が現れることになります）。

**2.** $k \ge \mu$ で $y_{k+\lambda} = g(s_{k+\lambda}) = g(s_k) = y_k$ なので、$(y_k)$ は最終的に周期的です。その周期を $\lambda_y$ とします。$\lambda = q\lambda_y + r$（$0 \le r < \lambda_y$）と割ると、十分大きい $k$ で $y_{k+r} = y_{k+r+q\lambda_y} = y_{k+\lambda} = y_k$ となります。$r \ge 1$ なら $r$ が $\lambda_y$ より小さい周期を与えて最小性に反するので、$r = 0$、すなわち $\lambda_y \mid \lambda$ です。

**3.** $T$ が全単射で $\mu \ge 1$ だとすると、$T(s_{\mu-1}) = s_\mu = s_{\mu+\lambda} = T(s_{\mu+\lambda-1})$ と単射性から $s_{\mu-1} = s_{\mu+\lambda-1}$ となり、$\mu$ の最小性に反します。よって $\mu = 0$ で、出力の列もすべての $k \ge 0$ で $y_{k+\lambda} = y_k$ を満たします。
</Proof>

この命題は、生成器の設計の出発点です。状態が $b$ ビットなら周期は $2^b$ 以下なので、$2^{64}$ 個を超える乱数を使う計算に 64 ビットの状態の生成器は使えません。また、遷移を全単射にしておけば、前周期の間だけ現れて二度と現れない状態がなくなり、すべての状態がいずれかの閉じた軌道（サイクル）に属します。

## 3. 線形合同法の周期

<Definition id="def-lcg" title="線形合同法・乗算型合同法">
整数 $m \ge 2$（法）、$a$（乗数）、$c$（増分）に対し、$S = O = \{0, 1, \ldots, m-1\}$、
$$
T(x) = (ax + c) \bmod m, \qquad g(x) = x
$$
で定まる生成器を**線形合同法**という。$c \equiv 0 \pmod m$ のものを**乗算型合同法**という。ある種について状態の列の周期が $m$ に等しいとき、**全周期**をもつという。
</Definition>

列 $x_{k+1} = T(x_k)$ は $a$ と $c$ の $m$ を法とする類だけで決まるので、以下では $a \ge 1$、$c \ge 0$ と取ります（必要なら $m$ の倍数を足します）。全周期をもつなら、<Ref to="prop-eventually-periodic" /> の 1 から $\mu = 0$ で、$x_0, \ldots, x_{m-1}$ は $\{0, \ldots, m-1\}$ のすべてを 1 回ずつ取ります。したがって $T$ は全射、有限集合の上なので全単射であり、どの種から出発しても同じ長さ $m$ のサイクルを回ります。全周期は種の取り方によらない性質です。

### 3.1. 閉じた形と等比和の付値

整数 $a \ge 1$ と $k \ge 0$ に対して、等比和を
$$
S_k(a) = 1 + a + a^2 + \cdots + a^{k-1} \quad (k \ge 1), \qquad S_0(a) = 0
$$
と書きます。$(a - 1)S_k(a) = a^k - 1$ です。

<Lemma id="lem-closed-form" title="線形合同法の閉じた形">
整数の列を $X_0 = x_0$、$X_{k+1} = aX_k + c$ で定めると（余りを取らない）、すべての $k \ge 0$ で
$$
X_k = a^kx_0 + c\,S_k(a), \qquad X_k - x_0 = S_k(a)\bigl((a-1)x_0 + c\bigr)
$$
が成り立つ。線形合同法の状態は $x_k \equiv X_k \pmod m$ を満たす。
</Lemma>

<Proof of="lem-closed-form">
$k$ についての帰納法で示します。$k = 0$ では $X_0 = x_0 = a^0x_0 + c\cdot 0$ です。$k$ で成り立つとすると、
$$
X_{k+1} = a\bigl(a^kx_0 + c\,S_k(a)\bigr) + c = a^{k+1}x_0 + c\bigl(aS_k(a) + 1\bigr) = a^{k+1}x_0 + c\,S_{k+1}(a)
$$
です。最後に $aS_k(a) + 1 = 1 + a + \cdots + a^k = S_{k+1}(a)$ を使いました。第 2 の式は $a^k - 1 = (a-1)S_k(a)$ から
$$
X_k - x_0 = (a^k - 1)x_0 + c\,S_k(a) = (a-1)S_k(a)x_0 + c\,S_k(a)
$$
となって従います。$x_k \equiv X_k$ は、合同式の和と積（<Ref to="mathematics/number-theory/congruences-and-fermat#prop-arithmetic" />）を使った帰納法です。$x_{k+1} \equiv ax_k + c \equiv aX_k + c = X_{k+1} \pmod m$。
</Proof>

「$a - 1$ で割る」形 $c(a^k-1)/(a-1)$ で書かれることが多い式ですが、$a - 1$ は $m$ と互いに素とは限らないので、合同式の中で割り算はできません。整数 $S_k(a)$ のまま扱うのが安全です。周期の問題は「$S_k(a)$ がいつ $m$ で割り切れるか」に帰着し、それに答えるのが次の補題です。素数 $p$ についての付値 $v_p$（<Ref to="mathematics/number-theory/elementary-number-theory#def-valuation" />）を使います。

次の補題のねらいは、「$S_k(a)$ が $m$ で割り切れる」という条件を、素数 $p$ ごとの「何回割れるか」に分けて調べられるようにすることです。$v_p(n)$ は $n$ を素数 $p$ で割れる回数で、たとえば $v_2(12) = 2$、$v_3(12) = 1$ です。補題は、条件のもとで $S_k(a)$ の割れる回数が $k$ の割れる回数とちょうど等しいと言っています。証明は、$k$ が $p$ の倍数でない場合、$k = p$ の場合、一般の $k$ の場合の順に進みます。

<Lemma id="lem-geometric-valuation" title="等比和の p 進付値">
$p$ を素数、$a \ge 1$ を整数とし、$a \equiv 1 \pmod p$ とする。$p = 2$ のときはさらに $a \equiv 1 \pmod 4$ とする。このとき、すべての整数 $k \ge 1$ について
$$
v_p\bigl(S_k(a)\bigr) = v_p(k)
$$
が成り立つ。
</Lemma>

<Proof of="lem-geometric-valuation">
3 段階に分けます。仮定を満たす $a$ の全体を $\mathcal{A}_p$ と書きます。$a \in \mathcal{A}_p$ なら、合同式の積（<Ref to="mathematics/number-theory/congruences-and-fermat#prop-arithmetic" />）から任意の $j \ge 1$ で $a^j \in \mathcal{A}_p$ です。

**(i) $p \nmid k$ のとき。** $a^i \equiv 1 \pmod p$ なので $S_k(a) \equiv k \not\equiv 0 \pmod p$。よって $v_p(S_k(a)) = 0 = v_p(k)$ です。

**(ii) $k = p$ のとき。** $p = 2$ なら $S_2(a) = 1 + a \equiv 2 \pmod 4$ なので、$v_2(S_2(a)) = 1$ です。$p$ が奇素数のとき、$a = 1 + pb$（$b \ge 0$ は整数）と書くと、二項定理から $a^i = (1 + pb)^i \equiv 1 + ipb \pmod{p^2}$ です（$p^2$ 以上の冪を含む項を落としました）。したがって
$$
S_p(a) = \sum_{i=0}^{p-1} a^i \equiv p + pb\sum_{i=0}^{p-1} i = p + pb\cdot\frac{p(p-1)}{2} \pmod{p^2}.
$$
$p$ は奇数なので $(p-1)/2$ は整数で、第 2 項は $p^2$ の倍数です。よって $S_p(a) \equiv p \pmod{p^2}$ となり、$v_p(S_p(a)) = 1$ です。

**(iii) 一般の $k$。** 正の整数 $j, k$ について
$$
S_{jk}(a) = \sum_{l=0}^{j-1}\sum_{i=0}^{k-1} a^{kl+i} = \Bigl(\sum_{l=0}^{j-1} (a^k)^l\Bigr)\Bigl(\sum_{i=0}^{k-1} a^i\Bigr) = S_j(a^k)\,S_k(a)
$$
です。まず $e \ge 0$ についての帰納法で $v_p(S_{p^e}(a)) = e$ をすべての $a \in \mathcal{A}_p$ について示します。$e = 0$ では $S_1(a) = 1$ です。$e \ge 1$ では $S_{p^e}(a) = S_p(a^{p^{e-1}})\,S_{p^{e-1}}(a)$ で、$a^{p^{e-1}} \in \mathcal{A}_p$ なので (ii) から第 1 因子の付値は $1$、帰納法の仮定から第 2 因子の付値は $e - 1$ です。付値は積について加法的なので（<Ref to="mathematics/number-theory/elementary-number-theory#prop-valuation-rules" /> の 1）、$v_p(S_{p^e}(a)) = e$ です。

一般の $k = p^ek'$（$p \nmid k'$）では $S_k(a) = S_{k'}(a^{p^e})\,S_{p^e}(a)$ で、(i) から第 1 因子の付値は $0$、第 2 因子の付値は $e$ です。よって $v_p(S_k(a)) = e = v_p(k)$ です。
</Proof>

$p = 2$ で $a \equiv 3 \pmod 4$ のときは、$S_2(a) = 1 + a$ が $4$ で割り切れるので補題は成り立ちません。これが、定理の第 3 の条件が現れる理由です。

### 3.2. 全周期の条件

<Theorem id="thm-hull-dobell" title="Hull–Dobell の定理">
$m \ge 2$、$a \ge 1$、$c \ge 0$ を整数とする。線形合同法 $x_{k+1} = (ax_k + c) \bmod m$ が全周期 $m$ をもつための必要十分条件は、次の 3 つがすべて成り立つことである。

1. $\gcd(c, m) = 1$。
2. $m$ を割り切るすべての素数 $p$ について $a \equiv 1 \pmod p$。
3. $4 \mid m$ ならば $a \equiv 1 \pmod 4$。
</Theorem>

<Proof of="thm-hull-dobell">
**十分性。** 条件 1〜3 を仮定し、任意の $x_0$ を取ります。$w = (a-1)x_0 + c$ と置くと、<Ref to="lem-closed-form" /> から $x_k - x_0 \equiv S_k(a)\,w \pmod m$ です。

まず $\gcd(w, m) = 1$ を示します。$m$ の素因数 $p$ について、条件 2 から $a - 1 \equiv 0$ なので $w \equiv c \pmod p$ です。条件 1 から $p \nmid c$ なので $p \nmid w$ です。$m$ のどの素因数も $w$ を割らないので、$\gcd(w, m) = 1$ です。したがって、互いに素な因子を落とす補題（<Ref to="mathematics/number-theory/elementary-number-theory#lem-euclid-general" />）から
$$
x_k \equiv x_0 \pmod m \iff m \mid S_k(a)\,w \iff m \mid S_k(a)
$$
です。

次に、$k \ge 1$ について $m \mid S_k(a) \iff m \mid k$ を示します。整除は素数ごとの付値の比較と同値です（<Ref to="mathematics/number-theory/elementary-number-theory#prop-valuation-rules" /> の 2）。これは、$m$ を素数冪の積に分け、各素数冪を法として別々に調べてよいということで、中国剰余定理（<Ref to="mathematics/number-theory/congruences-and-fermat#thm-crt" />）の整除についての言い換えです。$m$ の素因数 $p$ ごとに、$v_p(m) \le v_p(S_k(a)) \iff v_p(m) \le v_p(k)$ を確かめます。

- $p$ が奇素数のとき、および $p = 2$ で $4 \mid m$ のとき: 条件 2・3 から <Ref to="lem-geometric-valuation" /> が使えて $v_p(S_k(a)) = v_p(k)$ なので、両辺は同じ条件です。
- $p = 2$ で $v_2(m) = 1$ のとき: 条件 2 から $a$ は奇数なので $S_k(a) \equiv k \pmod 2$ です。よって $v_2(S_k(a)) \ge 1 \iff v_2(k) \ge 1$ です。

以上から $x_k \equiv x_0 \pmod m \iff m \mid k$ が任意の種 $x_0$ について成り立ちます。$0 \le i < j \le m - 1$ で $x_i = x_j$ だとすると、$x_i$ を種とみて $x_{i + (j-i)} = x_i$ から $m \mid j - i$ となりますが、$1 \le j - i \le m - 1$ なので矛盾です。ゆえに $x_0, \ldots, x_{m-1}$ は互いに異なり、$x_m = x_0$ なので、周期はちょうど $m$ です。

**必要性。** ある種で周期が $m$ だとします。定義の後で述べたとおり、$T$ は全単射で、列 $(x_k)$ は純周期的に $\{0, \ldots, m-1\}$ のすべての値を取り続けます。

*条件 1.* $x_{k_0} = 0$ となる $k_0$ を取ると、<Ref to="lem-closed-form" /> を種 $0$ に使って $x_{k_0+j} \equiv c\,S_j(a) \pmod m$ です。$g = \gcd(c, m)$ は右辺と $m$ を割るので、$j \ge 0$ で現れる値はすべて $g$ の倍数です。この列は値 $1$ も取るので $g \mid 1$、すなわち $g = 1$ です。

*条件 2.* $p$ を $m$ の素因数とし、$y_k = x_k \bmod p$ と置きます。$p \mid m$ なので $y_{k+1} \equiv ay_k + c \pmod p$ です。$a \not\equiv 1 \pmod p$ だと仮定します。$\mathbb{Z}/p\mathbb{Z}$ は体なので（<Ref to="mathematics/number-theory/congruences-and-fermat#cor-field" />）、$(a - 1)u \equiv 1 \pmod p$ となる整数 $u$ があります。$\xi = -cu$ と置くと $a\xi + c = \xi + (a-1)\xi + c \equiv \xi - c + c = \xi \pmod p$ で、$\xi$ は法 $p$ での不動点です。$(x_k)$ はすべての値を取るので $y_{k_1} \equiv \xi$ となる $k_1$ があり、帰納的にすべての $k \ge k_1$ で $y_k \equiv \xi$ です。一方 $(x_k)_{k \ge k_1}$ は $0$ と $1$ の両方を取るので、$(y_k)_{k \ge k_1}$ は $p$ を法として異なる 2 つの値を取ります。これは矛盾です。

*条件 3.* $4 \mid m$ とし、$z_k = x_k \bmod 4$ と置きます。条件 2（$p = 2$）から $a$ は奇数です。$a \equiv 3 \pmod 4$ だと仮定すると、$a^2 \equiv 1$、$a + 1 \equiv 0 \pmod 4$ なので
$$
z_{k+2} \equiv a^2z_k + c(a + 1) \equiv z_k \pmod 4
$$
です。すると $(z_k)$ が取る値は $z_0, z_1$ の高々 2 つですが、$(x_k)$ はすべての値を取るので $(z_k)$ は $4$ つの値を取るはずです。矛盾なので $a \equiv 1 \pmod 4$ です。
</Proof>

定理が述べているのは周期だけです。$a = 1$、$c = 1$ は 3 条件を満たし全周期をもちますが、出力は $0, 1, 2, \ldots$ で、乱数としては使えません。全周期は必要条件の 1 つにすぎず、組の並び方（§4）を別に調べる必要があります。

<Example id="ex-small-moduli" title="法 16 と法 31 の全周期">
$m = 16$ の素因数は $2$ だけで $4 \mid 16$ なので、条件は「$c$ が奇数、$a \equiv 1 \pmod 4$」です。$a \in \{1, 5, 9, 13\}$ と $c \in \{1, 3, \ldots, 15\}$ の $4 \times 8 = 32$ 組が全周期をもちます。$a = 5$、$c = 3$、$x_0 = 0$ では
$$
0, 3, 2, 13, 4, 7, 6, 1, 8, 11, 10, 5, 12, 15, 14, 9, 0, \ldots
$$
と 16 個の値がすべて現れます。$a = 3$（$a \equiv 3 \pmod 4$）では条件 3 が破れ、$c = 3$ のとき $0, 3, 12, 7, 8, 11, 4, 15, 0$ と周期は $8$ です。

$m = 31$（素数）では条件 2 が $a \equiv 1 \pmod{31}$ を要求するので、全周期をもつのは $a = 1$、$c \ne 0$ の 30 組だけで、どれも等差数列です。素数の法では、全周期の線形合同法は役に立たず、次の乗算型合同法（周期 $p - 1$）を使います。
</Example>

### 3.3. 乗算型合同法と下位ビット

$c = 0$ のとき、<Ref to="thm-hull-dobell" /> の条件 1 は成り立たず、周期は $m$ より小さくなります（$0$ が不動点だからです）。このときの最大の周期を、よく使われる 2 つの法について求めます。整数 $a$ の $(\mathbb{Z}/m\mathbb{Z})^\times$（<Ref to="mathematics/number-theory/congruences-and-fermat#def-units" />）での位数（<Ref to="mathematics/algebra/groups#def-order" />）を $\mathrm{ord}_m(a)$ と書きます。

次の命題は、$c = 0$ の場合に周期がどこまで長くなれるかを答えます。素数の法では $p - 1$ まで伸ばせ、法が $2^e$ では 4 分の 1 までしか伸ばせない、というのが結論です。2 番目の主張の証明は、$a$ を 4 で割った余りが 1 の場合と 3 の場合に分け、3 の場合は $a$ を 2 乗して 1 の場合に戻します。

<Proposition id="prop-multiplicative-period" title="乗算型合同法の周期">
乗算型合同法 $x_{k+1} = ax_k \bmod m$ について次が成り立つ。

1. $\gcd(a, m) = \gcd(x_0, m) = 1$ ならば、列 $(x_k)$ は純周期的で周期は $\mathrm{ord}_m(a)$ である。
2. $m = 2^e$（$e \ge 3$）で $a$ が奇数のとき、$\mathrm{ord}_{2^e}(a) \le 2^{e-2}$ であり、$a \equiv 3$ または $5 \pmod 8$ ならば等号が成り立つ。$e \ge 4$ のときは、等号が成り立つのはこの場合に限る。
3. $m = p$ が素数で $p \nmid a$ のとき、$\mathrm{ord}_p(a)$ は $p - 1$ の約数であり、$p - 1$ に等しいのは $a$ が $p$ を法とする原始根のとき、かつそのときに限る。原始根は必ず存在する。
</Proposition>

<Proof of="prop-multiplicative-period">
**1.** $x_k \equiv a^kx_0 \pmod m$ です（<Ref to="lem-closed-form" /> で $c = 0$）。$\gcd(a, m) = 1$ なので $x \mapsto ax$ は $\mathbb{Z}/m\mathbb{Z}$ の全単射で（<Ref to="mathematics/number-theory/congruences-and-fermat#prop-unit-criterion" /> により $a$ は単元）、<Ref to="prop-eventually-periodic" /> の 3 から列は純周期的です。$\gcd(x_0, m) = 1$ なので消去法則（<Ref to="mathematics/number-theory/congruences-and-fermat#prop-cancellation" />）が使えて、$a^kx_0 \equiv x_0 \iff a^k \equiv 1 \pmod m$ です。これを満たす最小の $k \ge 1$ が $\mathrm{ord}_m(a)$ です。

**2.** 位数は $a$ の $2^e$ を法とする類で決まるので、$1 \le a < 2^e$ とします。$a = 1$ の位数は $1$ なので $a \ge 3$ とし、$a$ を $4$ で割った余りで場合を分けます。

*$a \equiv 1 \pmod 4$ のとき。* <Ref to="lem-geometric-valuation" />（$p = 2$）と $a^k - 1 = (a-1)S_k(a)$、付値の加法性（<Ref to="mathematics/number-theory/elementary-number-theory#prop-valuation-rules" /> の 1）から
$$
v_2(a^k - 1) = v_2(a - 1) + v_2(k) \qquad (k \ge 1)
$$
です。$t = v_2(a-1)$ と置くと、$0 < a - 1 < 2^e$ から $2 \le t < e$ です。$a^k \equiv 1 \pmod{2^e}$ は $t + v_2(k) \ge e$ と同値で、これを満たす最小の $k \ge 1$ は $2^{e-t}$ です。すなわち
$$
\mathrm{ord}_{2^e}(a) = 2^{e - v_2(a-1)} \qquad (a \equiv 1 \pmod 4,\ 1 < a < 2^e).
$$
$a \equiv 5 \pmod 8$ なら $t = 2$ で位数は $2^{e-2}$、$a \equiv 1 \pmod 8$ なら $t \ge 3$ で位数は $2^{e-3}$ 以下です。

*$a \equiv 3 \pmod 4$ のとき。* $k$ が奇数なら $a^k \equiv 3 \pmod 4$ なので $a^k \not\equiv 1 \pmod{2^e}$ です。よって $a^k \equiv 1$ となる $k$ は偶数で、$k = 2j$ と書くと条件は $(a^2)^j \equiv 1$ です。したがって $\mathrm{ord}_{2^e}(a) = 2\,\mathrm{ord}_{2^e}(a^2)$ です。$b = a^2$ は $b \equiv 1 \pmod 4$、$b > 1$ を満たします。$t' = v_2(b - 1) = v_2(a-1) + v_2(a+1) = 1 + v_2(a+1)$ と置きます。$b$ を $2^e$ で割った余りを $b'$ とすると、$t' \ge e$ なら $b' = 1$ で $b$ の位数は $1$、$t' < e$ なら $v_2(b' - 1) = t'$ なので前の場合から $b$ の位数は $2^{e-t'}$ です。まとめると
$$
\mathrm{ord}_{2^e}(a) = 2\cdot 2^{\max(e - t',\,0)}, \qquad t' = 1 + v_2(a+1).
$$
$a \equiv 3 \pmod 8$ なら $v_2(a+1) = 2$、$t' = 3 \le e$ で、位数は $2\cdot 2^{e-3} = 2^{e-2}$ です。$a \equiv 7 \pmod 8$ なら $v_2(a+1) \ge 3$、$t' \ge 4$ で、位数は $2\cdot 2^{\max(e-4, 0)}$ 以下です。これは $e \ge 4$ なら $2^{e-3} < 2^{e-2}$、$e = 3$ なら $2 = 2^{e-2}$ です。

以上で、すべての場合に位数は $2^{e-2}$ 以下、$a \equiv 3, 5 \pmod 8$ で等号、$e \ge 4$ では $a \equiv 1, 7 \pmod 8$ のとき $2^{e-3}$ 以下であることが分かりました。

**3.** フェルマーの小定理（<Ref to="mathematics/number-theory/congruences-and-fermat#thm-fermat-little" />）から $a^{p-1} \equiv 1 \pmod p$ なので、位数による冪の判定（<Ref to="mathematics/algebra/groups#thm-order-divides" />）から $\mathrm{ord}_p(a) \mid p - 1$ です。$\varphi(p) = p - 1$ なので、$\mathrm{ord}_p(a) = p - 1$ は $a$ が原始根であることの定義そのものです（<Ref to="mathematics/number-theory/quadratic-reciprocity#def-primitive-root" />）。原始根の存在は <Ref to="mathematics/number-theory/quadratic-reciprocity#thm-primitive-root-prime" /> です。
</Proof>

法 $2^e$ の乗算型合同法は、種を奇数に限っても状態の 4 分の 1 しか回りません。法が素数 $p$ なら $0$ 以外の全状態を回れますが、$p$ で割った余りを取る演算は $2^e$ で割るより手間がかかります。$p = 2^{31} - 1$（Mersenne 素数）と $a = 16807 = 7^5$ の組は「最小標準」の生成器として長く使われました。

<Proposition id="prop-low-bits" title="法 2^e の線形合同法の下位ビット">
$m = 2^e$ の線形合同法の状態 $x_k$ の下位 $r$ ビット $x_k \bmod 2^r$（$1 \le r \le e$）の列は、最終的に周期的で、周期は $2^r$ 以下である。全周期をもつ場合は、純周期的で周期がちょうど $2^r$ である。特に最下位ビットは $0, 1, 0, 1, \ldots$ と交互に現れる。
</Proposition>

<Proof of="prop-low-bits">
$z_k = x_k \bmod 2^r$ と置きます。$2^r \mid 2^e$ なので、$x_{k+1} \equiv ax_k + c \pmod{2^e}$ から $z_{k+1} \equiv az_k + c \pmod{2^r}$ です。つまり $(z_k)$ は法 $2^r$、同じ $a, c$ の線形合同法の状態の列です。状態数は $2^r$ なので、<Ref to="prop-eventually-periodic" /> の 1 から最終的に周期的で、周期は $2^r$ 以下です。

法 $2^e$ で全周期をもつとき、<Ref to="thm-hull-dobell" /> から $c$ は奇数、$a$ は奇数で、$e \ge 2$ なら $a \equiv 1 \pmod 4$ です。これらは法 $2^r$ についての同じ定理の 3 条件を含むので（$r = 1$ のときは条件 3 は不要）、$(z_k)$ は法 $2^r$ で全周期 $2^r$ をもち、純周期的です。$r = 1$ では周期 $2$ で、$0$ と $1$ が交互に現れます。
</Proof>

<Ref to="ex-small-moduli" /> の $a = 5$、$c = 3$、$m = 16$ の列で確かめると、下位 2 ビットは $0, 3, 2, 1, 0, 3, 2, 1, \ldots$ と周期 $4$ です。64 ビットの状態をもつ線形合同法でも、最下位ビットは交互に並び、下位 8 ビットは 256 個ごとに繰り返します。「乱数を 2 で割った余りで表裏を決める」という使い方は、この種の生成器では正しく働きません。対策は 2 つあります。上位ビットだけを出力に使うこと、そして §6 で述べるように出力関数 $g$ で状態をかき混ぜることです。

次のコードは、小さな法について全周期をもつ $(a, c)$ を全数で数え、定理の判定と照合します。

```python
# Python 3.9.6（標準ライブラリのみ）
from math import gcd

def period(a, c, m, x0=0):
    seen, x, k = {}, x0, 0
    while x not in seen:
        seen[x] = k
        x = (a * x + c) % m
        k += 1
    return k - seen[x]          # 周期（前周期を除く）

def prime_factors(m):
    ps, p = [], 2
    while p * p <= m:
        if m % p == 0:
            ps.append(p)
            while m % p == 0:
                m //= p
        p += 1
    return ps + ([m] if m > 1 else [])

def hull_dobell(a, c, m):
    return (gcd(c, m) == 1
            and all((a - 1) % p == 0 for p in prime_factors(m))
            and (m % 4 != 0 or (a - 1) % 4 == 0))

for m in (16, 31, 36, 100, 360):
    full = agree = 0
    for a in range(m):
        for c in range(m):
            is_full = (period(a, c, m) == m)
            full += is_full
            agree += (is_full == hull_dobell(a, c, m))
    print(f"m={m:4d}  全周期の (a,c) = {full:5d} 組 / {m*m:6d}   判定と一致 = {agree:6d}")
```

```text
m=  16  全周期の (a,c) =    32 組 /    256   判定と一致 =    256
m=  31  全周期の (a,c) =    30 組 /    961   判定と一致 =    961
m=  36  全周期の (a,c) =    36 組 /   1296   判定と一致 =   1296
m= 100  全周期の (a,c) =   200 組 /  10000   判定と一致 =  10000
m= 360  全周期の (a,c) =   576 組 / 129600   判定と一致 = 129600
```

5 つの法のすべての $(a, c)$ で、実際の周期による判定と定理の 3 条件による判定が一致しています。$m = 16$ の 32 組、$m = 31$ の 30 組は <Ref to="ex-small-moduli" /> の数え上げと同じです。

## 4. 格子構造

周期が長くても、連続する出力の間に強い規則があれば乱数としては使えません。線形合同法では、その規則が格子として完全に記述できます。

<Definition id="def-equidistribution" title="k 次元の均等分布">
出力が $v$ ビットの整数（$O = \{0, \ldots, 2^v - 1\}$）で、出力の列が純周期的で周期 $P$ の生成器を考える。1 周期にわたる連続する $k$ 個の組 $(y_i, y_{i+1}, \ldots, y_{i+k-1})$（$i = 0, \ldots, P-1$、添字は周期的に延長）の中に、$O^k$ のすべての元が同じ回数ずつ現れるとき、この生成器は（$v$ ビット精度で）**$k$ 次元均等分布**であるという。全零の組だけ 1 回少ないことを許す流儀もある。
</Definition>

全周期の線形合同法は、$O = \{0, \ldots, m-1\}$ として 1 次元では均等分布です（各値がちょうど 1 回現れます）。しかし 2 次元では、$m^2$ 通りの組のうち $m$ 通りしか現れないので、均等分布ではありません。現れる $m$ 個の組がどう並ぶかを述べるのが次の定理です。$\mathbb{R}^d$ の格子と余体積の定義は <Ref to="mathematics/number-theory/class-group-and-minkowski-bound#def-lattice" /> に従います。ベクトルの成分の添字は $0$ から $d - 1$ とします。

格子とは、いくつかのベクトルの整数倍の和で表される点の集まりで、平面なら方眼紙の交点を斜めにゆがめたような等間隔の点の並びです。次の定理は、出力の組が必ずそのような並びの上にあることを示します。並びが格子なら、点をすべて通る等間隔の平行線（3 次元なら平面）の族が取れます。

<Theorem id="thm-lcg-lattice" title="線形合同法の格子構造">
線形合同法 $x_{k+1} = (ax_k + c) \bmod m$ と $d \ge 1$ について、
$$
\boldsymbol{u}_k = \frac{1}{m}(x_k, x_{k+1}, \ldots, x_{k+d-1}) \in [0,1)^d, \qquad
\boldsymbol{g} = \frac{1}{m}(1, a, \ldots, a^{d-1}), \qquad
\boldsymbol{v}_0 = \frac{c}{m}\bigl(S_0(a), S_1(a), \ldots, S_{d-1}(a)\bigr)
$$
と置き、$L = \mathbb{Z}\boldsymbol{g} + \mathbb{Z}^d$ とする。

1. $L$ は余体積 $1/m$ の格子で、すべての $k \ge 0$ について $\boldsymbol{u}_k \in (\boldsymbol{v}_0 + L) \cap [0,1)^d$ である。
2. 全周期をもつとき、$\{\boldsymbol{u}_k : 0 \le k < m\} = (\boldsymbol{v}_0 + L) \cap [0,1)^d$ である。
</Theorem>

<Proof of="thm-lcg-lattice">
標準基底を $\boldsymbol{e}_0, \ldots, \boldsymbol{e}_{d-1}$ とします。

**$L$ が格子であること。** $m\boldsymbol{g} = \boldsymbol{e}_0 + \sum_{i \ge 1} a^i\boldsymbol{e}_i$ なので $\boldsymbol{e}_0 = m\boldsymbol{g} - \sum_{i\ge1} a^i\boldsymbol{e}_i$ です。よって $L$ は $\boldsymbol{g}, \boldsymbol{e}_1, \ldots, \boldsymbol{e}_{d-1}$ の整数係数一次結合の全体に等しく、これらを列とする行列は下三角で対角成分が $1/m, 1, \ldots, 1$ なので、行列式は $1/m \ne 0$ です。したがって一次独立で、$L$ は余体積 $1/m$ の格子です（<Ref to="mathematics/number-theory/class-group-and-minkowski-bound#lem-fundamental-domain" /> の 3）。

**1.** $x = x_k$ と置きます。<Ref to="lem-closed-form" /> を種 $x$ に使うと、$j = 0, \ldots, d-1$ について $x_{k+j} \equiv a^jx + cS_j(a) \pmod m$ なので、整数 $n_j$ があって $x_{k+j} = a^jx + cS_j(a) + mn_j$ です。$m$ で割って並べると
$$
\boldsymbol{u}_k = x\boldsymbol{g} + \boldsymbol{v}_0 + (n_0, \ldots, n_{d-1}) \in \boldsymbol{v}_0 + L
$$
です。各成分は $x_{k+j}/m \in [0, 1)$ です。

**2.** 1 から「$\subseteq$」は済んでいます。逆を示します。$\boldsymbol{w} \in (\boldsymbol{v}_0 + L) \cap [0,1)^d$ を取り、$\boldsymbol{w} = \boldsymbol{v}_0 + x\boldsymbol{g} + \boldsymbol{n}$（$x \in \mathbb{Z}$、$\boldsymbol{n} \in \mathbb{Z}^d$）と書きます。$x = x' + mq$（$0 \le x' < m$）と割ると、$mq\boldsymbol{g} = q(1, a, \ldots, a^{d-1}) \in \mathbb{Z}^d$ なので、$\boldsymbol{n}$ を取り替えて $0 \le x < m$ としてよいです。全周期なので $x_k = x$ となる $k$（$0 \le k < m$）があります。1 の計算から $\boldsymbol{u}_k - \boldsymbol{w} \in \mathbb{Z}^d$ で、$\boldsymbol{u}_k$ も $\boldsymbol{w}$ も $[0,1)^d$ に属するので、差の各成分は絶対値が $1$ 未満の整数、すなわち $0$ です。よって $\boldsymbol{w} = \boldsymbol{u}_k$ です。
</Proof>

格子の点は、平行で等間隔な超平面の族の上に並びます。その超平面の枚数が少ないと、単位立方体の中に点のない厚い層ができます。次の定理は、どんなに乗数をうまく選んでも、枚数は $m^{1/d}$ の程度でしかないことを示します。

<Theorem id="thm-marsaglia-planes" title="Marsaglia の定理">
線形合同法 $x_{k+1} = (ax_k + c) \bmod m$ と $d \ge 1$ について、次を満たす零でない整数ベクトル $\boldsymbol{c} = (c_0, \ldots, c_{d-1})$ が存在する。
$$
\sum_{i=0}^{d-1} c_ia^i \equiv 0 \pmod m, \qquad \ell := \sum_{i=0}^{d-1} |c_i| \le (d!\,m)^{1/d}.
$$
このような $\boldsymbol{c}$ について、実数 $\theta$ があって、すべての $k$ で $\langle \boldsymbol{c}, \boldsymbol{u}_k\rangle \in \theta + \mathbb{Z}$ となる。したがって点 $\boldsymbol{u}_k$ はすべて、平行な超平面 $\{\boldsymbol{u} : \langle\boldsymbol{c}, \boldsymbol{u}\rangle = \theta + j\}$（$j \in \mathbb{Z}$）のうち高々 $\ell$ 枚の上に載る。特に、枚数は $(d!\,m)^{1/d}$ 以下である。
</Theorem>

<Proof of="thm-marsaglia-planes">
**$\boldsymbol{c}$ の存在。** $\Lambda = \{\boldsymbol{c} \in \mathbb{Z}^d : \sum_i c_ia^i \equiv 0 \pmod m\}$ と置きます。$\boldsymbol{c} \in \Lambda$ は「$c_0 = -\sum_{i \ge 1} c_ia^i + mq$ となる整数 $q$ がある」ことと同値なので、
$$
\boldsymbol{w}_0 = m\boldsymbol{e}_0, \qquad \boldsymbol{w}_i = -a^i\boldsymbol{e}_0 + \boldsymbol{e}_i \quad (1 \le i \le d-1)
$$
として $\Lambda = \{q\boldsymbol{w}_0 + \sum_{i\ge1} c_i\boldsymbol{w}_i\}$ です。これらを列とする行列は上三角で対角成分が $m, 1, \ldots, 1$ なので、$\Lambda$ は余体積 $m$ の格子です（<Ref to="mathematics/number-theory/class-group-and-minkowski-bound#lem-fundamental-domain" /> の 3）。

$r > 0$ に対して $X_r = \{\boldsymbol{t} \in \mathbb{R}^d : \sum_i |t_i| \le r\}$ と置きます。これはコンパクトな凸集合で $X_r = -X_r$ です（凸性は三角不等式から）。体積を求めます。$X_r$ は座標の符号による $2^d$ 個の合同な部分に分かれ、その 1 つは単体 $\Delta_d(r) = \{t_i \ge 0,\ \sum t_i \le r\}$ です。$\mathrm{vol}(\Delta_d(r)) = r^d/d!$ を $d$ についての帰納法で示します。$d = 1$ では長さ $r$ です。$d \ge 2$ では、最後の座標 $t_{d-1} = t$ を固定した切り口が $\Delta_{d-1}(r - t)$ なので、Fubini の定理（<Ref to="mathematics/real-analysis/product-measures-fubini#thm-fubini" />）から
$$
\mathrm{vol}(\Delta_d(r)) = \int_0^r \frac{(r-t)^{d-1}}{(d-1)!}\,dt = \frac{r^d}{d!}
$$
です。よって $\mathrm{vol}(X_r) = (2r)^d/d!$ です。$r = (d!\,m)^{1/d}$ と取ると $\mathrm{vol}(X_r) = 2^dm = 2^d\,\mathrm{covol}(\Lambda)$ なので、Minkowski の格子点定理のコンパクトな場合（<Ref to="mathematics/number-theory/class-group-and-minkowski-bound#thm-minkowski" />）から、$X_r$ は $\Lambda$ の零でない点 $\boldsymbol{c}$ を含みます。これが求めるベクトルです。

**超平面に載ること。** <Ref to="thm-lcg-lattice" /> の 1 の証明から $\boldsymbol{u}_k = x_k\boldsymbol{g} + \boldsymbol{v}_0 + \boldsymbol{n}$（$\boldsymbol{n} \in \mathbb{Z}^d$）です。$\langle\boldsymbol{c}, \boldsymbol{g}\rangle = \frac1m\sum_i c_ia^i \in \mathbb{Z}$、$\langle\boldsymbol{c}, \boldsymbol{n}\rangle \in \mathbb{Z}$ なので、$\theta = \langle\boldsymbol{c}, \boldsymbol{v}_0\rangle$ と置けば $\langle\boldsymbol{c}, \boldsymbol{u}_k\rangle \in \theta + \mathbb{Z}$ です。

**枚数。** $P = \sum_{c_i > 0} c_i$、$N = \sum_{c_i < 0} |c_i|$ と置くと $P + N = \ell$ です。$\boldsymbol{u} \in [0,1)^d$ について、$c_i > 0$ の項は $0 \le c_iu_i < c_i$、$c_i < 0$ の項は $c_i < c_iu_i \le 0$ を満たすので、$\langle\boldsymbol{c}, \boldsymbol{u}\rangle$ は区間 $I$ に属します。ここで $I$ は、$P, N > 0$ なら $(-N, P)$、$N = 0$ なら $[0, P)$、$P = 0$ なら $(-N, 0]$ です。どの場合も $I$ は長さ $\ell$ の閉区間でない区間なので、$I$ の 2 点の差は $\ell$ より真に小さくなります。$\theta + \mathbb{Z}$ の点が $I$ の中に $\ell + 1$ 個あったとすると、最大と最小の差は $\ell$ 以上となり矛盾します。よって $I$ に入る $\theta + j$ は高々 $\ell$ 個で、点が載りうる超平面は高々 $\ell$ 枚です。
</Proof>

$\boldsymbol{c}$ の存在に使った Minkowski の定理は、基本領域への折り返しに鳩の巣原理を使って証明されます。したがってこの定理は、「$m$ 個しかない剰余に対して候補のベクトルが多すぎる」という数え上げの帰結です。$m = 2^{32}$ のとき上界は $d = 3$ で約 $2953$ 枚、$d = 10$ では約 $41$ 枚です。高次元の組を使う計算に、法の小さい線形合同法を使えない理由がここにあります。

<Figure caption="m = 64、a = 37、c = 1 の線形合同法の連続する組 (x_k, x_{k+1})/64 の 64 点。すべての点が −3u + 7v = 7/64 + j（j = −2, …, 6）の 9 本の平行な直線に載る。関係式は −3 + 7·37 = 256 ≡ 0 (mod 64) から。">
<svg viewBox="0 0 316 316" width="100%" role="img" aria-label="線形合同法の連続する 2 つの出力の組が 9 本の平行な直線の上に並ぶ様子" style="max-width: 420px">
  <rect x="30" y="30" width="256" height="256" fill="none" stroke="currentColor" stroke-width="1" />
  <g stroke="var(--sl-color-accent)" stroke-width="0.8" opacity="0.7">
    <line x1="191.3" y1="286.0" x2="286.0" y2="245.4"/><line x1="106.0" y1="286.0" x2="286.0" y2="208.9"/>
    <line x1="30.0" y1="282.0" x2="286.0" y2="172.3"/><line x1="30.0" y1="245.4" x2="286.0" y2="135.7"/>
    <line x1="30.0" y1="208.9" x2="286.0" y2="99.1"/><line x1="30.0" y1="172.3" x2="286.0" y2="62.6"/>
    <line x1="30.0" y1="135.7" x2="276.7" y2="30.0"/><line x1="30.0" y1="99.1" x2="191.3" y2="30.0"/>
    <line x1="30.0" y1="62.6" x2="106.0" y2="30.0"/>
  </g>
  <g fill="currentColor">
    <circle cx="30" cy="282" r="2.2"/><circle cx="34" cy="134" r="2.2"/><circle cx="38" cy="242" r="2.2"/><circle cx="42" cy="94" r="2.2"/>
    <circle cx="46" cy="202" r="2.2"/><circle cx="50" cy="54" r="2.2"/><circle cx="54" cy="162" r="2.2"/><circle cx="58" cy="270" r="2.2"/>
    <circle cx="62" cy="122" r="2.2"/><circle cx="66" cy="230" r="2.2"/><circle cx="70" cy="82" r="2.2"/><circle cx="74" cy="190" r="2.2"/>
    <circle cx="78" cy="42" r="2.2"/><circle cx="82" cy="150" r="2.2"/><circle cx="86" cy="258" r="2.2"/><circle cx="90" cy="110" r="2.2"/>
    <circle cx="94" cy="218" r="2.2"/><circle cx="98" cy="70" r="2.2"/><circle cx="102" cy="178" r="2.2"/><circle cx="106" cy="286" r="2.2"/>
    <circle cx="110" cy="138" r="2.2"/><circle cx="114" cy="246" r="2.2"/><circle cx="118" cy="98" r="2.2"/><circle cx="122" cy="206" r="2.2"/>
    <circle cx="126" cy="58" r="2.2"/><circle cx="130" cy="166" r="2.2"/><circle cx="134" cy="274" r="2.2"/><circle cx="138" cy="126" r="2.2"/>
    <circle cx="142" cy="234" r="2.2"/><circle cx="146" cy="86" r="2.2"/><circle cx="150" cy="194" r="2.2"/><circle cx="154" cy="46" r="2.2"/>
    <circle cx="158" cy="154" r="2.2"/><circle cx="162" cy="262" r="2.2"/><circle cx="166" cy="114" r="2.2"/><circle cx="170" cy="222" r="2.2"/>
    <circle cx="174" cy="74" r="2.2"/><circle cx="178" cy="182" r="2.2"/><circle cx="182" cy="34" r="2.2"/><circle cx="186" cy="142" r="2.2"/>
    <circle cx="190" cy="250" r="2.2"/><circle cx="194" cy="102" r="2.2"/><circle cx="198" cy="210" r="2.2"/><circle cx="202" cy="62" r="2.2"/>
    <circle cx="206" cy="170" r="2.2"/><circle cx="210" cy="278" r="2.2"/><circle cx="214" cy="130" r="2.2"/><circle cx="218" cy="238" r="2.2"/>
    <circle cx="222" cy="90" r="2.2"/><circle cx="226" cy="198" r="2.2"/><circle cx="230" cy="50" r="2.2"/><circle cx="234" cy="158" r="2.2"/>
    <circle cx="238" cy="266" r="2.2"/><circle cx="242" cy="118" r="2.2"/><circle cx="246" cy="226" r="2.2"/><circle cx="250" cy="78" r="2.2"/>
    <circle cx="254" cy="186" r="2.2"/><circle cx="258" cy="38" r="2.2"/><circle cx="262" cy="146" r="2.2"/><circle cx="266" cy="254" r="2.2"/>
    <circle cx="270" cy="106" r="2.2"/><circle cx="274" cy="214" r="2.2"/><circle cx="278" cy="66" r="2.2"/><circle cx="282" cy="174" r="2.2"/>
  </g>
  <text x="158" y="306" text-anchor="middle" font-size="12" fill="currentColor">u = x_k / 64</text>
  <text x="14" y="158" text-anchor="middle" font-size="12" fill="currentColor" transform="rotate(-90 14 158)">v = x_(k+1) / 64</text>
</svg>
</Figure>

図の例では $\boldsymbol{c} = (-3, 7)$、$\ell = 10 \le (2\cdot64)^{1/2} \approx 11.3$ で、定理の上界 10 本に対して実際に点が載るのは 9 本です。

<Example id="ex-randu" title="RANDU の 3 つ組は 15 枚の平面に載る">
RANDU は $m = 2^{31}$、$a = 65539 = 2^{16} + 3$、$c = 0$ の乗算型合同法です。$2^{32} \equiv 0 \pmod{2^{31}}$ なので
$$
a^2 = 2^{32} + 6\cdot2^{16} + 9 \equiv 6(a - 3) + 9 = 6a - 9 \pmod{2^{31}},
$$
すなわち $9 - 6a + a^2 \equiv 0 \pmod m$ です。$\boldsymbol{c} = (9, -6, 1)$ は <Ref to="thm-marsaglia-planes" /> の合同式を満たし、$c = 0$ なので $\theta = 0$ です。よって連続する 3 つ組 $(u, v, w) = (x_k, x_{k+1}, x_{k+2})/m$ はすべて
$$
9u - 6v + w \in \mathbb{Z}
$$
を満たします。$u, v, w \in [0, 1)$ では $-6 < 9u - 6v + w < 10$ なので、取りうる整数は $-5, \ldots, 9$ の $15$ 個です。3 次元の点は 15 枚の平面の上にしかありません。隣り合う平面の間隔は $1/\sqrt{9^2 + 6^2 + 1^2} = 1/\sqrt{118} \approx 0.092$ で、立方体の中に厚さ約 $0.09$ の、点が 1 つもない層が並びます。定理の上界 $(3!\cdot2^{31})^{1/3} \approx 2344$ と比べると、RANDU の乗数は可能な範囲の中で極端に悪い選び方だったことが分かります。
</Example>

```python
# Python 3.9.6（標準ライブラリのみ）
a, m = 65539, 2**31
print("a^2 - 6a + 9 mod m =", (a * a - 6 * a + 9) % m)

x, planes, bad = 1, set(), 0          # 種 x_0 = 1
x1 = a * x % m
for _ in range(10**6):
    x2 = a * x1 % m
    s = 9 * x - 6 * x1 + x2           # 整数として厳密
    if s % m: bad += 1
    planes.add(s // m)
    x, x1 = x1, x2
print("線形関係が破れた回数:", bad)
print("現れた平面 9u - 6v + w = j の j:", sorted(planes))
print("枚数:", len(planes), "  Marsaglia の上界 (3! m)^(1/3) =", round((6 * m) ** (1 / 3), 1))
```

```text
a^2 - 6a + 9 mod m = 0
線形関係が破れた回数: 0
現れた平面 9u - 6v + w = j の j: [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
枚数: 15   Marsaglia の上界 (3! m)^(1/3) = 2344.4
```

$10^6$ 個の 3 つ組のすべてで線形関係が整数として成り立ち、15 枚の平面がすべて現れています。

<Remark id="rem-spectral-test" title="スペクトル検定">
乗数 $a$ の良し悪しは、各次元 $d$ について、点を覆う平行な超平面の族のうち間隔が最大のものの間隔 $1/\nu_d$（$\nu_d$ は <Ref to="thm-marsaglia-planes" /> の $\Lambda$ の最短ベクトルの Euclid ノルム）で測ります。これを**スペクトル検定**といい、線形合同法の乗数はこの値がよいものを計算機で探して選びます。各次元で達成できる最良の定数（Hermite 定数）や乗数の表は Knuth の本の 3.3.4 節と L'Ecuyer の論文 (1999) にあります。この章では最良の定数を定理としては扱いません。
</Remark>

## 5. $\mathbb{F}_2$ 上の線形生成器

ビット列を、2 元体 $\mathbb{F}_2 = \{0, 1\}$ 上のベクトルとみます。加法は排他的論理和（XOR）です。

<Definition id="def-f2-linear" title="F_2 上の線形生成器">
$n \ge 1$ とし、$A$ を $\mathbb{F}_2$ の元を成分とする $n$ 次正方行列とする。状態集合 $S = \mathbb{F}_2^n$、遷移 $T(x) = Ax$ の生成器を **$\mathbb{F}_2$ 上の線形生成器**という。次の 2 つが代表的である。

- **線形帰還シフトレジスタ**（LFSR）: ビットの列を漸化式 $b_{k+n} = \gamma_{n-1}b_{k+n-1} + \cdots + \gamma_0b_k$（$\gamma_i \in \mathbb{F}_2$）で作る。状態は $(b_k, \ldots, b_{k+n-1})$。
- **xorshift**: 状態を $n$ ビットの整数とみて、「左に $s$ ビットずらしたものとの XOR」「右に $s$ ビットずらしたものとの XOR」を何回か合成する。左シフトの行列を $L$、右シフトの行列を $R$ とすると、これらは $I + L^s$、$I + R^s$ という行列の積である。
</Definition>

零ベクトルは $A\boldsymbol{0} = \boldsymbol{0}$ で不動点なので、零でない種からの周期は高々 $2^n - 1$ です。この最大値が達成される条件は、$A$ の特性多項式 $\chi_A(t) = \det(tI - A) \in \mathbb{F}_2[t]$ で決まります。

<Definition id="def-primitive-f2" title="F_2 上の原始多項式">
次数 $n \ge 1$ のモニック多項式 $f \in \mathbb{F}_2[t]$ が**原始多項式**であるとは、$f$ が $\mathbb{F}_2[t]$ の既約元（<Ref to="mathematics/algebra/polynomial-rings-field-extensions#def-irreducible-prime" />）であり、かつ $f \mid t^k - 1$ となる最小の整数 $k \ge 1$ が $2^n - 1$ であることをいう。
</Definition>

<Aside type="caution">
同じ「原始多項式」という語が、係数の最大公約数が 1 の多項式（<Ref to="mathematics/algebra/polynomial-rings-field-extensions#def-primitive" />）にも使われます。この章の原始多項式は有限体の理論の用語で、それとは別の概念です。
</Aside>

次の定理のねらいは、「周期が最大になるか」という数列の問題を、「特性多項式が原始多項式か」という多項式の問題に言い換えることです。証明では、$A$ の多項式の全体を一つの環として扱い、周期が最大なら零でない元がすべて $A$ の冪で書けて逆元をもつ、という点から既約性を引き出します。

<Theorem id="thm-f2-linear-period" title="F_2 線形生成器の最大周期">
$A$ を $\mathbb{F}_2$ 上の $n$ 次正方行列とする。次の 2 条件は同値である。

1. 零でないある種 $x_0 \in \mathbb{F}_2^n$ について、状態の列 $x_k = A^kx_0$ が純周期的で周期が $2^n - 1$ である。
2. 特性多項式 $\chi_A$ が $\mathbb{F}_2$ 上の原始多項式である。

このとき、零でないすべての種について列は純周期的で周期は $2^n - 1$ である。
</Theorem>

<Proof of="thm-f2-linear-period">
$P = 2^n - 1$ と置き、$R = \{g(A) : g \in \mathbb{F}_2[t]\}$ とします。$R$ は $A$ の多項式の全体で、行列の和と積で閉じた可換な環です。Cayley–Hamilton の定理（<Ref to="mathematics/algebra/modules-over-pid#cor-cayley-hamilton" />）から $\chi_A(A) = 0$ です。

**2 ⇒ 1（およびすべての種についての主張）。** $f = \chi_A$ が原始多項式だとします。$f \mid t^P - 1$ なので $t^P - 1 = fq$ と書けて、$A^P - I = f(A)q(A) = 0$ です。よって $A$ は可逆で、$T$ は全単射、すべての列は純周期的です（<Ref to="prop-eventually-periodic" /> の 3）。$x_0 \ne 0$ と $1 \le k < P$ を取り、$A^kx_0 \ne x_0$ を示します。$g = t^k - 1$ と置くと、原始多項式の定義から $f \nmid g$ です。$\mathbb{F}_2[t]$ は単項イデアル整域なので（<Ref to="mathematics/algebra/ideals-and-quotient-rings#prop-pid-z-kx" /> の 2）、$f$ と $g$ が生成するイデアルは $(h)$ と書けます。$h \mid f$ で $f$ は既約なので、$h$ は単元か、$f$ の単元倍です。後者なら $f \mid g$ となるので、$h$ は単元、すなわち $uf + vg = 1$ となる多項式 $u, v$ があります。$A$ を代入すると $v(A)g(A) = I$ なので、$g(A) = A^k - I$ は可逆です。$x_0 \ne 0$ から $(A^k - I)x_0 \ne 0$、すなわち $A^kx_0 \ne x_0$ です。$A^Px_0 = x_0$ と合わせて、周期はちょうど $P$ です。

**1 ⇒ 2。** $x_0 \ne 0$ から始まる列が純周期的で周期 $P$ だとします。$x_0, \ldots, x_{P-1}$ は互いに異なります。また、ある $j$ で $x_j = 0$ なら以後ずっと $0$ で、純周期性から $x_0 = 0$ となってしまうので、どの $x_k$ も零ではありません。零でないベクトルは $P$ 個なので、
$$
\{x_0, x_1, \ldots, x_{P-1}\} = \mathbb{F}_2^n \setminus \{0\}
$$
です。$\mathbb{F}_2$ 線形写像 $\phi : R \to \mathbb{F}_2^n$ を $\phi(B) = Bx_0$ で定めます。

*$\phi$ は全単射。* 像は $\phi(A^k) = x_k$ と $\phi(0) = 0$ を含むので全射です。$Bx_0 = 0$ なら、$R$ が可換であることから $Bx_k = BA^kx_0 = A^kBx_0 = 0$ がすべての $k$ で成り立ち、$B$ はすべてのベクトルを $0$ に写すので $B = 0$ です。よって単射です。

*$R$ の零でない元はすべて $A$ の冪で、可逆。* $B \in R$、$B \ne 0$ なら $\phi(B) \ne 0$ なので $\phi(B) = x_k = \phi(A^k)$ となる $k$ があり、単射性から $B = A^k$ です。また $\phi(A^P) = x_P = x_0 = \phi(I)$ から $A^P = I$ なので、$A^k$ は逆元 $A^{P-k}$（$0 \le k < P$ に取る）をもちます。特に、$R$ の零でない 2 元の積は零ではありません。

*$g(A) = 0 \iff \chi_A \mid g$。* $\chi_A(A) = 0$ から $A^n$ は $I, A, \ldots, A^{n-1}$ の一次結合で、帰納的にすべての冪がそうです。よって $R$ は $I, \ldots, A^{n-1}$ で張られます。$\phi$ が同型なので $\dim_{\mathbb{F}_2}R = n$ で、この $n$ 個は一次独立です。つまり、次数が $n$ 未満の零でない多項式 $g$ は $g(A) \ne 0$ です。$J = \{g : g(A) = 0\}$ は $\mathbb{F}_2[t]$ のイデアルで、$\chi_A$ は $J$ の零でない元のうち次数が最小のものの 1 つなので、$J = (\chi_A)$ です（<Ref to="mathematics/algebra/ideals-and-quotient-rings#prop-pid-z-kx" /> の 2）。

*$\chi_A$ は既約。* $\chi_A$ は次数 $n \ge 1$ なので、零でも単元でもありません。$\chi_A = gh$ で $g, h$ がどちらも単元でない（次数が $1$ 以上）とすると、$g(A)h(A) = 0$ から $g(A) = 0$ または $h(A) = 0$ ですが、$g, h$ の次数は $n$ 未満なので矛盾です。

*位数。* $\chi_A \mid t^k - 1 \iff A^k = I \iff A^kx_0 = x_0$（$\phi$ の単射性）です。$x_0, \ldots, x_{P-1}$ が互いに異なり $x_P = x_0$ なので、これを満たす最小の $k \ge 1$ は $P$ です。以上で $\chi_A$ は原始多項式です。
</Proof>

最大周期の生成器がどの $n$ にも存在することは、有限体の乗法群が巡回群であることから従います。

<Proposition id="prop-primitive-exists" title="原始多項式の存在">
すべての $n \ge 1$ について、$\mathbb{F}_2$ 上の次数 $n$ の原始多項式が存在する。したがって、周期 $2^n - 1$ の $\mathbb{F}_2$ 線形生成器が存在する。
</Proposition>

<Proof of="prop-primitive-exists">
元の個数が $2^n$ の体 $K$ を取ります（<Ref to="mathematics/algebra/polynomial-rings-field-extensions#thm-finite-fields" /> の 2）。標数は $2$ なので $K$ は $\mathbb{F}_2 = \{0, 1\}$ を含み、$\mathbb{F}_2$ 上のベクトル空間です。その次元を $n'$ とすると元の個数は $2^{n'}$ なので、$n' = n$ です。乗法群 $K^\times$ は位数 $2^n - 1$ の巡回群なので（<Ref to="mathematics/algebra/polynomial-rings-field-extensions#thm-unit-cyclic" />）、生成元 $\alpha$ を取れます。

$K$ の基底を 1 つ固定し、$\mathbb{F}_2$ 線形写像 $x \mapsto \alpha x$ の表現行列を $A$ とします。種を $x_0 = 1$（に対応する座標ベクトル）とすると $x_k = \alpha^k$ で、$\alpha$ の位数が $2^n - 1$ なので、$x_0, \ldots, x_{2^n-2}$ は互いに異なり $x_{2^n-1} = x_0$ です。すなわち列は純周期的で周期は $2^n - 1$ です。<Ref to="thm-f2-linear-period" /> の 1 ⇒ 2 から、$\chi_A$ は次数 $n$ の原始多項式です。
</Proof>

<Example id="ex-lfsr-4" title="4 ビットの LFSR">
漸化式 $b_{k+4} = b_{k+1} + b_k$ の状態 $(b_k, b_{k+1}, b_{k+2}, b_{k+3})$ の遷移行列は
$$
A = \begin{pmatrix} 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ 1&1&0&0 \end{pmatrix}, \qquad \chi_A(t) = \det(tI - A) = t^4 + t + 1
$$
です（第 1 列についての余因子展開。$\mathbb{F}_2$ では $-1 = 1$）。種 $(1,0,0,0)$ からビットの列は
$$
1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1,\ 1, 0, 0, 0, \ldots
$$
と周期 $15 = 2^4 - 1$ で、$t^4 + t + 1$ は原始多項式です。一方、漸化式 $b_{k+4} = b_{k+3} + b_{k+2} + b_{k+1} + b_k$（特性多項式 $t^4 + t^3 + t^2 + t + 1$）では $1, 0, 0, 0, 1,\ 1, 0, 0, 0, 1, \ldots$ と周期 $5$ です。この多項式は $\mathbb{F}_2$ 上既約ですが、$(t - 1)(t^4 + t^3 + t^2 + t + 1) = t^5 - 1$ を割り切るので原始多項式ではありません。既約性だけでは最大周期になりません。
</Example>

次のコードは、16 ビットの xorshift について、零でない $65535$ 個の状態を軌道に分けて長さを調べます。

```python
# Python 3.9.6（標準ライブラリのみ）
def step(x, a, b, c):                 # 16 ビットの xorshift
    x ^= (x << a) & 0xFFFF
    x ^= x >> b
    x ^= (x << c) & 0xFFFF
    return x

def orbit_lengths(abc):               # 零でない全状態を軌道に分け、長さの集合を返す
    seen, lengths = bytearray(1 << 16), set()
    for x0 in range(1, 1 << 16):
        if seen[x0]:
            continue
        x, k = x0, 0
        while not seen[x]:
            seen[x] = 1
            x = step(x, *abc)
            k += 1
        lengths.add(k)
    return sorted(lengths)

for abc in ((7, 9, 8), (4, 3, 7), (1, 1, 1)):
    print(abc, "軌道の長さ:", orbit_lengths(abc))
```

```text
(7, 9, 8) 軌道の長さ: [65535]
(4, 3, 7) 軌道の長さ: [65535]
(1, 1, 1) 軌道の長さ: [1, 2, 4, 8, 16]
```

ずらし幅 $(7, 9, 8)$ と $(4, 3, 7)$ では零でない全状態が 1 つの軌道をなし、<Ref to="thm-f2-linear-period" /> からその遷移行列の特性多項式は原始多項式です。$(1, 1, 1)$ では軌道が短く分かれ、乱数には使えません。

<Remark id="rem-mersenne-twister" title="Mersenne Twister と線形性の限界">
Mersenne Twister（MT19937）は状態が $19937$ ビットの $\mathbb{F}_2$ 線形生成器で、周期は $2^{19937} - 1$、32 ビット精度で 623 次元均等分布です。$2^{19937} - 1$ は素数なので、原始性の判定が既約性の判定に帰着します。これらの証明は計算機による多項式の計算を含むので、ここでは定理として扱いません（Matsumoto–Nishimura の論文を参照してください）。一方、$\mathbb{F}_2$ 線形生成器の出力ビットは状態の線形関数なので、ビット列の線形複雑度や 2 元行列の階数を調べる検定には構造的に合格しません（L'Ecuyer–Simard の TestU01 の報告）。周期の長さと均等分布は、欠陥がないことを意味しません。
</Remark>

## 6. 飛び越しと並列ストリーム

### 6.1. 出力の攪拌と飛び越し

<Definition id="def-output-permutation" title="出力関数による攪拌">
遷移 $T$ は単純なまま、出力関数 $g$ で状態のビットを混ぜる設計を、出力の**攪拌**という。**PCG**（permuted congruential generator）は、法 $2^{64}$ または $2^{128}$ の全周期の線形合同法を遷移とし、状態の上位ビットで決まる量だけ下位側をずらす・回転させる関数（たとえば「上位と下位の XOR を取り、最上位の数ビットが指す量だけ回転する」）を出力関数とする生成器の族である。**カウンタ型**の生成器（Philox、Threefry）は、遷移を $T(x) = x + 1$ とし、出力関数にブロック暗号に似た全単射を使う。
</Definition>

PCG では、状態の列は線形合同法そのものなので、周期は <Ref to="thm-hull-dobell" /> で決まります。出力は上位ビットを使って下位ビットを選び直すので、<Ref to="prop-low-bits" /> の下位ビットの短い周期は出力に現れません。ただし、出力の統計的な性質がよいことは定理ではなく、検定による経験的な事実です。

<Definition id="def-jump-stream" title="飛び越しとストリーム">
生成器 $(S, T, g)$ について、状態 $s$ と整数 $k \ge 0$ から $T^k(s)$ を、$T$ を $k$ 回施すより少ない手間で求めることを**飛び越し**（jump-ahead）という。開始状態 $s$ と長さ $L$ に対して、状態の並び $s, T(s), \ldots, T^{L-1}(s)$（およびそこからの出力）を、$s$ から始まる長さ $L$ の**ストリーム**という。
</Definition>

<Proposition id="prop-jump-ahead" title="飛び越しの計算量">
1. 線形合同法について、$T^k(x) = (a_kx + c_k) \bmod m$ となる組 $(a_k, c_k) = (a^k \bmod m,\ cS_k(a) \bmod m)$ は、$m$ を法とする乗算と加算を合わせて高々 $6\lfloor\log_2 k\rfloor$ 回（$k \ge 1$）で求まる。
2. $\mathbb{F}_2$ 上の $n$ 次元の線形生成器について、$A^k$ は高々 $2\lfloor\log_2 k\rfloor$ 回の行列の積で求まる。
</Proposition>

<Proof of="prop-jump-ahead">
アフィン写像 $x \mapsto \alpha x + \gamma$ を組 $(\alpha, \gamma)$ で表します。$(\alpha, \gamma)$ の後に $(\alpha', \gamma')$ を施すと $\alpha'(\alpha x + \gamma) + \gamma' = \alpha'\alpha x + (\alpha'\gamma + \gamma')$ なので、合成は
$$
(\alpha', \gamma') \circ (\alpha, \gamma) = (\alpha'\alpha,\ \alpha'\gamma + \gamma') \pmod m
$$
で、乗算 2 回と加算 1 回の計 3 回の演算です。写像の合成は結合的なので、$T = (a, c)$ の $k$ 回の合成 $T^k$ について $T^{2q} = T^q \circ T^q$、$T^{2q+1} = T \circ (T^q \circ T^q)$ が成り立ちます。$T^k$ を求めるのに要る合成の回数を $M(k)$ とすると、$M(1) = 0$、$M(k) \le M(\lfloor k/2\rfloor) + 2$ なので、$k$ についての帰納法で $M(k) \le 2\lfloor\log_2 k\rfloor$ です（$\lfloor\log_2\lfloor k/2\rfloor\rfloor = \lfloor\log_2 k\rfloor - 1$ を使います）。演算の回数は $3M(k) \le 6\lfloor\log_2 k\rfloor$ です。$T^k = (a^k, cS_k(a))$ であることは <Ref to="lem-closed-form" /> です。

2 も同じ漸化式 $A^{2q} = A^qA^q$、$A^{2q+1} = A\cdot A^qA^q$ によります。
</Proof>

$\mathbb{F}_2$ 線形生成器で $n$ が大きいとき（MT19937 では $n = 19937$）、行列の積は重すぎます。その場合は多項式 $t^k \bmod \chi_A(t) = \sum_{i<n} r_it^i$ を繰り返し二乗で求め、$A^kx = \sum_i r_iA^ix$ を $T$ の $n$ 回の適用で計算します（$\chi_A(A) = 0$ なので $A^k = \sum_i r_iA^i$ です）。

### 6.2. ストリームの重なり

並列計算では、各ワーカーが別々のストリームを使います。同じ状態を 2 つのワーカーが使うと、その部分の標本は完全に相関し、独立な標本としての誤差評価が崩れます。ストリームの開始点の決め方は 2 通りあります。

<Figure caption="状態を円周上に並べたサイクルと、長さ L の区間（太い弧）で表した s = 3 本のストリーム。左は開始点が 120 度ずつ離れていて弧が重ならない。右は 1 本目と 2 本目の開始点が 60 度しか離れておらず、弧が重なる。">
<svg viewBox="0 0 460 190" width="100%" role="img" aria-label="円周上の 3 本の弧が、左では重ならず右では 2 本が重なる様子" style="max-width: 460px">
  <g fill="none" stroke="currentColor" stroke-width="1" opacity="0.5">
    <circle cx="110" cy="90" r="50"/><circle cx="350" cy="90" r="50"/>
  </g>
  <g fill="none" stroke="var(--sl-color-accent)" stroke-width="9" opacity="0.5" stroke-dasharray="87.27 400">
    <circle cx="110" cy="90" r="50" transform="rotate(0 110 90)"/>
    <circle cx="110" cy="90" r="50" transform="rotate(120 110 90)"/>
    <circle cx="110" cy="90" r="50" transform="rotate(240 110 90)"/>
    <circle cx="350" cy="90" r="50" transform="rotate(0 350 90)"/>
    <circle cx="350" cy="90" r="50" transform="rotate(60 350 90)"/>
    <circle cx="350" cy="90" r="50" transform="rotate(200 350 90)"/>
  </g>
  <text x="110" y="173" text-anchor="middle" font-size="12" fill="currentColor">重ならない</text>
  <text x="350" y="173" text-anchor="middle" font-size="12" fill="currentColor">重なる（濃い部分）</text>
</svg>
</Figure>

<Theorem id="thm-stream-overlap" title="ストリームの重なりの確率">
遷移 $T$ が状態を長さ $P$ の 1 つのサイクル $C$ として巡回させる生成器を考える。$s \ge 2$ 本のストリームの長さを $L$ とし、$2L - 1 \le P$ とする。

1. 開始状態 $\sigma_1, \ldots, \sigma_s$ が $C$ 上の一様分布に従う互いに独立な確率変数ならば、ある 2 本のストリームが同じ状態を含む確率 $p_{\mathrm{overlap}}$ は $p_{\mathrm{overlap}} \le \dfrac{s(s-1)(2L-1)}{2P} \le \dfrac{s(s-1)L}{P}$ を満たす。
2. 開始状態を $\sigma_j = T^{jJ}(\sigma_0)$（$j = 0, \ldots, s-1$）と取り、$J \ge L$ かつ $sJ \le P$ ならば、どの 2 本のストリームも同じ状態を含まない。
</Theorem>

<Proof of="thm-stream-overlap">
$C$ の 1 点 $s^*$ を固定し、$C$ の状態 $T^i(s^*)$ を $i \in \mathbb{Z}/P\mathbb{Z}$ と同一視します。$T$ は $i \mapsto i + 1$ になり、開始位置 $\sigma$ のストリームは $\{\sigma, \sigma+1, \ldots, \sigma + L - 1\}$ です。開始位置 $\sigma, \sigma'$ の 2 本が共通の状態をもつのは、$\sigma + i = \sigma' + i'$ となる $0 \le i, i' < L$ があるとき、すなわち
$$
\sigma' - \sigma \in D := \{-(L-1), \ldots, -1, 0, 1, \ldots, L-1\} \pmod P
$$
のときです。$D$ の元は $2L - 1 \le P$ 個以下です。

**1.** $i \ne j$ とします。$\sigma_i, \sigma_j$ は独立で一様なので（独立性は <Ref to="mathematics/probability/random-variables#def-independence" />）、各 $\delta \in \mathbb{Z}/P\mathbb{Z}$ について
$$
\mathbb{P}(\sigma_j - \sigma_i = \delta) = \sum_{\tau} \mathbb{P}(\sigma_i = \tau)\,\mathbb{P}(\sigma_j = \tau + \delta) = P\cdot\frac1P\cdot\frac1P = \frac1P
$$
です。よって 2 本が重なる確率は $|D|/P \le (2L-1)/P$ です。組 $\{i, j\}$ は $s(s-1)/2$ 通りあるので、確率の劣加法性から、どれかの組が重なる確率は $\frac{s(s-1)}{2}\cdot\frac{2L-1}{P}$ 以下です。$2L - 1 < 2L$ から第 2 の不等式が従います。

**2.** $i < j$ について $\sigma_j - \sigma_i = (j - i)J$ で、$J \le (j-i)J \le (s-1)J \le P - J$ です。$L \le J$ なので、整数として $L \le (j-i)J \le P - L$ です。$D$ の元を $0$ 以上 $P$ 未満の代表で表すと $\{0, \ldots, L-1\} \cup \{P - L + 1, \ldots, P - 1\}$ なので、$(j-i)J \bmod P \notin D$ です。
</Proof>

1 の論法は、ハッシュ値の衝突についての誕生日境界（<Ref to="computer-science/software-engineering/version-control-git#prop-birthday" />）と同じ「組の数 × 1 組あたりの確率」です。違いは、点の一致ではなく長さ $L$ の区間の重なりを数えるので、因子 $2L - 1$ が付くことです。$P = 2^{128}$、$s = 2^{20}$ 本、各 $L = 2^{40}$ 個のとき、上界は
$$
\frac{s(s-1)L}{P} < \frac{2^{20}\cdot2^{20}\cdot2^{40}}{2^{128}} = 2^{-48} \approx 3.6\times10^{-15}
$$
です。一方 $P = 2^{64}$ で同じ規模の計算をすると上界は $1$ を超え、何も保証しません。並列計算に周期 $2^{64}$ の生成器を使えない理由です。

<Definition id="def-seed-sequence" title="種の展開">
利用者が与える短い種（整数や整数の列）から、生成器の状態全体をハッシュに似た混合関数で作ることを**種の展開**という。NumPy の `SeedSequence` は、種と「何番目の子か」を表す番号の列を混ぜて状態を作り、`spawn` は番号の違う子を返す。
</Definition>

種の展開で作った開始状態は、サイクル上の独立で一様な点のように振る舞うと期待され、そのモデルのもとで <Ref to="thm-stream-overlap" /> の 1 が重なりの確率を評価します（モデルが成り立つこと自体は定理ではありません）。これに対して飛び越しは、2 のとおり重なりが起きないことを保証します。

```python
# Python 3.9.6, NumPy 2.0.2
import numpy as np
from fractions import Fraction
from math import log2

ss = np.random.SeedSequence(20261008)
children = ss.spawn(3)                               # 種の展開によるストリーム
print([round(np.random.default_rng(s).random(), 6) for s in children])

bg = np.random.PCG64(20261008)
jumped = [bg.jumped(j) for j in range(3)]            # 飛び越しによるストリーム
print([round(np.random.Generator(b).random(), 6) for b in jumped])

P, J = 2**128, 210306068529402873165736369884012333109   # PCG64 の周期と jumped の歩数
starts = sorted(j * J % P for j in range(2**10))
gap = min((b - a) % P for a, b in zip(starts, starts[1:] + starts[:1]))
print(f"jumped 1024 本の開始点の最小間隔 = 2^{log2(gap):.2f}")

def bound(logP, logs, logL):                         # s(s-1)L/P を有理数で計算
    s, L, P = 2**logs, 2**logL, 2**logP
    return float(Fraction(s * (s - 1) * L, P))
for logP, logs, logL in ((64, 10, 30), (64, 20, 40), (128, 20, 40), (128, 30, 64)):
    b = bound(logP, logs, logL)
    print(f"P=2^{logP:<4d} s=2^{logs:<3d} L=2^{logL:<3d} 上界 = {b:.3e}" + ("  (1 を超え、情報がない)" if b >= 1 else ""))
```

```text
[0.953464, 0.617841, 0.494969]
[0.125636, 0.622824, 0.591673]
jumped 1024 本の開始点の最小間隔 = 2^116.89
P=2^64   s=2^10  L=2^30  上界 = 6.098e-05
P=2^64   s=2^20  L=2^40  上界 = 6.554e+04  (1 を超え、情報がない)
P=2^128  s=2^20  L=2^40  上界 = 3.553e-15
P=2^128  s=2^30  L=2^64  上界 = 6.250e-02
```

NumPy の `PCG64.jumped(j)` は、状態を $jJ$ 歩（$J \approx 2^{128}/\varphi$、$\varphi$ は黄金比。値は NumPy の文書に記載）進めます。等間隔ではないので <Ref to="thm-stream-overlap" /> の 2 の仮定とは形が違いますが、証明の本質は開始点どうしのサイクル上の距離が $L$ 以上であることなので、1024 本なら各ストリームが $2^{116}$ 個以下の乱数を使うかぎり重なりません（出力の 3 行目）。

<Remark id="rem-empirical" title="統計的な検定と、定理にならない性質">
この章の定理が保証するのは周期・格子構造・重なりの確率だけです。「連続する出力が独立な一様乱数と区別できない」ことは、TestU01 の BigCrush のような検定群に合格するという形でしか確かめられず、合格は欠陥がないことの証明ではありません。PCG・Philox・xoshiro などの現在の標準的な生成器の統計的な性質は、それぞれの設計者の報告と独立な検定の結果によります。Monte Carlo 計算では、別の族の生成器で計算をやり直して結果が標準誤差の範囲で一致することを確かめるのが、実務上の検算です。
</Remark>

<Remark id="rem-cryptographic" title="暗号論的擬似乱数との違い">
この章の生成器は、出力の一部から状態を復元できます。線形合同法は連続する数個の出力から $a, c$ と状態が求まり、$\mathbb{F}_2$ 線形生成器は $n$ ビット程度の出力から連立一次方程式で状態が求まります。暗号に使う擬似乱数には、出力をいくら観測しても次の出力を多項式時間で予測できないこと（予測不能性）が要求され、これは計算量的な仮定に基づく別の理論です。鍵やトークンの生成には、オペレーティングシステムが提供する暗号論的な乱数源（Python では `secrets` モジュール）を使ってください。
</Remark>

## 7. 演習

<Exercise id="exr-hull-dobell-apply" difficulty="易">
$m = 2^3\cdot 3^2\cdot 5 = 360$ の線形合同法が全周期をもつ乗数 $a$（$1 \le a \le 360$）をすべて求め、その個数を答えよ。また、全周期をもつ組 $(a, c)$（$0 \le c < 360$）の個数を求めよ。
<Solution>
<Ref to="thm-hull-dobell" /> の条件 2・3 は $a \equiv 1 \pmod 3$、$a \equiv 1 \pmod 5$、$a \equiv 1 \pmod 4$（これは $a \equiv 1 \pmod 2$ を含みます）です。$3, 4, 5$ は対ごとに互いに素なので、これは $a \equiv 1 \pmod{60}$ と同値です（<Ref to="mathematics/number-theory/elementary-number-theory#cor-crt-multi" />）。$1 \le a \le 360$ では $a = 1, 61, 121, 181, 241, 301$ の $6$ 個です。条件 1 は $\gcd(c, 360) = 1$ で、そのような $c$ は $\varphi(360) = 360\cdot\frac12\cdot\frac23\cdot\frac45 = 96$ 個です。組の個数は $6 \times 96 = 576$ で、§3 のコードの出力と一致します。
</Solution>
</Exercise>

<Exercise id="exr-lcg-pairs" difficulty="標準">
$e \ge 4$ を偶数とし、$m = 2^e$、$a = 2^{e/2} + 1$、$c$ を奇数とする。(1) この線形合同法が全周期をもつことを示せ。(2) 連続する 3 つ組 $(x_k, x_{k+1}, x_{k+2})/m$ が、$c$ によらず高々 $4$ 枚の平行な平面に載ることを示せ。
<Solution>
(1) $a - 1 = 2^{e/2}$ は $e \ge 4$ から $4$ の倍数で、$c$ は奇数なので、<Ref to="thm-hull-dobell" /> の 3 条件が成り立ちます。

(2) $a^2 = 2^e + 2\cdot2^{e/2} + 1 \equiv 2(a - 1) + 1 = 2a - 1 \pmod{2^e}$ なので、$1 - 2a + a^2 \equiv 0 \pmod m$ です。$\boldsymbol{c} = (1, -2, 1)$ は <Ref to="thm-marsaglia-planes" /> の合同式を満たし、$\ell = 4$ です。同じ定理の後半から、3 つ組はすべて $u - 2v + w \in \theta + \mathbb{Z}$ を満たし、載りうる平面は高々 $4$ 枚です。乗数 $a$ が $2$ の冪に近いと、低い次数の関係式 $(a - 1)^2 \equiv 0$ が生じて 3 次元の格子がつぶれます。RANDU（<Ref to="ex-randu" />）の $(a - 3)^2 \equiv 0$ と同じ仕組みです。全周期であること（1 次元の性質）は、3 次元の並びについて何も保証しません。
</Solution>
</Exercise>

<Exercise id="exr-primitive-degree-3" difficulty="標準">
$\mathbb{F}_2$ 上の次数 $3$ のモニック多項式のうち原始多項式であるものをすべて求め、そのそれぞれについて LFSR の漸化式と、種 $(1, 0, 0)$ からの 1 周期分のビット列を書け。
<Solution>
次数 3 のモニック多項式は $t^3 + \alpha t^2 + \beta t + \gamma$ の $8$ 個です。既約であるためには根 $0, 1$ をもたないこと（3 次式が可約なら 1 次の因子をもつので、これで十分）、すなわち $\gamma = 1$ かつ $1 + \alpha + \beta + 1 \ne 0$、つまり $\alpha + \beta = 1$ が必要十分です。該当するのは $t^3 + t + 1$ と $t^3 + t^2 + 1$ の 2 つで、残りの 6 個は既約でないので原始多項式ではありません。

漸化式 $b_{k+3} = \gamma_2b_{k+2} + \gamma_1b_{k+1} + \gamma_0b_k$ の状態 $(b_k, b_{k+1}, b_{k+2})$ の遷移行列は、第 3 行が $(\gamma_0, \gamma_1, \gamma_2)$ の同伴行列で、特性多項式は $t^3 + \gamma_2t^2 + \gamma_1t + \gamma_0$ です（<Ref to="ex-lfsr-4" /> と同じ余因子展開）。2 つの多項式に対応する列を種 $(1, 0, 0)$ から計算します。

- $t^3 + t + 1$: $b_{k+3} = b_{k+1} + b_k$。$1, 0, 0, 1, 0, 1, 1,\ 1, 0, 0, \ldots$
- $t^3 + t^2 + 1$: $b_{k+3} = b_{k+2} + b_k$。$1, 0, 0, 1, 1, 1, 0,\ 1, 0, 0, \ldots$

どちらも 8 番目から最初の状態 $(1, 0, 0)$ に戻り、それまでの 7 個の状態は互いに異なるので、列は純周期的で周期は $7 = 2^3 - 1$ です。<Ref to="thm-f2-linear-period" /> の 1 ⇒ 2 から、2 つとも原始多項式です。7 ビットの中に $1$ が $4$ 個、$0$ が $3$ 個あり、零でない 3 ビットの状態 $7$ 通りがちょうど 1 回ずつ現れます。
</Solution>
</Exercise>

<Exercise id="exr-overlap-design" difficulty="易">
周期 $P = 2^{128}$ の生成器で、$s$ 本のストリームを種の展開で作り、各ストリームで $L = 2^{50}$ 個の乱数を使う。重なりの確率の上界 $s(s-1)L/P$ を $10^{-12}$ 以下にしたい。$s$ はどこまで増やせるか。
<Solution>
$s(s-1) < s^2$ なので $s^2 \cdot 2^{50}/2^{128} \le 10^{-12}$、すなわち $s^2 \le 10^{-12}\cdot2^{78}$ であれば十分です。$2^{78} \approx 3.02\times10^{23}$ なので $s^2 \le 3.02\times10^{11}$、$s \le 5.5\times10^5 \approx 2^{19}$ です。約 50 万本までなら、<Ref to="thm-stream-overlap" /> の 1 のモデルのもとで重なりの確率は $10^{-12}$ 以下です。これを超える規模では、飛び越しで分ける（<Ref to="thm-stream-overlap" /> の 2）か、状態のより大きい生成器を使います。
</Solution>
</Exercise>

## 参考文献

1. D. E. Knuth, *The Art of Computer Programming, Volume 2: Seminumerical Algorithms*, 3rd ed., Addison-Wesley, 1997 — 第 3 章（3.2 節 線形合同法の周期、3.3.4 節 スペクトル検定）。
2. T. E. Hull, A. R. Dobell, "Random number generators", *SIAM Review* 4 (1962), 230–254.
3. G. Marsaglia, "Random numbers fall mainly in the planes", *Proceedings of the National Academy of Sciences* 61 (1968), 25–28.
4. M. Matsumoto, T. Nishimura, "Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator", *ACM Transactions on Modeling and Computer Simulation* 8 (1998), 3–30.
5. P. L'Ecuyer, R. Simard, "TestU01: A C library for empirical testing of random number generators", *ACM Transactions on Mathematical Software* 33 (2007), Article 22.
6. P. L'Ecuyer, "Tables of linear congruential generators of different sizes and good lattice structure", *Mathematics of Computation* 68 (1999), 249–260.
7. G. Marsaglia, "Xorshift RNGs", *Journal of Statistical Software* 8 (2003), Issue 14.
8. J. K. Salmon, M. A. Moraes, R. O. Dror, D. E. Shaw, "Parallel random numbers: as easy as 1, 2, 3", *Proceedings of SC '11* (2011).
9. NumPy の文書 [Random sampling: Parallel random number generation](https://numpy.org/doc/stable/reference/random/parallel.html) — `SeedSequence.spawn` と `jumped` の使い方。