キャリー付き乗算

キャリー付き乗算 (キャリーつきじょうさん、英: multiply-with-carry, MWC) は、George Marsagliaにより開発された整数疑似乱数生成用の手法である。乱数種には2〜数千の値を必要とする。MWC の主な長所は、単純な整数演算からなっており非常に高速に動作するという点と、2⁶⁰ - 2^2000000 という非常に長周期であるという点である。

基本理論

MWC 列は b を法とする剰余からなる。b = 2³² とするのが普通だが、これはコンピュータで扱う整数が通常そのようになっているからである。ただ、b = 2³² − 1とすることもある。これは、2³² − 1 を法とする演算は 2³² を法とする演算を少し変更するだけで済み、さらに、b = 2³² の MWC の理論にはいくつか厄介な問題があるが、b = 2³² − 1 ではそれを回避できるからである。

その最も一般的な形では、lag-r MWC 生成器は基数 b、乗数 a、そして r + 1 個の乱数種を必要とする。乱数種は r 個の b の剰余

x₀, x₁, x₂ ,..., x_r−1

と最初のキャリー c_r−1 < a である。

そして、lag-r MWC 列は

$x_{n}=(ax_{n-r}+c_{n-1})\,{\bmod {\,}}b,\ c_{n}=\left\lfloor {\frac {ax_{n-r}+c_{n-1}}{b}}\right\rfloor ,\ n\geq r$

で定義される x_n, c_n のペアの数列であり、MWC 生成器の出力は以下の x の列となる。

x_r , x_r+1 , x_r+2, ...

lag-r MWC 生成器の周期は ab^r − 1 を法とする数の乗法群における b の位数である。通例、a は b の位数を大きくできるよう、p = ab^r − 1 が素数になるように選ぶ。b = 2³² は p = ab^r − 1 の原始根とはならないので、b = 2³² を基数とした MWC 生成器の周期は、MWC の最大周期 p = ab^r − 2 になることはない。これが、b = 2³² − 1 の方が有利になる点の1つである。

Couture と l'Ecuyer (1997) により、MWC 生成器には最上位ビットが少し偏っているという理論的な問題があると指摘された。ただし、この問題は相補的キャリー付き乗算により解決されている。「相補的 MWC を使えば、最上位ビットは均一に出ることが分かるだろう。つまり、全周期中で0と1とが同等の頻度で現れ、その傾向も MWC 生成器間に関連性はない。」彼らはビットの偏り具合についてそれ以上詳しくは語っていないようである。相補的 MWC 生成器は計算時間が若干増えるため、実装の要求によってどちらを使うか決めるといいだろう。

線形合同法との比較

線形合同法は通常以下で定義される。

$x_{n}=ax_{n-1}+c\ {\bmod {\,}}2^{32}$

ほとんどの演算プロセッサでは乗数 a と現在の x は32ビットレジスタに置け、積は64ビットの連結レジスタ上に求められる。積の下位32ビット、すなわち

$a\times x\ {\bmod {\,}}2^{32}$

は、右側のレジスタを参照するだけで得られる。同様に、これに32ビットの c を足してやれば自然に ax+c mod 2³² が求まる。もし a mod 8 が 1 か 5 で c が奇数なら、2³² を法とする合同数列の周期は 2³² となる。

一方 lag-1 MWC は、線形合同法とほぼ同じ演算を使うだけで約 2⁶² の周期を実現する。違うのは、MWC では64ビット積の上半分も利用している点である。これは 次の キャリー c として使われる。一方、線形合同法では c は固定値である。ax+c を64ビット値として求め、上半分を次の c として使い、下半分を x として使う訳である。

乗数 a が与えられれば、x, c のペアを入力として、新たなペアが作成される。

$x\leftarrow (ax+c)\,{\bmod {\,}}2^{32},\ \ c\leftarrow \left\lfloor {\frac {ax+c}{2^{32}}}\right\rfloor$

