J. Plasma Fusion Res. Vol.96, No.7 (2020)３
６７‐３７１

講座 シミュレーションのための乱数入門

３．乱数発生の高速化と並列化技法
3. Techniques of Acceleration and Parallelizing the Random Number Generation
佐竹真介
SATAKE Shinsuke
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年３月１
７日）

現代のスーパーコンピュータを用いた大規模シミュレーションは，多数の計算ノードを用いた並列計算が標準的になっています．本章では，MPI 分散並列プログラムにおける疑似乱数の並列発生方法について，数値計算
ライブラリ KMATH̲RANDOM を例に解説します．また，乱数発生の高速化チューニングの実例や，物理乱数発生器を使う場合の高速化，並列化についてのアイデアについても紹介します．
Keywords:
pseudo random number, physical random number, parallel computing

・はじめに依存しますが，基本的に分散並列の各プロセスは異なる乱

本章では，大規模シミュレーションでの乱数の利用を想

数列を利用しなければならないということです．例えば，定し，乱数生成を効率的に行うためのプログラミング上の技

分子動力学（MD）シミュレーションでタンパク質の折り畳法について解説します．現代の大規模シミュレーションではみ構造を調べる計算を，１０プロセスの分散並列で行うとし

特に，MPI や OpenMP といった分散並列，スレッド並列をましょう．異なる初期配置を持ったタンパク質分子を，乱駆使したプログラムを利用することが一般的になっていま

数を利用して多数用意する際に，それぞれのプロセスで同

すから，そのような並列化プログラムの中で乱数を利用すじ疑似乱数発生法，例えば MT 法で，第２章で説明した

る上での留意点や工夫について紹介したいと思います．な

seed の値も同じに取ってしまったら，１０個のプロセスで全

お，本 章 で は 疑 似 乱 数 発 生 法 と し て 第２章 で 紹 介 し た

く同じ計算を開始してしまうことになり，並列化の意味がMersenne-Twister（メルセンヌ‐ツイスター）法（MT 法）を 対ありません．そんなことはわかりきったことであり，seed

象として説明しますが，物理乱数器の利用を想定した場合の値をプロセス毎に変えればよい，と考える人もいるかもの並列化や高速化に関するアイデアも紹介したいと思います．

知れませんが，それはお勧めできません．疑似乱数の seed

まず乱数の並列化の説明に先立って，MPI による分散並の値は，利用する乱数列の開始点を一周期のどこに取るか

列と，OpenMP（あるいはコンパイラの自動並列化機能にを決めるものだと２．
１節で説明しましたが，例えば seed のよる）スレッド並列の違いを明確にしておきましょう．分

値を１０個のプロセスででたらめに（人手で決める，時計を散並列は，計算を実行する CPU がそれぞれ物理的，あるい使う，他の疑似乱数を使うなど）与えても，それらの部分は論理的に分割された異なるメモリ領域に配列や変数を記

列の開始点が互いに十分離れているという保証は全くあり

憶し，シミュレーションに必要な計算を各CPUが分担してません．seed の選び方がたまたま悪くて２つのプロセスの行い，必要に応じて他のCPUと明示的な通信命令によって乱数列の開始点が数個しかずれてなかった，ということもデータをやり取りする並列計算です．一方，スレッド並列

考えづらいですが，起こりえない話ではありません．各プ

は，CPU の中にある複数のコアが，同一のメモリ上のデー

ロセスが参照する乱数の部分列が確実に重複しないようにタを使った演算を分担するやり方で，基本的には Fortran

選ぶことは，プロセス間で擬似乱数によるおかしな相関関の do loop や C 言語の for loop などの繰り返し計算を分割し係が生じる懸念を排除するために必要なことです．それでてそれぞれのコアが並列計算する方法です．この両者を組は，具体的にどのように乱数発生を並列していけばよいのみ合わせたハイブリッド並列も大規模シミュレーションでか見ていきましょう．

は一般的に行われています．以下，分散並列における並列・Mersenne-Twister 法による乱数生成の並列化

化の区切りを「プロセス」，スレッド並列におけるそれを分散並列計算における MT 法の乱数発生の並列化につい「スレッド」と呼んで区別します．
さて，並列シミュレーションにおいて乱数を利用する上ては，以前にもプラズマ・核融合学会誌の講座で紹介させで留意しなければならないことは，プログラムの詳細にもていただいたことがあります［１］．その時に紹介したのは，

National Institute for Fusion Science, Toki, GIFU 509-5292, Japanauthor’s e-mail:satake@nifs.ac.jp

