本論文は,二元体F2 上の線形擬似乱数発生法に対して,その評価方法の一つとして知られている高次元均等分布性における均等分布の次元計算アルゴリズムの高速化について論じる.近年,大規模なコンピュータ・シミュレーションを行う際,並列化・分散化が行われるようになり,複数のCPU上に,異なるパラメータセットを有する擬似乱数発生法を割り当てて使用する機会が増加している.そこで,非常に多数のパラメータ探索を行い,それぞれのパラメータに対して,乱数性を評価する必要が生じる.既存の均等分布の次元計算アルゴリズムは計算時間を要し,探索のボトルネックとなるため,高速化が望まれていた.

本論文では,F2-線形擬似乱数発生法と呼ばれる次のモデルを扱う.状態空間をS := Fp2,状態遷移関数をf : S → S,出力の集合をO := Fw2 (wはコンピュータのワードサイズ),出力関数をo : S → Oとおく.fとoは高速に計算可能なF2-線形写像を熟考の上に選ぶ.初期状態s0 ∈ S が与えられたとき,単位時間ごとに漸化式

si+1=f(si) (i=0; 1; 2;...)

を計算して,出力列

o(s0); o(s1); o(s2); ... ∈ O

を得るものとする.これらの出力を符号なしw ビット二進整数と同一視し,擬似乱数として用いる.

fとoを決める際には,擬似乱数発生法の一様性の評価規準として知られている高次元均等分布性が用いられる.この概念は,出力列の上位vビットのみに着目して得られる,vビット精度の均等分布の次元k(v)の大きさに基づく.k(v)は次に述べる二元体係数形式的べき級数体K := F2((t-1))上の簡約基底法を用いて,具体的に計算することができる.状態から出力の上位vビットのみの出力を与える出力関数をov : S o→O→Fv2とし,上位v ビット出力列ov(s0); ov(s1); ov(s2);...∈ Fv2に対して,以下のF2-線形生成母関数xv:S→Kvを考える:

このχv(s0)とv 次元単位ベクトルe1,...,evを生成集合とするF2[t] 上の加群

を考えると,v 次元F2[t]-格子となる.ここで,F2[t]-格子の基底として,tの次数によって定義されたウルトラノルムに関して最短のものを簡約基底と呼ぶ.状態遷移関数fの特性多項式が既約の場合,Λvの簡約基底における最長ベクトルのノルムがk(v)を与えることが知られている.

擬似乱数発生法の評価のためには,すべてのk(v) (v=1,...,w)を計算することが必要となる.Couture-L'Ecuyer(2000)は,Λvの双対格子とLenstra(1985)による格子基底簡約アルゴリズムに基づいて,k(v-1)を計算する際に求めた簡約基底を利用してk(v)を効率的に計算する方法を提案した.この方法は既存の最も高速な計算法として知られていたが,次数の高い多項式を成分に持つベクトルを扱う必要があり,改善の余地があった.

そこで,本論文では,第I部と第II部とに分けて,F2-線形擬似乱数発生法に対して,均等分布の次元k(v)の計算アルゴリズムの高速化を提案する.

