J. Plasma Fusion Res. Vol.96, No.6 (2020)２
９２‐２９９講座 シミュレーションのための乱数入門

２．乱数発生の原理
2. Random Number Generation
佐 竹 真 介，菅 野 龍 太 郎
SATAKE Shinsuke and KANNO Ryutaro
自然科学研究機構 核融合科学研究所（原稿受付：２
０
２
０年３月１
７日）

乱数は，発生法により大きく２つに分けることができ，１つは，決定論的な演算による疑似乱数，もう１つは，ランダムな自然現象を乱数源として利用する物理乱数です．乱数を利用するシミュレーションは，当然のことながら利用した乱数に，その計算結果が影響されます．そのため，利用する乱数について，事前に理解しておくことが大切です．本章では，まずは，疑似乱数，物理乱数それぞれの発生法について概説します．
Keywords:
pseudo random number, physical random number, quantum random number generator

２．
１ 疑似乱数の発生法

では我々がコンピュータプログラムで利用している乱数

コンピュータシミュレーションでは，様々な用途で乱数とは何なのでしょうか？本節ではまず，我々が通常プログが必要となります．必要な数が１０個程度であれば，人の手ラムで利用している「疑似乱数」の発生原理について説明

でサイコロを振るなり乱数表を引くなりすればよいです

します．疑似乱数とは決定論的なアルゴリズムから生成さが，乱数を使うシミュレーションで要求される乱数の数はれる「乱数のように見える」数列のことです．本節ではい

一般的に数万から数億個といった膨大な量になります．こくつかの代表的な疑似乱数の生成方法と，それらの利点・

れは乱数を利用するシミュレーションの多くが，ある種の欠点を紹介します．なお本節の内容は［１‐３］を参考に書か

ランダムな事象をプログラムで模擬し，その統計的な性質れています．一方２．
２節では，まだあまり一般的に普及してに基づいて解を求めるという考え方，言い換えれば，あるいませんが，疑似乱数とは異なり現実の確率的事象の観測

種のアンサンブル平均を数値計算で評価することによって値から乱数列を取得する「物理乱数」について紹介します．

解を得るという考え方に基づいて作られているからです．

また，疑似乱数の発生法の高速化・並列化といった数値的

一般的に，確率変数のアンサンブル平均の統計誤差は標本な技法や，疑似乱数の品質，つまり「乱数らしさ」をどう

数 ! に対し "#!
( に比例して小さくなります．しかし，も定義し，評価するかについてはそれぞれ第３章，第４章でし乱数と思って使っている数列に予期せぬ偏りや相関があ

詳しく解説します．

る場合，それを使っていくら多くの標本を用意したとして!線形合同法も，得られる答え（統計量）の信頼性は上がりません．よって，性質の良い乱数を効率よく発生させることは，乱数をさて，決定論的な計算機の演算から乱数のように見える

使うシミュレーションの規模が大きくなるほど，またより数列を生み出す方法ですが，一般的に使われる方法は全て高い精度の計算が要求されるほど重要になると言えるで何らかの漸化式の計算によって生成します．また，コン

しょう．

ピュータの内部では数字は２進数で表現されますから，こところで，我々が現在利用できるコンピュータは，内部

こで説明する疑似乱数とは４バイトや８バイトの２進数で的にはデジタルに表現される数値の四則演算と論理演算表現された整数だと考えて下さい．単純で古くからよく知（AND，OR など），ビット操作（２進数で表現された数値

ら れ た 疑 似 乱 数 は「線 形 合 同 法」と 呼 ば れ る も の で，を全体に右または左にシフトする，特定のビットの値を変

: %#"
!
!
#
) を疑似乱数列とした時，一般的に以下のよ'&%'更するなど）を組み合わせた手続きによって動いており，

うな漸化式で表現されるものを指します．

完全に決定論的な演算しか行えません．一方で，乱数とい

'
$&%"%
$
!'
"'
""##
!""*"#
!&"$
%
%"#
%
&'
%

うのはその名の通りランダムで予測不可能な数列というこ（１）とですから，決定論的な計算機のアルゴリズムからこれをこ こ で パ ラ メ ー タ #&，$，"は い ず れ も 非 負 の 整 数 で生成するというのはおかしな話にも聞こえます．

#&""，$""であり，$&%"は（１）式右辺を "で割った

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

２９２!2020 The Japan Society of Plasma
Science and Nuclear Fusion ResearchLecture Note

2. Random Number Generation

S. Satake and R. Kanno