３６７!2020 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.96, No.7 July 2020

Dynamic Creator と呼ばれる，MT 法の考案者の松本眞氏

コンパイラとコンパイルオプションの選択，NTL をインス

ら自身が開発した技法［２］で，ソースコードも公開されてトールしたディレクトリの指定です．random/arch/以下

い ま す［３］．Dynamic Creator で は，第２章 で 説 明 し た，にいくつかのコンピュータにおけるサンプルがあるので，

MT 法の漸化式の形を決めるパラメータのうち，ベクトル

それを参考にご自身の環境に合わせて設定します．なお，!!
"!
#$を変えることで，特性多項式が互いに素で，最大#コンパイルするには MPI 環境も必須です．また，MT 法の!

周期長がある Mersenne 素数 !に対し " !!となり，かつ

周期長 "!!!はここでコンパイルオプション “-DDSFMT_

高次元に均等分布する疑似乱数の漸化式の型を多数用意す

MEXP=p”で指定します．make コマンドを実行し，コンパイる，という方法を取ります．数学的に，特性多項式が互いルが完了すると，random ディレクトリ下の c/，c++/，f90/に素な漸化式同士は相関を持たないとされているため，

以下に libkm_random.a というそれぞれの言語用のライブ

Dynamic Creator は言わば互いに独立な乱数表を複数用意

ラリができるので，自分のプログラムで利用するにはこれすることに相当します．ただし，特性多項式が互いに素なをリンクしてコンパイルすることになります．

パラメータ #!!
"!
#$の探索はトライ＆エラー的になされる

次に，ジャンプファイルという N 個飛ばしの疑似乱数列

ため，数千並列の並列計算のための準備には非常に長い時の開始点の状態を記録したファイルを生成します．まず，間がかかります（筆者が十数年前に１０２４個用意するのに，

make ptool を実行すると，random/ptool/に km_rand_gen

ワークステーション１台を１か月以上稼働し続ける必要が_jump ができるので，これを実行します．

ありました）．

（例）km_rand_gen_jump -seed S1 S2 -max_ranks R -rand

そこで，超多並列計算のための実用的な MT 法の並列化_range M

として今回紹介したいのは，Jump Method と呼ばれるもの

ここで，-max_ranks の値Rは想定される並列計算の最大プです［４］．MT 法の特徴である長大な周期長を活かし，１

ロセス数，-rand_range の値 M はジャンプの間隔 N=2M

つの疑似乱数列の中の，遠く離れた複数の点を並列計算での指数部分を指定します．-seed の値 S1，S2 は，いわゆる乱数を利用する開始点として明示的に指定する方法です．

通常の意味での seed の値を S1，S1＋1，…，S2 と変えたもこれは，先に述べたような seed の値を単に変化させるのとのを（S2−S1＋1）個生成するということです．これは N 個は異なり，開始点間の距離を #"!個ずつ飛ばしに指定す飛ばしで R 個用意する開始点全体の起点が異なるジャンプ

ることができます．N 個飛ばしした後の漸化式の内部状態

ファイルを複数作成することを意味します．通常，１つのを並列プロセス数 P だけ事前に用意するには，原理的にはシミュレーションに使うジャンプファイルはどれか１つをN 個乱数を発生させては内部状態を記録する（第２章のサ

選べばよいですが，異なるジャンプファイルは，乱数をンプルプログラム に お け る”mt_rand.txt”の デ ー タに相使ったシミュレーションを複数回行って統計平均を取りた当）という操作を P 回繰り返せば可能です．ですが，P 個い場合などに利用します．試しに S1＝1，S2＝2，R＝200，の乱数列の部分区間が重ならないように N を十分大きく

M＝100で実行すると，jump/ディレクトリに file_00001，（１００兆など）取り，かつ P が数千となると，実用的な方法

00002 の２つのジャンプファイルのセットができます．ことは言えません．しかし，Jump Method に使われる SIMD-

れで，MT 法の疑似乱数に対し，2100 飛ばしの開始点の情報

oriented Fast Mersenne Twister（SFMT）という新しいタが２００個×２系列作られました．この計算は例えば最近のイプの MT 法［５］では，［６］に示される原理に基づいて N

Xeon CPU 上では数秒で終わってしまいます．これは逐次個飛ばし後の MT 法の漸化式の内部状態に一気にジャンプ

計算で 2100×200 個×２回，通常の MT 法で乱数を発生させできる計算ルーチンが用意されているので，それを利用す