もし x と c が両方とも 0 でなければ、MWC 列の周期は ab − 1 を法とする剰余の乗法群における b = 2³² の位数、すなわち、b³²ⁿ = 1 mod (ab − 1) となる最小の n になる。もし 28〜31ビットの a を ab − 1 が安全素数となるよう選ぶと、ab − 1 と ab/2 − 1 は両方とも素数になり、周期は ab/2 − 1 となり、2⁶² に近づく。実際にはとりうる32ビット値のペア x, c の数となるだろう。

以下は、上記の安全素数の条件を満たす a の最大値の例である。

a のビット数	b	ab-1 が安全素数となる a の最大値	周期
15	216	31,743	1,040,154,623
16	216	64,545	2,115,010,559
31	232	2,147,483,085	4,611,684,809,394,094,079
32	232	4,294,966,893	*

次に、10で割った余りという単純な例を使って線形合同法と MWC の比較を行う。ここではレジスタは10進数1桁のサイズであるとし、連結レジスタは10進数2桁であるとする。

$x_{0}=1$ から始めると、合同数列

$x_{n}=7x_{n-1}+3\,{\bmod {\,}}10$

は連結レジスタ上で

$10,03,24,31,10,03,24,31,10,\ldots$

の数列を持ち、x の出力列（右のレジスタ）の周期は4となる。

$0,3,4,1,0,3,4,1,0,3,4,1,\ldots$

$x_{0}=1$, $c_{0}=3$ から始めると、MWC 列

$x_{n}=(7x_{n-1}+c_{n-1})\,{\bmod {\,}}10,\ c_{n}=\left\lfloor {\frac {7x_{n-1}+c_{n-1}}{10}}\right\rfloor ,$

は連結レジスタ上で

$31,10,01,07,49,67,55,40,04,28,58,61,13,22,16,43,25,37,52,19,64,34,\,31,10,01,07,...$

の数列を持ち、x の出力列（右のレジスタ）の周期は22となる。

$1,0,1,7,9,7,5,0,4,8,8,1,3,2,6,3,5,7,2,9,4,4,\,1,0,1,7,9,7,5,0,...$

x の繰り返し部分を逆順にすると

$449275\cdots 97101\,449275\cdots 9710144\cdots$

となり、これは a=7, b=10, j=31 とした時の j/(ab−1) の値

${\frac {31}{69}}=.4492753623188405797101\,4492753623\ldots$

の小数部と一致する点に注目したい。

これは一般的にもなりたち、lag-r MWC により生成された x の列

$x_{n}=(ax_{n-r}+c_{n-1}){\bmod {\,}}b\,,\ \ c_{n}=\left\lfloor {\frac {ax_{n-r}+c_{n-1}}{b}}\right\rfloor ,$

は、逆順にすると、ある 0 < j < ab^r において j/(ab^r − 1) の基数 b での小数部と一致する。

さらに、合同数列を

$x_{n}=7x_{n-1}\,{\bmod {\,}}69$

で定義した場合、$x_{0}=34$ から始めると以下の周期22の数列を得るが、

$31,10,1,7,49,67,55,40,4,28,58,61,13,22,16,43,25,37,52,19,64,34,\,31,10,1,7,\ldots$

これを10で割った余りを求めると

$1,0,1,7,9,7,5,0,4,8,8,1,3,2,6,3,5,7,2,9,4,4,\,1,0,1,7,9,7,5,0,\ldots$

となり、上記の MWC 列に一致する。

これは一般的にも成り立ち（ただし、これは lag-1 MWC でしか成り立たないが）、初期値が $x_{0}$, $c_{0}$ の lag-1 MWC 列

$x_{n}=(ax_{n-1}+c_{n-1})\,{\bmod {b}}\,,\ \ c_{n}=\left\lfloor {\frac {ax_{n-1}+c_{n-1}}{b}}\right\rfloor$

により得られる数列 $x_{1},x_{2},\ldots$ は、合同数列 y_n = ay_n − 1 mod(ab − 1) を b で割った余りに一致する。