余りを取る，という意味です．したがって疑似乱数 .
%はの 実 数 に 変 換 し た 場 合，(次 元 空 間 に"より小さい非負の整数になります．

，'#!
!
!
!
"
#
) と点を並べていく
+
+
+
*!
%
('""!
('"#!
(%
'""&&

これがどのような疑似乱数となるかはパラメータと，結晶のように規則的な格子状に点が並んでしまうこと$!
"&の選び方に強く依存します．疑似乱数のアルゴリ#&!
%を指します．例えば混合合同法の場合，２次元平面に並べズムを考える上で重要となる性質はその「統計的性質」，

るとこの結晶構造が現れます．これは係数%
#!!
$
!
"&をどの

「周期性」，「再現性」の３つです．統計的性質とは生成さように選んでも現れます．また，"は周期長が最大になるれた性質がいかに乱数っぽく振る舞うかのことで，詳しいように，疑似乱数のバイト長と同じ（４バイト整数なら

説明は第４章で行いますが，ここでは私たちが一般的な確"##$#とするなど）に取ることが一般的ですが，"が偶数

率的事象に対して期待するような振る舞いを，疑似乱数がの場合，一番下位のビットが必ず 0，1 が交互に現れる，つどれだけ上手く模擬できているかのことだと考えて下さまり偶数，奇数が交互に現れるといった性質もあります．

い．周期性についてはこの線形合同法に限らず，有限の有乱数を用いたシミュレーションの中には，このような規則

効桁数で演算されるコンピュータの漸化式による疑似乱数

性が計算結果に影響を与えてしまうものもあるため，疑似発生法には必ず固有の周期長がある，というのが重要な点

乱 数 と し て 望 ま し い 性 質 と は 言 え ま せ ん．周 期 長 がです．つまり，（１）式のような漸化式で数列を生成してい#$"$#!億というのも，１００万個の粒子を使ったモンテカル

くと，どこかでまた同じ数列の繰り返しになるということ

ロシミュレーションを２０００ステップ進めただけで１周期使です．また，決定論的なアルゴリズムなので，同じ初期条

い切ってしまうということで，現在の大規模計算において

件'
.
.
.
.
*!
!&!
!&
""!
!"!
%
%
%
%(を 与 え れ ば 同 じ 乱 数 列は十分に長い周期とは言えなくなってしまいました．

.
.
* が毎回得られる，つまり疑似乱数には再現性
""!
"#!
%
%

そこで，線形合同法とは異なる疑似乱数がい ろいろがあります．乱数なのに再現性があるというのもおかしくと 考 案 さ れ て き ま し た．こ こ で は そ の 代 表 例 と し て聞こえますが，同じ計算を同じ条件で走らせた場合に同じ

Tausworthe（トーズワース）の方法と，その発展形である

結果が得られるかどうかは，シミュレーションコードの開

Mersenne-Twister（メルセンヌ・ツイスター）法について

発時のバグ取りや検証のためには必要不可欠な性質です．

簡単に紹介します．

線形合同法のうち，特に &#!で，.
+-,"&
!.
""##
%
%%
!Tausworthe 法と い う 最 も 簡 単 な 形 を と っ た も の は 乗 積 合 同 法，.
+-,"&の形のものは混合合同法とも呼ば
%
!.
""##
%
%"$
れ，一つ前の値 .
%のみ記憶していればよいので，メモリや造を示さない疑似乱数として知られており，M 系列と呼ば計算速度に制限の多かった古い時代からよく使われた疑似れる疑似乱数の仲間の一つです．!%を -ビットの２進数を乱数です．なお，現代では混合合同法のことを線形合同法

表す横ベクトルとするとき，以下の漸化式

Tausworthe 法は上に述べた多次元分布における結晶構と呼ぶことが多いようです．

!%") #!%"*! !%C 言語や JAVA，Visual Basic，Excel や UNIX OS の標準（ただし )$*$!）（２）

的なライブラリーに入っている乱数ルーチンの中身は，現によって生成されます．ここで，! は各ビットごとに排他在でも混合合同法であることが多いようです．こうした標的論理和（XOR）を取ることを意味します．M 系列乱数に準化された乱数ルーチンでは，疑似乱数列の初期化のためついて詳しく説明するには，'!
!
"(である２つの要素に対に seed という値を最初にユーザが指定することになってするガロア体 GF（２）における議論が必要になりますの

いますが，これは乗積合同法や混合合同法における最初ので，詳しくは［１］を参照ください．ここではごく簡単な説明

.
#!!
$
!
"&は，発生する!を与えているに過ぎません．係数 %にとどめておきます．
（２）式に対して GF（２）上の多項式

乱数列の周期長がなるべく長く，乱数としての性質がよく.)".*""を特性多項式と呼びます．もし，.,!"という式なるように選ばれていますが，原理的にその周期長は "が ,##)!"の場合この多項式で割り切れ，,"#)!"に対を超えることはありません．例えば，C 言語の stdlib に実装

しては割り切れない場合，この特性多項式を特に )次の原されている RAND
（ ）関数は４バイト整数の疑似乱数ルー

始多項式と呼びます．このとき，（２）式から生成される数

チ ン で す が，そ の パ ラ メ ー タ は #!#""!$&"&#%&，

列の各ビットに着目すると，それぞれの '!
!
"(の並びのパ$#

ターンの周期が #)!"となることが知られています．原理$"$#"#$%&，"## であり，周期長は # です．
さて，古くから使われてきた線形合同法ですが，その乱的に，（２）式の漸化式が作る数列の周期は，#)!"より長数としての性質にはいくつかの問題があることが知られてくなることはないので，原始多項式から生成される数列はいます．「乱数らしさ」とは何かについては，また第４章

最大周期長を持つことになります．

で詳しく述べますが，線形合同法に関する問題点として

線形合同法と同様に Tausworthe 法でも，原始多項式とは，２進数で表現した乱数の下位ビットに，乱数列全体となるような %)!
*&の組み合わせを上手く選ぶ必要がありましての周期よりも短い周期性が現れること，「結晶構造」