るのとは比較にならないほど速いです．

ることで乱数列の遠く離れた先の開始点から始めるための

これで前準備は完了です．次にこのジャンプファイ

内部状態を高速に取得できます．SFMT は名前から推察で

ルを使って実際に並列計算で SFMT を利用する方法をきるように，高速化チューニングされた乱数発生ルーチン

見ていきましょう．第２章と同様に，サンプルプログラム

です．従って Jump Method の利用に限らず，逐次的に MT

rand_test.f90 と km_rand_wrapper.f90 を GitHub に挙げ

法で疑似乱数を発生させる用途としても，従来のバージョておきます［８］．ifort コンパイラを利用する場合のコンパンよりも SFMT を用いた方が高速実行が期待できます．

イル用 makefile のサンプルが用意してあるので，ご自分

さて，実際の利用法についてですが，ここでは［５］のオリの 環 境 に 合 わ せ て 修 正 し て make し，実 行 プ ロ グ ラ ム

ジナルのソースではなく，C，C++，Fortran９０での実装がkmrand_time を作成してください．実行方法ですが，この用意されている，理化学研究所の大規模並列数値計算技術

サンプルプログラムは MPI とスレッド並列のハイブリッ

研究チームが整備・公開しているライブラリ"KMATH̲

ドコードになっているので，環境変数 OMP_NUM_THREADS にRANDOM"［７］をベースに紹介したいと思います．このラ

スレッド並列数を設定し，また KMATH_RAND_JUMP_FILE_PATH

イブラリの利用許諾についてはマニュアルの記載事項を参に先ほど作成した file_00001 等のジャンプファイルが格照してください．まず前準備として，NTL という数値計算

納されているディレクトリを絶対 path で設定します．そしライブラリを別途ダウンロード，インストールする必要がて，”mpirun -np N kmrand_time”のような MPI 実行コマンありますが，［８］に説明があるので詳細は省きます．次に，ドで計算を開始してください．ここで N は試したい MPIKMATH̲RANDOM を 展 開 し て，random/Makefile.ma-

並列数（N はジャンプファイル作成時の-max_ranks の値 R

chine というファイルを編集します．必要なのは，利用する以下）です．
３６８

Lecture Note

3. Techniques of Acceleration and Parallelizing the Random Number Generation

S. Satake

このテストプログラムでは，rand_test.f90 の冒頭に指る場合，ある程度の量の乱数列に対する変換をスレッド並定された，１プロセス当たり ntest 個の乱数を，ntime 回発

列で一気に行った方が効率が良さそうです．更に，シミュ生させてその平均時間を計測します．そして，そのテスト

レーションの中で不定長の乱数を発生させるためにを nseed 個の異なるジャンプファイルに対して繰り返しまKMATH_Random_Vector を高頻度に呼ぶより，ある程度の大す．計測結果は"log.txt"に書き出されます．rand_test.f９０

きさの，固定長 nbuff の乱数列を一気に生成し，必要な分の主要部分を抜粋すると以下のようになっています．

だけ切り崩して使っていき，足りなくなったらまた一気に-----------------------------------------

作る，というやり方にします．井戸水を使う人は，使うた

call MPI_INIT(IERR)

びに必要な量を井戸まで汲みに行くのではなく，一定量水

call MPI_COMM_SIZE(MPI_COMM_WORLD,nmpi,ierr)

瓶にため込んでおいてそこから必要な量だけ柄杓ですくっcall MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)て使うのと同じように，ある一定量の乱数を貯めこんでか

call KMATH_Random_Init(handle1,MPI_COMM_WORLD, ierr)ら使う方法を，私は「水瓶方式」と呼んでいます．そのような一連の作業をしているのが grnd_dsfmt になります．
ここで，最も効率が良くなる nbuff を選ぶことで高速化

do k=1,nseed
call KMATH_Random_Seed(handle1, k , ierr)につながりますが，その値は計算機の環境によります．

do j=1,ntime

図１は，IFERC-CSC の ス ー パ ー コ ン ピ ュ ー タ JFRS1 上
で，8 MPI×5 スレッド並列で kmrand_time を様々な nbuff

call grnd_dsfmt(rnd1, ntest)と ntest に対して実行した際の，乱数１００万個当たりの平均

end do

発生時間を表しています．ntest が大きくても小さくても，

end donbuff=64000 の時に最も乱数発生時間が短くなりました．
call KMATH_Random_Finalize(handle1, ierr)