第I部では,双対格子に変換せず,Λvの簡約基底を直接求める方針を取った.そして,従来用いられていたLenstra アルゴリズムの代わりに,Schmidt(1991)による格子基底簡約アルゴリズムを採用し,入力が生成集合に対応するように修正し,SGR アルゴリズムと名づけた.これは,双対格子の場合の計算順序とは逆に,k(v)をv=w;w-1,...,1の順に帰納的に計算するためである.すなわち,p: Kv → Kv-1をv 番目の座標を削除する射影とすると,p(Λv)=Λv-1 が示されている(Couture-'Ecuyer, 2000).そこで,Λvの簡約基底を求めておき,各々のベクトルにpを施すとΛv-1の生成集合が得られる. 入力が基底に制限されているLenstra アルゴリズムでは適用できないが,SGR アルゴリズムではこの生成集合から簡約基底を求めることが可能となる.この射影を用いた計算法をインダクティブ・プロジェクションと名付けた.

更に,Λvの点を状態空間Sの元を用いて表現する方法を提案した.状態遷移関数fの特性多項式が既約で,χv が非零であると仮定する.このとき,F2-ベクトル空間の直和としてΛv=F2[t]v +χv(S)と分解でき,任意の格子点ω∈ Λvは,ω=poly + χv(s) なる一意的な表現(poly; s) ∈ F2[t]v×Sを有する.この(poly; s)を!の状態表現と呼ぶ.SGR アルゴリズムで必要となる格子点同士の和は状態表現の和に,t倍の操作ω ・ tは(poly・ t + ov(s); f(s))と状態表現され,状態表現内部の代数的な演算に帰着可能であることを示した.この表現により,格子点の成分に現れるべき級数を回避し,結果として,格子簡約アルゴリズム内で扱うΛvの生成集合に対して,各ベクトルのメモリ使用量をdim(S) ビットに抑えることに成功した.

以上,SGR アルゴリズム,インダクティブ・プロジェクション,状態表現の一連の手法をSIS 法と名づけ,ビット演算量の最悪値解析を行った.その結果,Couture-L'Ecuyer(2000) よりも弱い仮定の下で,すべてのk(v) (v=w,...,1)を求めるために要するビット演算量の上限

を得た.Couture-L'Ecuyer(2000)の方法と比較し,dim(S)2の項におけるwのオーダーが減少しており,実用上,dim(S) がwに対して非常に大きいため,SIS 法の方が高速であると予想される.実際,上述のアルゴリズムを代表的なF2-線形擬似乱数発生法であるMT 法(松本-西村, 1997)とWELL 法(Panneton-L'Ecuyer-松本, 2006)に適用し,計算機実験を行った.すべてのk(v)を求めるのかかったCPU時間を測定したところ,既存の双対格子による方法よりも10倍程度の高速化に成功した.

第II部では,SIS 法の改良を行い,更なる高速化を図った.SIS 法において採用したSGR アルゴリズムは,係数ベクトル間のF2-線形関係式を求めるために,行列の掃き出し計算を繰り返す必要が生じる.格子の次元が大きくなると,この部分に起因する計算量の増加は無視できない.そこで,掃き出し計算を回避するために,SGR アルゴリズムの代わりに,Mulders-Storjohann(2003)による簡約アルゴリズムを採用することを提案した.特に,Λvの生成集合の簡約では,Wang-Zhu-Pei(2004)による並べ替えの手法を用いることにより,係数ベクトルを一定の条件(三角条件)に保ったまま,簡約を行うことが可能となる.これらの方法を使いやすい形に修正し,ピボット格子簡約アルゴリズムとして与え,インダクティブ・プロジェクションおよび状態表現と組み合わせ,PIS 法と名付けた.

PIS 法の計算量として,すべてのk(v) (v=w,...,1)の計算に要するビット演算量の上限

を導出し,SIS 法におけるdim(S)1の係数並びにdim(S)を含まない項におけるwのオーダーが共に1 次減少することを示した.これは,w が増加した場合にPIS 法の方が有効であることを示唆している.実際,SIS 法とPIS 法による計算機実験を行い,MT 法およびWELL 法に対して,すべてのk(v)を求めるのにかかったCPU時間を計測した.dim(S)=19937の場合,w=32の例ではSIS法の約1:5倍,w=64の場合には約3倍高速化された.

一方,異なった方向性の高速化技法として,初期状態s0 ∈ Sの取り方について考察した.Panneton-L'Ecuyer-松本(2006)は,MT 法の初期状態に,極端に0の多い零超過状態を与えると,長時間に渡って,出力の0と1のバランスに偏りが見られるMT 法の欠点を述べている.本論文では,この現象を逆手にとり,零超過状態s0を用いてΛvの生成集合の一元χv(s0)を与えると,MT 法の格子簡約に現れる状態表現に1 が疎な状態が続くことから,非常に高速に簡約基底が求まることを発見した.PIS 法と組み合わせ,実験を行ったところ,MT 法においては,通常の初期化を用いたSIS 法に比べて10倍以上の高速化を得た.また,零超過状態による初期化の有無にかかわらず,いずれの場合も,PIS 法はSIS 法よりも高速となった.