す．代 表 的 な 例 と し て，%
，%
)
!
*&
!
!
'!(
#($&
"
#(*
%")&
#%

と呼ばれる分布の偏りが見られること，そして最近使われ等が知られています．さらに Tausworthe 法では )個のるようになった他の新しい疑似乱数ルーチンに比べて周期-ビットベクトル !%を初期条件として与える必要があり，長が短いことが挙げられます．
「結晶構造」とは，例えば漸化式を回すために )個分の過去の情報を保存しておく必! により%!".
!
!
"&区間%"! の整数疑似乱数列を +%#.
%#要があります．一つ留意しなければならないのは，原始多２９３Journal of Plasma and Fusion Research Vol.96, No.6 June 2020作を行ったベクトル 'を出力として取ります．

項式によって生成される数列が最大周期長を与えるのは.$"ビットの場合であり，.'#ビットになった場合に

，
&
%&-)
'(
"$%

は，%$ベクトルの初期条件を上手く設定しないと，. ビッ

，
&
&
)& #)
#$&
"'(
"%+

トベクトル %$全体のビットパターンの周期長が#(!"より

，
&
&
)& $)
$$&
#'(
#%,

短くなってしまう場合がある点です．幸い適切な初期ベク

（４）

，
'$&
&
)
$'(
$&%

トルの与え方は Tausworthe 法のサンプルプログラムに付随していますので，ユーザは自分で適当に初期ベクトルを

ここで -，+，,，%はある整数定数，(#!
$)は . ビットの定

定義せずにそれを利用しましょう．サンプルプログラムは数ベクトル，そして &-，%%はそれぞれ -ビット右シフ［３‐５］などで見つけることができます．

ト，%ビット左シフトを表します．
（３）式と（４）式は一見

Tausworthe 法の疑似乱数の性質については，多次元分

複雑な演算を行っているように見えますが，XOR，AND，

布における一様性が良好であることが特徴です．また，線ビットシフトなど簡素な２進数演算のみで構成されてお"%形合同法の周期長 "はプログラムが扱える最大の整数（#

り，コード化すると思いの外短く記述できます．