この例では，チューニングの効果はあまり大きくありませcall MPI_Finalize(IERR)

んが，他の計算機環境では，特に ntest が小さい時に nbuff--------------------------------の選び方で数十％の速度の差が出ることもありました．な

また，上のコードの流れは，以下のようになっています．

おこの水瓶方式ですが，grnd_dsfmt は MPI の各プロセスで１．MPI_INIT で MPI 並列計算を起動

独立に動くので，それぞれのプロセスで乱数の発生数やタ２．KMATH_Random_Init で KMATH_RANDOM の利用準備

イミングがずれていても問題ありません（並列化乱数発生

３．KMATH_Random_Seed で k 番目のジャンプファイルを読法とは当然そうあるべきものですが）
．また，
２．
２節で紹介

み込む

したような物理乱数を利用する際も，物理乱数のデバイス

４．grnd_dsfmt(rnd1, ntest)で配列 rnd1 に ntest 個の乱にプログラムからアクセスして乱数データを取得するオー数を発生させる

バーヘッドタイムを考慮すると，必要に応じて乱数を取得

５．KMATH_Random_Finalize で KMATH_RANDOM の利用終了

するよりも，水瓶方式を用いてある程度の量の乱数を一気

６．MPI_Finalize で MPI 並列計算の終了に取得する方が効率がよいと考えられます．

このプログラムでは nseed 個のジャンプファイルについ・物理乱数発生器を利用した並列化て，各プロセスで発生した乱数の最初の ntest 個を，”rand_s.txt”に全プロセス分をまとめて書き出します（全部をさて物理乱数の話が出たので，次に分散並列計算を行い出力すると膨大になるので 1/100 に縮めてあります）．プ
ロセス毎，seed 毎に異なる乱数列が生成されることが確認できると思います．

・KMATH̲RANDOM の高速化チューニング
さて，このサンプルプログラムで発生速度が計測できる
ので，次に高速化について考えましょう．乱数の発生は
grnd_dsfmt というサブルーチンで行っていますが，これはkm_rand_wrapper.f90 の 中 に 定 義 さ れ て い ま す．実 は，KMATH_RANDOM が 用 意 す る 乱 数 発 生 ル ー チ ン KMATH_Random_Vector は，高速化のために区間 !
!
"#"の８バイト実数を，３８６個以上の偶数個発生するという少し変わった作りになっています．なので，１回の必要な発生数が少ない場
合や，奇数の場合の対処法が必要です．また，通常シミュ
レーションで使われる乱数は !
!
!
""区間に発生させた一様乱数を利用することが前提になっていることが多いですか

図１

ら，KMATH_Random_Vector で発生させた乱数を!
!
!
""区間に
変換する操作が必要です．スレッド並列が有効になってい３６９

水瓶方式でバッファの大きさ nbuff を変えた場合の，乱数の発生時間の変化（百万個当たり）．nbuff ＝ 64000 で極小値となります．

Journal of Plasma and Fusion Research Vol.96, No.7 July 2020

たいが物理乱数発生器が１つしかない場合，乱数の並列化をどうするか考えてみましょう．図２のように，計算の１ステップが Calc. A, B, C, …といったルーチンからなり，Calc.B において 乱 数 が 必 要 で あ る 場 合 を 想 定しま
す．この時，例えば量子乱数発生器（QRNG）が rank0 の
MPI プロセスが走るノードにのみあるとします．考えられる並列化の方法としては，rank0 は物理乱数の取得と，それを他のプロセスに分配することに専念し，他のプロセス
図２

が計算を担当するというやり方です．この際，「１ステッ

MPI 並列プログラムで量子乱数発生器（QRNG）が１つしかない場合の「源泉かけ流し方式」におけるプログラムの流れの模式図．rank0 のプロ セ ス は 乱 数 生 成 に 専 念 し，Calc.A が終わった他の MPI プロセスは乱数を受信し，乱数を利用する Calc.B に備えます．

プ当たり必要とされる乱数の最大個数（プロセス毎に異なる場合も含め）」×「計算に従事するプロセス数」の乱数を
rank0 のプロセスが乱数発生器から取得し，乱数が必要になる Calc.B の手前の Calc.A ルーチンの間に MPI_ISCATTER（あるいは ISEND/IRECV の組み合わせ）（＊）による非同期

integer :: n,nmpi

通信によって各プロセスに乱数を配信します．このやり方

real(kind=8) :: rand(n*nmpi), rand_local(n)では，乱数を利用しないルーチンにかかる時間が，必要な!$OMP PARALLEL