このアルゴリズムはMT型発生法のパラメータ探索ソフトウェアDynamic Creator 及びMTGPDCに組み込まれ,ホームページから配布されている.

擬似乱数発生法とは、乱数のように見える数列を、決定性のアルゴリズムにより計算機内で高速かつ再現性のあるように生成する方法を指す。計算機による科学シミュレーションにおいて、確率的な要素を含む部分には擬似乱数は欠かせない。また、近年の計算機の並列化に伴い、各計算処理単位システムごとに相異なる擬似乱数発生器を配置することが多くなり、高速で高品質な擬似乱数発生法を、非常に多数用意する必要性が増えてきた。

擬似乱数発生法として現在最も有力なものは、二元体F2 上のベクトル空間を状態空間とし、線形な状態遷移関数と出力関数により状態を変移させつつ出力列を得る、F2 線形擬似乱数発生法である。擬似乱数の品質を測る重要な評価基準としては、v ビット精度での高次元均等分布性が有力である。擬似乱数の出力列を[0,1] 実数列に正規化してk個ずつの組にわけ、k 次元単位立方体内に擬似ランダムな点を一周期に渡って生成した際、各座標に関してv ビット精度で一様に点が分布しているとき、この擬似乱数列はv ビット精度でk 次元一様分布しているという。そのようなkの最大値をk(v)で表す。k(v) が大きいほど擬似乱数の品質は良いと考えられるが、状態空間の次元dに対してk(v) ≦ d/v なる上限がある。

生成する擬似乱数自体の精度をw ビットとする。1 ≦ v ≦ wの各vについてk(v) が上限に近いとき、擬似乱数は(作為的なものを除けば)通常の統計的検定に合格することが経験的に知られており、この基準に基づくパラメータ探索が有効である。しかし、d が大きいとき、k(v)の計算には大きな時間がかかり、パラメータ検索におけるボトルネックとなっていた。本研究以前の最速アルゴリズムは次の通り:v ビット出力の無限列に対応するベクトル値の生成母関数を考えたとき、初期状態として全ての状態を動かして得られる母関数の集合が多項式係数格子を生成するが、その簡約基底のノルムの最大値がk(v)を与える(Couture-zuka-'Ecuyer, 1993)。それを、双対格子により計算する(Couture-L'Ecuyer, 2000)。

これに対し、本論文では、(1) 状態空間を用いて格子点を表現しメモリ使用量とビット演算数を減らす、(2)w 次元空間における簡約基底から、順番にw-1 次元、w-2 次元の簡約基底を計算することで、前計算を有効利用する、(3)Mulders-trjohann(2003)による高速簡約基底計算を利用する、という三つの工夫により、k(v) (1 ≦ v ≦ w)の計算に要するビット演算の回数を2wd2 + 12w2d + 12w2(w + 1)に減らすことに成功した。これは、前述のCouture-L'Ecuyer 双対格子法の約1/wである。実際、計算機実験により速度を測定し、w=32, d=19937では20倍程度の高速化になることを確認している。このアルゴリズムは、F2 線形擬似乱数発生法のパラメータを自動的に探索するソフトウェアであるDynamic Creator 及びMTGPDCに斎藤睦夫広島大学助教の手により組み込まれ、ホームページから配布されており、世界的に利用が進んでいる。この研究結果の一部は米国数学会出版のMathematics of omputationに出版されており、残りは投稿中である。

以上によって、論文提出者 原瀬 晋は、博士(数理科学)の学位を受けるにふさわしい十分な資格があると認める。