キャリー付き乗算のソースを表示
←
キャリー付き乗算
ナビゲーションに移動
検索に移動
あなたには「このページの編集」を行う権限がありません。理由は以下の通りです:
この操作は、次のグループに属する利用者のみが実行できます:
登録利用者
。
このページのソースの閲覧やコピーができます。
'''キャリー付き乗算''' (キャリーつきじょうさん、{{lang-en-short|multiply-with-carry, MWC}}) は、George Marsagliaにより開発された整数疑似乱数生成用の手法である。乱数種には2〜数千の値を必要とする。MWC の主な長所は、単純な整数演算からなっており非常に高速に動作するという点と、2<sup>60</sup> - 2<sup>2000000</sup> という非常に長周期であるという点である。 ==基本理論== MWC 列は ''b'' を法とする剰余からなる。''b'' = 2<sup>32</sup> とするのが普通だが、これはコンピュータで扱う整数が通常そのようになっているからである。ただ、''b'' = 2<sup>32</sup> − 1とすることもある。これは、2<sup>32</sup> − 1 を法とする演算は 2<sup>32</sup> を法とする演算を少し変更するだけで済み、さらに、''b'' = 2<sup>32</sup> の MWC の理論にはいくつか厄介な問題があるが、''b'' = 2<sup>32</sup> − 1 ではそれを回避できるからである。 その最も一般的な形では、'''lag-r''' MWC 生成器は'''基数''' ''b''、'''乗数''' ''a''、そして ''r + 1'' 個の'''乱数種'''を必要とする。乱数種は ''r'' 個の ''b'' の剰余 :''x''<sub>0</sub>, ''x''<sub>1</sub>, ''x''<sub>2 </sub>,..., ''x''<sub>''r''−1</sub> と最初の'''キャリー''' ''c''<sub>''r''−1</sub> < ''a'' である。 そして、lag-''r'' MWC 列は : <math> 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</math> で定義される ''x''<sub>''n''</sub>, ''c''<sub>''n''</sub> のペアの数列であり、MWC 生成器の出力は以下の ''x'' の列となる。 :''x''<sub>''r''</sub> , ''x''<sub>''r''+1</sub> , ''x''<sub>''r''+2</sub>, ... lag-''r'' MWC 生成器の周期は ''ab<sup>r</sup>'' − 1 を法とする数の乗法群における ''b'' の位数である。通例、''a'' は ''b'' の位数を大きくできるよう、''p = ab<sup>r</sup>'' − 1 が素数になるように選ぶ。''b'' = 2<sup>32</sup> は ''p = ab<sup>r</sup>'' − 1 の原始根とはならないので、''b'' = 2<sup>32</sup> を基数とした MWC 生成器の周期は、MWC の最大周期 ''p = ab<sup>r</sup>'' − 2 になることはない。これが、''b'' = 2<sup>32</sup> − 1 の方が有利になる点の1つである。 Couture と l'Ecuyer (1997) により、MWC 生成器には最上位ビットが少し偏っているという理論的な問題があると指摘された。ただし、この問題は相補的キャリー付き乗算により解決されている。「相補的 MWC を使えば、最上位ビットは均一に出ることが分かるだろう。つまり、全周期中で0と1とが同等の頻度で現れ、その傾向も MWC 生成器間に関連性はない。」彼らはビットの偏り具合についてそれ以上詳しくは語っていないようである。相補的 MWC 生成器は計算時間が若干増えるため、実装の要求によってどちらを使うか決めるといいだろう。 ==線形合同法との比較== [[線形合同法]]は通常以下で定義される。 :<math>x_n=ax_{n-1}+c\ \bmod\,2^{32}</math> ほとんどの演算プロセッサでは乗数 ''a'' と現在の ''x'' は32ビットレジスタに置け、積は64ビットの連結レジスタ上に求められる。積の下位32ビット、すなわち :<math>a\times x\ \bmod\,2^{32}</math> は、右側のレジスタを参照するだけで得られる。同様に、これに32ビットの ''c'' を足してやれば自然に ''ax''+''c'' mod 2<sup>32</sup> が求まる。もし ''a'' mod 8 が 1 か 5 で ''c'' が奇数なら、2<sup>32</sup> を法とする合同数列の周期は 2<sup>32</sup> となる。 一方 lag-1 MWC は、線形合同法とほぼ同じ演算を使うだけで約 2<sup>62</sup> の周期を実現する。違うのは、MWC では64ビット積の上半分も利用している点である。これは '''次の''' キャリー ''c'' として使われる。一方、線形合同法では ''c'' は固定値である。''ax''+''c'' を64ビット値として求め、上半分を次の ''c'' として使い、下半分を ''x'' として使う訳である。 乗数 ''a'' が与えられれば、''x'', ''c'' のペアを入力として、新たなペアが作成される。 :<math>x\leftarrow (ax+c)\,\bmod\,2^{32},\ \ c\leftarrow \left\lfloor\frac{ax+c}{2^{32}}\right\rfloor</math> もし ''x'' と ''c'' が両方とも 0 でなければ、MWC 列の周期は ''ab'' − 1 を法とする剰余の乗法群における ''b'' = 2<sup>32</sup> の位数、すなわち、''b''<sup>32''n''</sup> = 1 mod (''ab'' − 1) となる最小の ''n'' になる。もし 28〜31ビットの ''a'' を ''ab'' − 1 が[[安全素数]]となるよう選ぶと、''ab'' − 1 と ''ab''/2 − 1 は両方とも素数になり、周期は ''ab''/2 − 1 となり、2<sup>62</sup> に近づく。実際にはとりうる32ビット値のペア ''x'', ''c'' の数となるだろう。 以下は、上記の安全素数の条件を満たす ''a'' の最大値の例である。 {| class="wikitable" border="1" |- ! ''a'' のビット数 ! ''b'' ! ''ab''-1 が安全素数となる ''a'' の最大値 ! 周期 |- | 15 | 2<sup>16</sup> | 31,743 | 1,040,154,623 |- | 16 | 2<sup>16</sup> | 64,545 | 2,115,010,559 |- | 31 | 2<sup>32</sup> | 2,147,483,085 | 4,611,684,809,394,094,079 |- | 32 | 2<sup>32</sup> | 4,294,966,893 | <sup>*</sup> |} 次に、10で割った余りという単純な例を使って線形合同法と MWC の比較を行う。ここでは''レジスタ''は10進数1桁のサイズであるとし、連結レジスタは10進数2桁であるとする。 <math>x_0=1</math> から始めると、合同数列 :<math>x_n=7x_{n-1}+3\,\bmod\,10</math> は連結レジスタ上で :<math>10,03,24,31,10,03,24,31,10,\ldots</math> の数列を持ち、''x'' の出力列（右のレジスタ）の周期は4となる。 :<math>0,3,4,1,0,3,4,1,0,3,4,1,\ldots</math> <math>x_0=1</math>, <math>c_0=3</math> から始めると、MWC 列 :<math>x_n=(7x_{n-1}+c_{n-1})\,\bmod\,10,\ c_n=\left\lfloor\frac{7x_{n-1}+c_{n-1}}{10}\right\rfloor,</math> は連結レジスタ上で :<math>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,...</math> の数列を持ち、''x'' の出力列（右のレジスタ）の周期は22となる。 :<math>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,...</math> ''x'' の繰り返し部分を逆順にすると :<math>449275\cdots97101\,449275\cdots9710144\cdots</math> となり、これは ''a''=7, ''b''=10, ''j=31'' とした時の ''j''/(''ab''−1) の値 :<math>\frac{31}{69}=.4492753623188405797101\,4492753623\ldots</math> の小数部と一致する点に注目したい。 これは一般的にもなりたち、lag-''r'' MWC により生成された ''x'' の列 :<math> x_n=(ax_{n-r}+c_{n-1})\bmod\,b\,,\ \ c_n=\left\lfloor\frac{ax_{n-r}+c_{n-1}}{b}\right\rfloor,</math> は、逆順にすると、ある 0 < ''j'' < ''ab''<sup>''r''</sup> において ''j''/(''ab''<sup>''r''</sup> − 1) の基数 ''b'' での小数部と一致する。 さらに、合同数列を :<math>x_n=7x_{n-1}\,\bmod\,69</math> で定義した場合、<math>x_0=34</math> から始めると以下の周期22の数列を得るが、 :<math>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</math> これを10で割った余りを求めると :<math>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</math> となり、上記の MWC 列に一致する。 これは一般的にも成り立ち（ただし、これは lag-1 MWC でしか成り立たないが）、初期値が <math>x_0</math>, <math>c_0</math> の lag-1 MWC 列 :<math> x_n=(ax_{n-1}+c_{n-1})\,\bmod b\,,\ \ c_n=\left\lfloor\frac{ax_{n-1}+c_{n-1}}{b}\right\rfloor</math> により得られる数列 <math>x_1,x_2,\ldots</math> は、合同数列 ''y''<sub>''n''</sub> = ''ay''<sub>''n'' − 1</sub> mod(''ab'' − 1) を ''b'' で割った余りに一致する。 初期値 ''y''<sub>0</sub> の選び方は、単に ''x'' のサイクルをずらすだけにすぎない。 ==相補的なキャリー付き乗算== lag-''r'' MWC の周期を決めるのは、通常 ''p''=''ab''<sup>''r''</sup> − 1 が素数となる乗数 ''a'' の選び方である。もし ''p'' が[[安全素数]]なら、''b'' の位数は ''p'' − 1 か (''p'' − 1)/2 となる。ただ、''b'' mod ''p'' の位数を求めるには ''p'' − 1 を因数分解する必要があるだろうが、これを因数分解するのは難しいだろう。 しかし、''p'' = ''ab''<sup>''r''</sup> + 1 の形の素数であれば、''p'' − 1 = ''ab''<sup>''r''</sup> は簡単に因数分解できる。そのため、''p'' = ''ab''<sup>''r''</sup> + 1 を素数とする場合の ''b'' によって規定される MWC を使えば、MWC の周期を求めるのに必要な計算数論は著しく簡単になるだろう。 幸いな事に、MWC を少し変更するだけで、''ab''<sup>''r''</sup> + 1 の形の素数を作る事ができる。これを[[相補的なキャリー付き乗算]](CMWC)と呼ぶ。 必要な入力は lag-''r'' MWC と同じで、'''乗数''' ''a'' と '''基数''' ''b''、そして ''r'' + 1 個の乱数種 : ''x''<sub>0</sub>, ''x''<sub>1</sub>, ''x''<sub>2</sub>, ..., ''x''<sub>''r''−1</sub>, ''c''<sub>''r'' − 1</sub> である。新しいペア ''x'', ''c'' を生成する方法は、以下のように少し変更される。 :<math> 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</math> つまり、新しい ''x'' を作る際に、補数 (''b'' − 1) − ''x'' を使うのである。 CMWC 生成器により作られる数列 ''x'' の周期は ''ab''<sup>''r''</sup> + 1 を法とする剰余の乗法群における ''b'' の位数になり、''x'' を逆順にしたものは、ある 0 < ''j'' < ''ab''<sup>''r''</sup> において ''j''/(''ab''<sup>''r''</sup> + 1) の基数 ''b'' での小数部と一致する。 lag-''r'' CMWC を使うと、''r'' が 512, 1024, 2048 と大きくなっても、簡単に周期を求める事ができる。（''r'' を 2 の累乗にすると、''r'' 個の最新の ''x'' を保持する配列内の要素へのアクセスが若干簡単に（そして速く）なる。） 以下に、いくつかの例を示す。 ''b''=2<sup>32</sup> の時、lag-1024 CMWC : <math> 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</math> の周期は ''a''<math>\cdot</math>2<sup>327652</sup> となり、109111 or 108798 or 108517 の ''a'' に対して約 10<sup>9867</sup> となる。 ''b''=2<sup>32</sup> で ''a'' = 3636507990 の時、''p'' = ''ab''<sup>1359</sup> − 1 は[[安全素数]]となるので、MWC 列の周期は 3636507990<math>\cdot</math>2<sup>43487</sup> <math>\approx</math>10<sup>13101</sup> となる。 ''b'' = 2<sup>32</sup> の時、CMWC 生成器で現在見つかっている中で最大に近い周期を持つものに、素数 ''p'' = 15455296''b''<sup>42658</sup> + 1 を元にしたものがあり、''b'' の位数は 241489*2<sup>1365056</sup> すなわち約 10<sup>410928</sup> である。 ==関連項目== * [[疑似乱数]] ==参考文献== * 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 * {{Citation | last =Couture | first =Raymond | last2 =L'Ecuyer | first2 =Pierre | title =Distribution properties of Multiply-with-carry random number generators | journal =Mathematics of Computation | volume =66 | issue =218 | pages =591–607 | year =1997 | doi =10.1090/S0025-5718-97-00827-2}} ==外部リンク== *[http://www.agner.org/random/ C++による実装（英語）] *[http://code.activestate.com/recipes/576707/ Pythonによる実装（英語）] {{DEFAULTSORT:きやりいつきしようさん}} [[Category:乱数]] [[Category:数学に関する記事]]
このページで使用されているテンプレート:
テンプレート:Citation
(
ソースを閲覧
)
テンプレート:Lang-en-short
(
ソースを閲覧
)
キャリー付き乗算
に戻る。
ナビゲーション メニュー
個人用ツール
ログイン
名前空間
ページ
議論
日本語
表示
閲覧
ソースを閲覧
履歴表示
その他
検索
案内
メインページ
最近の更新
おまかせ表示
MediaWiki についてのヘルプ
特別ページ
ツール
リンク元
関連ページの更新状況
ページ情報