乱数の発生にかかる時間に比べて長ければ，ほぼ待ち時間!$OMP MASTER ! only master thread does the follow-

なしで並列計算が進められるため，計算の高速化にもつな

ing process

がります．各プロセスでは受け取った乱数を１ステップでif (myrank==0) then !assume MPI rank 0 has the

全部利用しない場合もありますが，コードの実装を単純に

QRNG board

するために，利用しなかった乱数は捨てて，次のステップ

call QRNG_GEN(rand,n*nmpi) !generate random

で新しい乱数を rank0 から受けとる形にします．この方法

numbers

を私は「源泉かけ流し方式」と名付けました．温泉で源泉

end ifから汲んだお湯を，下流の各浴槽の使用量にかかわらず常

call MPI_SCATTER(rand,n,mpi_real8,rand_local,n, &に溢れるぐらい十分な量を流し込み続け，使わなかったお&湯は再利用せずに捨ててしまう方式をイメージしたもので

mpi_real8,0,mpi_comm_world,ierr)!scatter random numbers ( n per each MPI rank)

す．Master-Slave 形式の並列計算（パラメータ並列計算!$ OMP END MASTER

で，Master プロセスが計算全体の管理を行い，手の空いた!$ OMP DO schedule (dynamic)

Slave プロセスに次の計算パラメータを順次伝える形式）!thread parallel

do j=1,mで初期条件に乱数が必要な場合にも，この源泉かけ流し方…. ! (Calculation A)

式で Slave プロセスのどれかが計算を終えるまでに，Mas-

end doter プロセスが次に使う乱数をため込んで準備しておく，!$ OMP END DOという形での利用が可能です．また，乱数発生法が物理乱!$ OMP END PARALLEL

数に限定されるわけでもなく，疑似乱数を使った並列計算!... Both communication and calculation A are com-

をしたいが，先に説明した Jump Method を準備するのが面

pleted here.

倒という場合にも，一つの MT 法で発生させた疑似乱数列

call sub_calc_B (n, a, b, c, •••, rand_local)

を源泉かけ流し方式で各プロセスに分配して使うことも可! Calc B which uses the random numbers

能です．

call sub_calc_C (.....)

次に MPI×OpenMP のハイブリッド並列の場合を考え•••ます．もし図２の Calc.A の計算時間が，QRNG で全 MPI

各 MPI プロセスには，0〜k-1 の k 個のスレッドがありまプロセスが利用する乱数を発生する時間より十分長い場

す．このうちの０番スレッドが Master スレッドであり，合，rank0 を乱数生成だけに専念させるのは CPU コアを無Master スレッドだけが，$OMP MASTER 〜 END MASTER の区間

駄に遊ばせていることになりますので，rank0 には計算にを実行します．Rank0 の MPI プロセスの Master スレッドはも従事してもらいたいです．そこで，Calc.A の並列実行の乱数を発生器から取得し，MPI_SCATTER で n 個ずつ各プロセ

仕方を図３のように変更してみましょう．この計算のプロ

ス に 分 配 し ま す．一 方，Rank0 以 外 の MPI プ ロ セ ス のグラムのイメージは以下のようになります（GitHub［８］にMaster スレッドは，乱数の受信に携わります．その間，も gensen_hybrid_img.f90 として上げてあります）．

全ての MPI プロセスにおいて，Master スレッド以外の（＊）図２では簡単のため MPI_ISCATTER を使っていますが，これでは rank0 自身にも乱数を分配することになります．この無駄を省くには，総プロセス数を nmpi，１プロセス当たりに分配する乱数の数を n とした場合，まず n×(nmpi−1)個の物理乱数を rank0
で生成し，それを配列 rand(n*nmpi)の n+1 番目以降に保存します．配列 rand の先頭の n 個の値は与えなくてよいです．そして，
この配列 rand を MPI_ISCATTER で配信すれば，rank0 自身には無意味なデータが渡り，rank1 以降に生成した乱数が n 個ずつ分配
されます．

３７０

Lecture Note

3. Techniques of Acceleration and Parallelizing the Random Number Generation

S. Satake