初期値 y₀ の選び方は、単に x のサイクルをずらすだけにすぎない。

相補的なキャリー付き乗算

lag-r MWC の周期を決めるのは、通常 p=ab^r − 1 が素数となる乗数 a の選び方である。もし p が安全素数なら、b の位数は p − 1 か (p − 1)/2 となる。ただ、b mod p の位数を求めるには p − 1 を因数分解する必要があるだろうが、これを因数分解するのは難しいだろう。

しかし、p = ab^r + 1 の形の素数であれば、p − 1 = ab^r は簡単に因数分解できる。そのため、p = ab^r + 1 を素数とする場合の b によって規定される MWC を使えば、MWC の周期を求めるのに必要な計算数論は著しく簡単になるだろう。

幸いな事に、MWC を少し変更するだけで、ab^r + 1 の形の素数を作る事ができる。これを相補的なキャリー付き乗算(CMWC)と呼ぶ。

必要な入力は lag-r MWC と同じで、乗数 a と 基数 b、そして r + 1 個の乱数種

x₀, x₁, x₂, ..., x_r−1, c_r − 1

である。新しいペア x, c を生成する方法は、以下のように少し変更される。

$x_{n}=(b-1)-(ax_{n-r}+c_{n-1})\,{\bmod {\,}}b,\ c_{n}=\left\lfloor {\frac {ax_{n-r}+c_{n-1}}{b}}\right\rfloor$

つまり、新しい x を作る際に、補数 (b − 1) − x を使うのである。

CMWC 生成器により作られる数列 x の周期は ab^r + 1 を法とする剰余の乗法群における b の位数になり、x を逆順にしたものは、ある 0 < j < ab^r において j/(ab^r + 1) の基数 b での小数部と一致する。

lag-r CMWC を使うと、r が 512, 1024, 2048 と大きくなっても、簡単に周期を求める事ができる。（r を 2 の累乗にすると、r 個の最新の x を保持する配列内の要素へのアクセスが若干簡単に（そして速く）なる。）

以下に、いくつかの例を示す。

b=2³² の時、lag-1024 CMWC

$x_{n}=(b-1)-(ax_{n-1024}+c_{n-1})\,{\bmod {\,}}b,\ c_{n}=\left\lfloor {\frac {ax_{n-1024}+c_{n-1}}{b}}\right\rfloor$

の周期は a$\cdot$2³²⁷⁶⁵² となり、109111 or 108798 or 108517 の a に対して約 10⁹⁸⁶⁷ となる。

b=2³² で a = 3636507990 の時、p = ab¹³⁵⁹ − 1 は安全素数となるので、MWC 列の周期は 3636507990$\cdot$2⁴³⁴⁸⁷ $\approx$10¹³¹⁰¹ となる。

b = 2³² の時、CMWC 生成器で現在見つかっている中で最大に近い周期を持つものに、素数 p = 15455296b⁴²⁶⁵⁸ + 1 を元にしたものがあり、b の位数は 241489*2^1365056 すなわち約 10⁴¹⁰⁹²⁸ である。

関連項目

• 疑似乱数

参考文献

• G. Marsaglia and A. Zaman, A new class of random number generators, Annals of Applied Probability V. 1, No. 3, 462-480
• G. Marsaglia, Random number generators, Journal of Modern Applied Statistical Methods,V. 2, May 2003. http://tbf.coe.wayne.edu/jmasm/vol2_no1.pdf
• G. Marsaglia, On the randomness of Pi and other decimal expansions, Interstat, October 2005, #5, http://interstat.statjournals.net/YEAR/2005/articles/0510005.pdf
• Couture, Raymond; L'Ecuyer, Pierre (1997), “Distribution properties of Multiply-with-carry random number generators”, Mathematics of Computation 66 (218): 591–607, doi:10.1090/S0025-5718-97-00827-2 

外部リンク

• C++による実装（英語）
• Pythonによる実装（英語）