や #$#）を 超 え ら れ ま せ ん が，Tausworthe 法 の 場 合，MT 法の特徴の１つはその周期長の長さにあります．適(!
))の上手い組み合わせさえ見つけられれば，(&.と取(切にパラメータ &，'，*，"を選ぶことにより，生成されれるため格段に長い周期長 #(!"を実現できます．XOR

る数列 %$の周期長が #(!"$#'.!*!"になることが数学

演算もCPUで高速実行されますし，使用するメモリ量も線的に証明されています［７］．先に挙げた Tausworthe 法と形合同法より少し増える程度です．しかしながら，多次元

同様，その肝は漸化式（３）を決めるパラメータの選び方に分布の均等性以外の乱数としての性質については必ずしもあります．その詳細については数学の専門家ではない筆者

よいというわけではありませんでした［６］．また，（２）式の理解を超えるため，考案者の講義資料［３］を参照いただ

において各 .ビットが独立に計算されるため，ビット間でくことにします．ごく簡単に説明すると，数列の漸化式の情報のやり取りがないことが欠点として指摘されています．

形を表す特性多項式が /(!"の原始多項式になる，そのようなパラメータの組み合わせを見つけ出すことが一般的に!Mersenne-Twister 法は非常に難しいのですが，(が Mersenne 数と呼ばれる，

そこで，更なる改良として考案されたのが Mersenne-($#'!"の形を取る特別な素数の場合には，特性方程式

Twister 法［７（］以下 MT 法と略す）です．これは，日本のが原始多項式か否かの判定が容易になるという性質を使っ数学者，松本 眞氏によって考案された方法で，非常に長いているそうです．Mersenne-Twister の名前もこれが由来

周期性を持つ疑似乱数列が生成でき，乱数としての性質もになっています．また，（４）式のような操作は乱数の高次

良いことから，今では世界中の乱数を使うシミュレーショ

元分布の均等性をより高めるための工夫です．
「高次元分

ンで幅広く利用されるようになりました．

布の均等性」の説明は第４章で改めてしますが，例えば下に示す MT 法のサンプルでは周期長が #"''$&!"で，６２３次

MT 法の漸化式は形式的に以下のように記述されます．
.!* *%$"' $%$"& '(
%$
%$"")
!
,

元空間に均等分布する疑似乱数を生成します．

（３）

MT 法 に つ い て は，様 々 な 人 が Fortran，C，Java，

ここで '"&""です．また，%$は（２）式と同様に . ビッPython など様々な言語で実装しており，MT 法の考案者のトの２進数を表す横ベクトルであり，!は *!
!
"+を要素と松本氏のwebページにはそれらへのリンクがまとめてありする .#. の正方行列で，次のような形を取ります．
"
!
! " !
&
%
&
%! ! "&%&%&%! ! !
&
!$%
&
%
&
%
&
%
&
%
"
$
#
#! #" ## - #.!"

ます［８］．なお，初期のバージョンの MT 法のソースコー
ドには，乱数の初期化の seed の与え方があまり良くないものが使われており，紹介されているコードの中にも旧い初期化の仕方のままになっているものもあるので注意が必要
です．C言語については，MT法の作者自身が２００２年に公開した修正版（mt19937ar）［９］か，更にSIMD化して改良された SFMT 版（［９］中のリンク）をダウンロードして利用するのがよいでしょう．

ここで，"はある定数ベクトルで，演算 %!は以下のような

ここでは Fortran 版［１０］の利用法について紹介します．

演算に相当します．

なお，他の疑似乱数と同様，MT 法でも漸化式の形を決めるパラメータ("!
#!
$)を適切に選ぶ必要があ
&!
'!
*
!
-!
+
!
,
!
%
!

%!
!
%を左に１ビットシフト （%の最下位ビットが０）＝$"%を左に１ビットシフト ' "（%の最下位ビットが１）$#

りますが，これらは定数としてサンプルコード中に指定されているのでユーザが特に気にする必要はありません．
文献［１０］の mt19937ar.f は初期化と発生のルーチンのみ含

.!* *
また，%$
%$""は %$の上位 .!*ビットと %$""の下位,

まれています の で，これらを利用するためのメインプ*ビットをつなぎ合わせて１つの .ビットベクトルにする

ログラムを書きましょう．以下にサンプルプログラムという操作です．MT 法の漸化式は（３）式ですが，乱数とmtsample.f90を記載します（GitHubにも同じコードを公開して利用する際には，更に以下のような一連のかき混ぜ操してあります［１１］）．
２９４

Lecture Note

2. Random Number Generation

S. Satake and R. Kanno

input seed number (0 to continue)
program mtsample
implicit noneと聞かれるので，適当な自然数（４５６７など）を入れます．次

integer(kind=4), parameter :: N=624に，

integer(kind=4) :: seed,i,nrnd
input number to be generated
real(kind=8),allocatable :: rnd(:)

character(len=1) :: yesnoと聞かれるので，例えば１００と入れてみましょう．すると即!
!
"'区間疑似乱数が生成さ座に１００個の倍精度実数 の &real(kind=8) :: genrand_real1

れ，”mt_rand.txt”に格納されます．最後に，integer(kind=4) :: mti,initialized, mt(0:N-1)

common /mt_state1/ mti,initializedrecord the MT status? (y/n)
common /mt_state2/ mtと聞かれるので，y と入力すると乱数列を次回の計算に継!
print *,”input seed number (0 to continue)”

続するための情報が”mt_cont.dat”に記録されます．前回

read(5,*) seedの続きから疑似乱数を発生させたい場合は，seed に 0 を入if(seed>0) then

れて下さい．なお，mt_rand.txt と mt_cont.dat は実行する
度に上書きされるので注意して下さい．

call init_genrand(seed)

ここで，「疑似乱数を継続する」というのは，例えばあelsecall mt_initln

る seed からスタートして１万個を一気に生成したのと，同open(7,file=”mt_cont.dat”,action=”read”)

じ seed から始めて５千個の乱数を２回発生させたので同read(7,*) mti,initialized,mt

じものになるようにするということです．これは，実際にclose(7)疑似乱数を使って，１回のジョブに実行時間制限のあるス
パコンなどで，非常に長いシミュレーションを複数回繋いend ifで計算する場合に必要となります．サンプルの MT 法の場!
print *,”input number to be generated”

合は，項数６２４の漸化式になっているのでその状態を記録

read(5,*) nrnd

して，次に実行する時に読み直すことで継続できます．
０個ず
図１は seed=4567 で１００個一気に生成した場合と，５allocate(rnd(nrnd))

つを２回継続で生成した場合の乱数の分布をプロットした!

ものですが，両者が完全に一致していることが確認できま!
!
"#す．なお，mt19937ar.fにはgenrand̲real1の他にも，"，

do i=1,nrnd
rnd(i)=genrand_real1() ! [0,1]-real8!
!
"#など様々な区間の乱数を生成する function が定義さ&

end doれているので，用途に応じて rnd(i)=genrand_real1() の所!

を切り替えてください．

open(8,file=”mt_rand.txt”,action=”write”)

なお，サンプル中で適当に与えた seed の値の役割です

write(8,’(f18.15)’)rnd(:)

が，これは乱数の初期値を設定するものです．Tausworthe!
print *,”record the MT status? (y/n)”法や MT 法では漸化式をスタートさせるために，数百の漸read(5,*) yesno

化式の初期値 $!!%を指定する必要があります．MT 法の初if(yesno==”y”) then期化では１項だけの漸化式によるシンプルな疑似乱数ルー

open(7,file=”mt_cont.dat”,action=”write”)

チンを使ってその初期値を与えています（線形合同法より

write(7,*) mti,initialized,mt

は凝ったものですが）．seed の値を変えることは漸化式のclose(7)初期値のセットを変えることに相当し，より具体的に言えば周期長#"&&$%!"の数列のどこを開始点にするかを疑似乱

end if!

数で決めることになります．Tausworthe 法の場合，初期

end program mtsample

値を上手く設定しないと乱数列の周期長が期待通りの長さにならない場合があると述べましたが，MT 法ではビット

実 行 す る た め に は，［１０］か ら ダ ウ ン ロ ー ド し て き た

間で情報をやり取りするために，初期値の設定と実現されmt19937ar.f とこの mtsample.f90 を合わせてコンパイルしる周期長に依存性はありません［６］．

てください．例えば gfortran なら

ところで，このサンプルプログラムで生成する数を百万など大きめに取ると，それなりに計算時間が掛かることがgfortan mt19937ar.f mtsample.f90

体感できると思います．実はこのサンプルでは，１回につとすれば良いです．できあがった a.out を実行すると，

き１つの疑似乱数を返す function を必要な回数だけ call
する形になっています．実用上は，配列 rnd(nrnd)に一気に nrnd 個の疑似乱数を生成できた方がコードも書きやす
２９５

Journal of Plasma and Fusion Research Vol.96, No.6 June 2020

射線や，ダイオードの出力雑音，量子ゆらぎなどの物理現象を利用しようというアイデアは，古くからあり，日本では，遅くとも１
９５７年には物理乱数発生器を計算機に接続し，モンテカルロ計算に物理乱数が利用されていたようで
す［１２］．以下では，定電圧ダイオード（ツェナーダイオー
ド）の出力信号のゆらぎ，およびレーザー光の偏光の量子
ゆらぎを利用した乱数発生の原理について説明します．
!ダイオードの出力雑音による乱数発生
図１

MT 法で生成した乱数のサンプル．１０
０個を１度に生成した
場合と，５０個ずつ２回継続して生成した場合をプロットしてあります．

まず，定電圧ダイオードを利用した乱数発生器についてですが，出力信号に含まれるわずかな（数 μV 程度の）ゆらぎ（ノイズ波）をノイズ源として利用します［１３］．このゆいですし，効率もよさそうですよね？また，現在では大規らぎを ADC（Analog to Digital Converter）の変換範囲まで

模シミュレーションは MPI などの並列環境で実行するの増幅し，ある一定間隔の離散時刻でノイズ波のデジタル値が当たり前になっていますが，疑似乱数の並列化はどのよを取得します．各時刻に得られたデジタル値を００００〜１１１１

うにすればよいのでしょうか？こうした効率化，並列化にまでの１６個の２進数のうちの１つに対応させることで，各関するテクニックについては次章で解説したいと思いま

時刻における４ビットの乱数（０または１の値で表現されす．

る二値乱数）が生成されます．２回のデジタル値の取得で８ビット（つまり１バイト）の乱数となり，これを繰り返

２．
２ 物理乱数すことで，必要なバイト数の二値乱数を作ることができま

本節においては，乱数の生成にランダムな物理現象を利す．Fortran や C などのプログラミング言語における単精!
!
"!に含まれる一様乱数を生成す
度浮動小数点型の区間 "用した「物理乱数」について紹介します．乱数の生成に放!"""""""""""""""""""

乱 数 茶 話!"""""""""""""""""""!"""""""""""""""""""""""""
少し余談になりますが，疑似乱数は数値シミュレー

数調整」と称して一見無意味な行動をしてからイベント

ションの世界で利用されるだけではありません．一般のを起こすと（宝箱を開けるなど）決まった結果が起こる，

人々が疑似乱数と係わるのは，例えば，コンピュータというようなことをたまに聞きますが，これらは乱数のゲームで遊んでいる時ではないでしょうか．敵キャラク

代わりに使われている上に述べたようなデータの隠れた

タの出現位置や行動などを，プレイヤーが飽きないよう

規則性を見抜いて，都合のよい結果になるように「調整」にランダムにするなどの用途に疑似乱数が使われていましているわけです．

す．ところで筆者（S.S.）が子どもの頃のコンピュータと今回の記事の執筆のために調べていて興味深かったの

言えば，８ビットの CPU に 32 kB のメモリという，非常が，８ビット PC の時代から線形合同法や XOR を使ったに貧弱な計算資源の中でゲームを動かしていました．当

M 系列の乱数も利用されていたことを知ったこ とで然，乱数発生にメモリを割いたり複雑なアルゴリズムをす．ただし８ビットなので周期が非常に短いですが．線採用したりする余裕はありません．当時のパソコンのメ

形合同法を使ったゲームでも，偶奇が交互に出る性質をモリ（RAM）は放っておくと電気的に記録されたビット

理解せずにくじのイベントの当たり判定に使ってしまっ

情報が減衰してしまうので，時々データを書き直す必要たために，プレイヤーに癖を見抜かれてしまったというがありました．このタイミングをコントロールするため事例もあったようです．筆者のいとこは以前パチンコ機に CPU 内部にリフレッシュカウンタというレジスタが器のプログラマーをしていましたが，少ないプログラム

あり，タイマのようにCPUのクロックごとに値が増えて行数で統計性のよい疑似乱数を発生させる方法はないか

いきました．当時のCPUはせいぜい数kHzで動作していと相談されたことがありました．パチンコでは不正防止

ましたが，それでも人間から見れば高速にレジスタの値のために，当たりの確率が規定された通りになっているが変わるので，このレジスタの値が疑似乱数としてよくことが認可を受けるために必要な一方，プログラムに使使われていたようです．他にも，CPU には，命令を実行

えるメモリ量にも制約があるため疑似乱数ルーチンの選するたびにその結果に応じてレジスタの特定のビットが択で悩んでいたようです．

０になったり１になったりする，フラグレジスタという

このようにゲームやギャンブルの世界では疑似乱数はものがあります．これは本来，直前に行った計算がオー広く使われていますが，人工知能が発達すると，やがてバーフローしたとか，論理演算の真偽値を確認するためは疑似乱数に潜む隠れた規則性をあっさり見抜いて，囲のものですが，このフラグレジスタの値を乱数として使碁や将棋だけでなくコンピュータゲームでも人間のプレ

うこともありました．古いゲームの攻略法などで，「乱イヤーを凌駕してしまうかもしれませんね．

!"""""""""""""""""""""""""
２９６

Lecture Note

2. Random Number Generation

S. Satake and R. Kanno

る場合は，例えば，４バイトの乱数を用いて得た正整数をして，一様乱数表を利用する方法をここでは紹介しました

その最大値に１を足した数で割ることで得られます．以上

が，別の方法もあります．詳しくは，文献［１３］や関連するは，４ビ ッ ト の 乱 数 を 基準 に 説 明 し ま し た が，一 般 に

特許情報などを参照ください．

! ビットの乱数としても同様です．

このように乱数発生原理の概略がわかると，次に乱数発

ところで，上記の乱数発生法を実際の装置で実現化する
際には，得られる二値乱数の一様性が気になるのではない

生の速度が気になるのではないかと思います．この発生速度については，近年，640 MB/s（#5.1 Gbps）程度になってでしょうか．そこで，乱数の一様性を向上させる方法についるようです．複数台の乱数発生器を用意できるのであれいて紹介します．乱数生成のスピードを落とさずに一様性

ば，大規模並列計算機におけるシミュレーションでの利用を向上させる手法として，排他的論理和による補正手法がに（必要とする乱数の数や呼び出す頻度にも依りますが）考え出されました［１３］．文 献［１３］に 従 い，こ の 手 法 をほぼ問題の無いレベルの発生速度であると言えます．

２ビットの乱数を例として説明します．ADC が２ビット

定電圧ダイオードを用いた物理乱数発生器による乱数をの変換値として０〜３を出力するとし，変換値 "の発生確利用したい場合は，例えば，統計数理研究所の乱数ライブ

率を #" とします．ここで，""!
!
!
!
"
#
$で，"$
" ""
""!#

ラリー［１４］における乱数取得サービスから乱数をダウンです．あらかじめ与えた２ビットの一様乱数表から等しい

ロードする方法があります．この記事を書いている時点で確率"で補正値０〜３を生成し，ADCの変換値との排他的は，東京エレクトロンデバイス製の物理乱数発生器による論理和を取ります．表１に ADC の変換値と補正値との排他的論理和を行ったときの補正後の値（上段）とその生成乱数を取得できるようです．取得できる乱数のデータ型!
!
"$の単は，区間［−231，231］の３２ビット符号付整数，区間%確率（下段）を示しました．ここで，排他的論理和は，２．
１

精度浮動小数点数および倍精度浮動小数点数です．

節でも使いましたが，以下のような演算（! と表記）であ!レーザー光の偏光の量子ゆらぎによる乱数発生

り，表１では，２ビットの値００〜１１について，ビットごとに演算して補正後の値を得ています．

次に，光学的な装置を使った「量子乱数」と呼ばれる物
理乱数発生器の原理について紹介します．量子力学的な不! ! ! " !

確定性を持つ物理過程の観測値を乱数のソースとして利用! ! " " "しようというアイデアは１
９５０年代には既にありましたが" ! ! " "［１４］，初めに考えられたのは放射性物質から出てくる放射" ! " " !

線をガイガー計測器で計測した時の，２つのパルスの発生

例えば，変換値１と補正値３で考えると，!"! """ "!

間隔が Poisson 分布の確率分布をとるという性質を利用すですので，補正後の値は，２となります．表１を見ると，補るものでした．しかし放射性物質を線源として内蔵するこ

正後の値が２となるケースは全部で４つあり，その生成確との問題や，ガイガー計測器の計測速度に原理的な限界が率の和は，#!"!#""!##"!#$"""となるので，補正後のあったこと（数 kbps 程度），小型のレーザー発振器が実用

二値乱数１０は確率 "で生成されることがわかります．他の化されたことなどから，現在では光学デバイスによる量子補正後の値についても，すべて同じ確率 "で生成され，こ

乱数発生器が主流になってきています．文献［１５］の Tableのことから，生成された２ビットの乱数の一様性が向上しII にいくつかの異なる手法による光学的な量子乱数発生器ていることが分かります．上記の説明は，簡単のためにの一覧が載っていますが，乱数の生成速度は数 Mbps から

２ビットの場合で説明しましたが，２を ! に一般化して，

数 Gbps に達しています．

! ビットの乱数とした場合も同様な補正を行えば，生成すレーザー光を使った量子乱数の例としてここでは，［１６］る乱数の一様性が向上します．また，補正値の生成方法とのレーザーパルスの位相差を利用したものを紹介します．
図２のようなレーザー光の回路を組み，発振したレーザー

表１

パルスを一度２つの光ファイバーのケーブルに分岐させ，

ADC の変換値と補正値との排他的論理和により与えられる補正後の値（上段）とその発生確率（下段）．例えば，変換値１と補正値３の補正後の値は，２行‐４列目の上段に表示され，２です．この表は，文献［１３］の表１を参考に作成
しました．

０
０
#!"

補正値
１
２
１
２
#!"
#!"

３
３
#!"

１

１
#""

０
#""

３
#""

２
#""

２

２
##"

３
##"

０
##"

１
##"

３

３
#$"

２
#$"

１
#$"

０
#$"

０

変換値

片方の経路を長く取ることで，一方をちょうど１パルス分ずらします．その後，２つのパルスを再び合流させてから

図２

２９７

レーザー光を使った量子乱数発生器の概念図（［１６］の Fig.2を引用）．Adapted with permission from [16] ! The Optical Society.

Journal of Plasma and Fusion Research Vol.96, No.6 June 2020

その光の強度を計測します．レーザー光のパルスは量子力

NIST テストと呼ばれる一連の乱数検定テストもパスして学的揺らぎにより，一つ一つ異なる位相 #"を持った偏光にいます．このケースでの量子乱数の発生速度は理論的にはなっています．したがって２つのパルスを合成した後の40 Gbps を超えるとされていますが，実際には，Hash 関数

レーザー光の強度は，干渉によって ##!!#!!!$の依存$%#の演算処理による時間を含めると実用上のスピードは性を持ちます．#!と #!!!が量子力学的な揺らぎから決まる

数 Mbps〜Gbps 程度のようです［１９］．

ため，#!と #!!!は互いに独立で一様分布する確率変数であ今回紹介したレーザー光の偏光を利用した手法の他にると期待できるので，そこから観測光の強度の確率分布がも，レーザー光を２つにスプリットしたものを再度合わせ計算できます．仮に観測光の強度を #"%!
!&に規格化し!!

て計測することで，観測に混じる量子力学的な真空のゆらたとすると，その確率分布は !
"!
' !#"に比例すると期待さぎを増幅して観測するという原理に基づいた光学的量子乱れます．しかし，実際の観測値には観測のノイズやパルス数発生器も開発されています［２０］．この手法も信号をデジ幅に関係する補正項が入るので，単純に理論通りに分布す

タルサンプリングした生データの情報の一部を捨てる代わるわけではなく，図３のようなひずんだ凹型の分布になりりに統計的性質を向上させる手法を使って，計測データにます．この分布は，２つのパルスに何らかの相関関係があ含まれる量子的ゆらぎ以外の計測ノイズの影響を落とす後るため生じているわけではなく，上に述べたパルス間の位処理を使います．この原理に基づく乱数発生器がオースト

相差がランダムであることに起因しています．しかし，こラリア国立大学の研究グループによって公開運用されておの測定値には計測誤差やレーザー光源の揺らぎなどの誤差

り［２１］，乱数発生速度は 5.7 Gbps に達するそうです．インによる期待値からのずれも含まれており，これらの誤差はターネットを通じて生成された乱数を利用でき，様々なプ

必ずしもランダムであるとは言えません．また，非一様なログラム言語からこの乱数サーバのデータを取得するルー

強度分布のままでは乱数として使うには適さないので，デチンが公開されているので，関心のある方は［２１］のリンク

ジタル化した計測データをWhirlpool-Hash関数［１７］というから情報をたどってみてください．

ものを使って変換したものを乱数として使います．これところで比較参考のために MT 法の発生速度を見る

は，入力されたデータをある決まった規則に従って変換と，２．
１節のサンプルプログラムを Xeon のワークステー

し，入力データか ら は 予 測 で き な い 固 定 長 の デ ー タ列ション（クロック周波数 2.5 GHz）で実行した場合，３２ビッ（Hash 値）を返すもので，暗号化の技術の一つとして使わ

ト乱数を１千万個生成するのに約 0.2 秒，発生速度に直す

れているものです．この際，ある確率分布を持つ確率変数と約 1.6 Gbps でした．筆者（S.S）は数年前に［１６］の量子乱の最低エントロピー状態という理論を基に，デジタルサン

数発生器の試作機のテストに参加させてもらいましたが，プリングした元データのビット数よ り 低 い ビ ッ ト数の

当時はまだ乱数発生器からコンピュータのメインメモリにHash 値を乱数として採用することで，元のデータの情報

データを転送するところが律速となり，理論性能ほどの速量を一部落とす代わりに一様性の高い乱数列を得るという度は出ていませんでした．しかしそれでも当時のワークス手法［１８］を取っています．具体的には，１４ビットでサンプ

テーション用CPUでMT法を走らせた場合の１０分の１程度

8

リングした 1.2×10 個の観測データを Whirlpool-Hash 関数の速度が出ていました．物理乱数発生器は，まだあまり世で 1.25×108 個の７ビットのバイナリデータに変換していの中に出回っていませんが，［１４］や［２１］の乱数サーバの実

ます．こうして得られた数列は非常によい一様性を示し，績値を見ると，物理乱数の発生速度は疑似乱数に比較しう
る時代になっていると言え，今後の展開が期待されるとこ
ろです．

謝

辞

第２．
２節における定電圧ダイオードを利用した物理乱数の発生原理の詳細についてご教示いただいた，田村義保氏（情報・システム研究機構 統計数理研究所 特任教授・名誉 教 授）に 感 謝 い た し ま す．ま た，第２．
２節 に 対 し て，
河村学思氏（自然科学研究機構 核融合科学研究所 助教）より，貴重なコメントをいただきました．ここに感謝申し上
げます．

参考文献
図３

［１］津田孝夫：モンテカルロ法とシミュレーション（三訂版）（培風館，１９９５）．
［２］伏見正則，逆瀬川浩孝（監訳）：モンテカルロ法ハンド
ブック（朝倉書店，２
０１
４）．
［３］松本 眞：有限体の擬似乱数への応用，
http : / / www. math. sci. hiroshima-u. ac. jp / 〜 m-mat /短，長２つの経路を通ったそれぞれのパルス強度（〇と□の印，左の縦軸）と，２つを合わせて観測した場合の強度（＋印，右の縦軸）の分布．２つのレーザーパルスを合成す
ると干渉によって広がりを持った分布になります（［１６］のFig.3（a）を引用）．Adapted with permission from [16] !
The Optical Society.

２９８Lecture Note

2. Random Number Generation

S. Satake and R. Kanno

たけ

しん

すけ

佐 竹 真 介自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授，２００
３年総合研究大学院大学 博士（学術）．
モンテカルロ法を使った３次元磁場配位中の新古典輸送現象，新古典粘性のシミュレーションや最適化配位の研究が主なテーマ．乱数については深い思い入れがありますが，特にギャンブル好きというわけではありません．

かん

の

りゅう

た

ろう

菅 野 龍 太 郎
自然科学研究機構 核融合科学研究所 ヘリ
カル研究部 核融合理論シミュレーション
研究系 准教授．乱数を用いた確率論的な計算手法であるモンテカルロ法全般に興味があります．最近の研究では，プラズマの衝突輸送現象に対するモンテカルロ法に基づいたドリフト運動論シミュレー
ションを行っています．

!!!!!!!!!!!!!!!!!!!!!!!!!!

さ!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!

［１２］統計数理研究所 統計科学技術センター 計算機の歴史
https://www.ism.ac.jp/computer̲system/jpn/outline/
hist.html［１３］田村義保 他：日本統計学会誌 35, 201 (2006).
［１４］統計数理研究所 乱数ライブラリー
http://random.ism.ac.jp/［１５］Miguel Herrero-Collantes et al., Rev. Mod. Phys. 89,
015004 (2017).
［１６］C. Abellán et al., Optics Express 22, 1645 (2014).
［１７］P.S.L.M. Barreto and V. Rijmen, The Whirlpool hashing
function, https://www.researchgate.net/publication/
228610491［１８］N. Nisan and A. Ta-Shma, J. Comput. System Sci. 58, 148(1999).
［１９］Carlos Abellán et al., Optica 3, 989 (2016).
［２０］T. Symul et al., Appl. Phys. Lett. 98, 231103 (2011).
［２１］The ANU Quantum Random Numbers Server,
https://qrng.anu.edu.au/!!!!!!!!!!!!!!!!

TOKYOHOMEPAGE/TEACH/1011.pdf［４］統計数理研究所 乱数ライブラリー 乱数について
http://random.ism.ac.jp/info01.html［５］宮村 修，牧野 純：ベクトル化乱数発生プログラム，大阪大学大型計算機センターニュース 75, 39 (1989).
http://hdl.handle.net/11094/65855［６］松本 眞，栗田良春：Twisted GFSR：新しい乱数発生法，京都大学数理解析研究所講究録 85, 86 (1993).
http://www.kurims.kyoto-u.ac.jp/〜kyodo/kokyuroku/
contents/pdf/0850-08.pdf［７］M. Matsumoto and T. Nishimura, ACM T MODEL COMPUT S 8, 3 (1998).
［８］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
VERSIONS/versions.html［９］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
MT2002/mt19937ar.html［１０］http://www.math.sci.hiroshima-u.ac.jp/ 〜 m-mat/MT/
VERSIONS/FORTRAN/mt19937ar.f［１１］https://github.com/satakeshinsuke/random̲number!!!!!!!!!!!!!!!!!!!!!!!!!!

２９９