は，MT 法の内部状態と水瓶方式の バ ッ フ ァ の 中 身 を
save，load す る た め の サ ブ ル ー チ ン save(load)_grnd_
dsfmt が用意してあります．その利用法のサンプル cont_test.f90 が GitHuB に上げてありますのでご参照ください［９］．は じ め に，input_cont.txt 中 の njob を０に し て
kmrand_cont を MPI 並列で実行すると，インデックスがiseed のジャンプファイルを使って新規に ntest 個の疑似乱数を各 MPI プロセスが生成し，”out_r**R_s**I_i000”に書き出します．ここで**R は MPI のランク，**I は選択した
ジャンプファイルの seed 番号です．また，乱数の継続用のファイル ”rand_cont_s**I_i000” ，”grnd_buff_ s**I_i000”も作られます．これらはそれぞれ，MT 法の内部状態と水瓶方式のバッファの状態を全 MPI ランク分まとめて書き出したものになって い ま す．次 に，njob=1 と し て 再 度
kmrand_cont を実行すると，上の継続用ファイル２つを読

み込み，疑似乱数列の続きを ntest 個作ります．以降，njob=2,3,…と す る ご と に，”rand_cont_s**I_i**N”,”grnd_
buff_ s**I_i**N” (N=njob-1)のデータを読んで疑似乱数列を一つ前のジョブの続きから生成します．ntest の数をいろいろ変えて実行し，各 MPI ランクが作る疑似乱数がそ
図３

れぞれ異なる系列の乱数列で，かつ各々はジョブの区切り

MPI-OpenMP ハイブリッドコードにおける源泉かけ流し
方式の模式図．各プロセスの Master スレッドは，乱数の送受信が終わったところから Calc.A のスレッド並列計算に途
中参加します．

を跨いで連続した乱数列になることを確認してください．
駆け足になりましたが，これで KMATH_RANDOM ライブラリを利用した，並列化 Mersenne-Twister 法の利用と高速化

1〜k-1 番スレッドは，Calc. A の Do loop を並列実行し始と，物理乱数発生器を並列プログラムで利用する方法の解めています．$OMP DO に SCHEDULE (dynamic)が指定されて説はおしまいです．本章の内容とサンプルプログラムが，いるため，Master スレ ッ ド が 乱 数 の 生 成 と 分 配 を終え乱数を利用した多並列シミュレーションコードを書く方々て，まだ他のスレッドが Calc.A を終えていなければ，図３の一助になれば幸いです．

に示すように Master スレッドは途中から Calc. A の並列計算に参加することになります．このようにすれば，源泉か

参考文献

け流し方式でも乱数発生器を持つプロセスも計算に参加［１］講座「核融合プラズマシミュレーションの技法 ５．粒子シミュレーションのコーディング技法」プラズマ・核融合学会誌 89, 245 (2013).
［２］M. Matsumoto and T. Nishimura, "Dynamic Creation of
Pseudorandom Number Generators", Monte Carlo andQuasi-Monte Carlo Methods 1998, Springer, 2000, pp. 5669.
［３］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
DC/dc.html［４］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
SFMT/JUMP/index-jp.html［５］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
SFMT/index-jp.html［６］H. Haramoto et al., "A Fast Jump Ahead Algorithm for
Linear Recurrences in a Polynomial Space", Lecture
Notes in Computer Science, vol. 5203, Springer, Berlin,
Heidelberg (2008).
［７］https : / / www. r-ccs. riken. jp / labs / lpnctrt / projects /kmath-random/［８］https://www.shoup.net/ntl/index.html［９］https://github.com/satakeshinsuke/random̲number/blob/master/Chapter3.tar

し，かつ乱数の分配のMPI通信時間を Calc.A の演算に隠ぺいすることが可能です．もちろん，この方式が効率よく稼働するには MPI 並列数とスレッド並列数を適切に選び，Calc.A にかかる時間より乱数生成時間が短くなるようにしなければならないので，どんなプログラム・計算環境にも使える万能な方法ではないことに留意してください．

・KMATH̲RANDOM における継続計算の方法
最後に，KMATH_RANDOM を使った計算を複数回継続する方法を説明します．第２章でもみたように，１回のジョブの最後に MT 法の内部状態をファイルに書き出し，継続ジョ
ブの最初にそれを読み直せばよいのですが，MPI 並列プロ
グラムではそれを各プロセスが行う必要があります．更に，grnd_dsfmt のように水瓶方式を使っている場合，「水瓶」の中身も保存し，継続計算に引き継がなければなりま
せん．そうでなければ，例えば１００ステップを３回つないだ乱数を使った計算と，１５０ステップを２回つないだ計算とで途中から利用する乱数がずれてしまいます．
それでは，サンプルプログラムを使って，継続の仕方を具体的に説明していきましょう．km_rand_wrapper.f90 に
３７１