キャリー付き乗算

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

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

その最も一般的な形では、**lag-r** MWC 生成器は**基数** *b*、**乗数** *a*、そして *r&nbsp;+&nbsp;1* 個の**乱数種**を必要とする。乱数種は *r* 個の *b* の剰余
:*x*0, *x*1, *x*2 ,..., *x**r*&minus;1
と最初の**キャリー** *c**r*&minus;1  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\ge r
で定義される *x**n*,&nbsp;*c**n* のペアの数列であり、MWC 生成器の出力は以下の *x* の列となる。
:*x**r* , *x**r*+1 , *x**r*+2, ...
lag-*r* MWC 生成器の周期は *ab<sup>r</sup>*&nbsp;&minus;&nbsp;1 を法とする数の乗法群における *b* の位数である。通例、*a* は *b* の位数を大きくできるよう、*p = ab<sup>r</sup>*&nbsp;&minus;&nbsp;1 が素数になるように選ぶ。*b* = 232 は *p = ab<sup>r</sup>*&nbsp;&minus;&nbsp;1 の原始根とはならないので、*b* = 232 を基数とした MWC 生成器の周期は、MWC の最大周期 *p = ab<sup>r</sup>*&nbsp;&minus;&nbsp;2 になることはない。これが、*b* = 232&nbsp;&minus;&nbsp;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 232 が求まる。もし *a* mod 8 が 1 か 5 で *c* が奇数なら、232 を法とする合同数列の周期は 232 となる。

一方 lag-1 MWC は、線形合同法とほぼ同じ演算を使うだけで約 262 の周期を実現する。違うのは、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*&nbsp;&minus;&nbsp;1 を法とする剰余の乗法群における *b* = 232 の位数、すなわち、*b*32*n* = 1 mod (*ab*&nbsp;&minus;&nbsp;1) となる最小の *n* になる。もし 28〜31ビットの *a* を *ab*&nbsp;&minus;&nbsp;1 が安全素数となるよう選ぶと、*ab*&nbsp;&minus;&nbsp;1 と *ab*/2&nbsp;&minus;&nbsp;1 は両方とも素数になり、周期は *ab*/2&nbsp;&minus;&nbsp;1 となり、262 に近づく。実際にはとりうる32ビット値のペア *x*, *c* の数となるだろう。

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

次に、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\cdots97101\,449275\cdots9710144\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 *r* において *j*/(*ab**r*&nbsp;&minus;&nbsp;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*&nbsp;&minus;&nbsp;1 mod(*ab*&nbsp;&minus;&nbsp;1) を *b* で割った余りに一致する。

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

## 相補的なキャリー付き乗算
lag-*r* MWC の周期を決めるのは、通常 *p*=*ab**r*&nbsp;&minus;&nbsp;1 が素数となる乗数 *a*  の選び方である。もし *p* が安全素数なら、*b* の位数は *p*&nbsp;&minus;&nbsp;1 か (*p*&nbsp;&minus;&nbsp;1)/2 となる。ただ、*b* mod *p* の位数を求めるには *p*&nbsp;&minus;&nbsp;1 を因数分解する必要があるだろうが、これを因数分解するのは難しいだろう。

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

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

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

: *x*0, *x*1, *x*2, ..., *x**r*&minus;1, *c**r*&nbsp;&minus;&nbsp;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*&nbsp;&minus;&nbsp;1)&nbsp;&minus;&nbsp;*x* を使うのである。

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

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

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

*b*=232 の時、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*\cdot2327652 となり、109111 or 108798 or 108517 の *a* に対して約 109867 となる。

*b*=232 で *a* = 3636507990 の時、*p* = *ab*1359&nbsp;&minus;&nbsp;1 は安全素数となるので、MWC 列の周期は 3636507990\cdot243487 \approx1013101 となる。

*b* = 232 の時、CMWC 生成器で現在見つかっている中で最大に近い周期を持つものに、素数 *p* = 15455296*b*42658&nbsp;+&nbsp;1 を元にしたものがあり、*b* の位数は 241489*21365056 すなわち約 10410928 である。

## 関連項目
* 疑似乱数

## 参考文献
* 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
* 

## 外部リンク
*C++による実装（英語）
*Pythonによる実装（英語）