TeX で乱数を使う

この記事は TeX & LaTeX Advent Calendar 2018 の19日目の記事です。18日目は keisuke495500 さんでした。 20日目は isaribi_saitoh さんです。

本記事では，TeX で乱数を取り扱う方法，および乱数の活用法を解説します。

線形合同法による擬似乱数生成

決定的な動きをする計算機で真の乱数を作り出すことは原理的に不可能なので，（熱雑音を利用するハードウェア乱数生成器など外部からの不規則な入力を用いてランダムネスを実現する方法はありますが）通常は見た目上ランダムに見える数列を計算式に従って規則的に作り出した擬似乱数を用います。擬似乱数を作り出す計算式としては，素朴な方法である線形合同法 (Linear Congruential Generator, LCG) が，実装の単純さおよび動作の高速さゆえに，昔から広く用いられています。

線形合同法では，漸化式 $$ X_{n+1} = (aX_n +c) \bmod m $$ によって0以上m未満の擬似乱数列 {$X_n$} を作り出します。このときパラメータ $a$，$c$，$m$ の選び方によって生成される擬似乱数の“質”が大きく左右されます。よく使われるパラメータのリストが Wikipedia (en) にあります。下手なパラメータを用いると偶奇が交互に繰り返す実装になる場合があります。かつて，Xbox 360 でこのトラップを踏んでしまったと思われる不具合によりゲームソフトの回収騒ぎが起こったことがありました。また，弾幕シューティングゲームなどでも，乱数の質が悪いと“安全地帯”が生まれる原因となりますので，質の高い乱数生成は，カジュアルな用途においても意外に重要です。（なお，ゲームと乱数の歴史についてはこの記事が参考になります。）

seed の与え方

線形合同法のように漸化式を用いて疑似乱数生成を行う場合，seed すなわち初期値 $X_0$ をいかに与えるかが問題となります。特定の数値に固定していると毎回同じ乱数列になってしまいます（出力を固定・再現したい場合はむしろその性質が重要になります）。「毎回異なるseed」を与える方法としては，「現在時刻」や「システムが起動してからの時間」など，時間に関連する値がよく用いられます。

TeX による時刻取得

TeX による時刻取得の方法は mod_poppoさんの記事 に詳しい解説があります。この記事中にありますように，元々の TeX には，時刻取得用プリミティブは次の4つしかありません。

• \year
• \month
• \day
• \time：その日の午前0時から経過した分（時*60+分）

つまり元々の TeX では「毎秒変わる値」を取得することができません。それゆえ，元々の TeX では「毎秒ごとにコンパイル結果が変わる」文書を作成することができないわけです。時刻に関する値を用いて seed を自動設定するタイプの擬似乱数生成機能を用いた場合，次の異なるコンパイル結果を得るまでに最大1分待つ必要があります。

一方，近年の拡張エンジンでは，PDF の「作成日時」に入れるための「現在時刻を秒単位で取得する拡張プリミティブ」が次のように用意されています。

• pdfTeX：\pdfcreationdate
• e-pTeX, e-upTeX：\pdfcreationdate（TeX Live 2014 以降）
• LuaTeX：\pdffeedback creationdate
• XeTeX：\creationdate（TeX Live 2019 以降）

これらは，例えば 2018年12月19日08時45分59秒（日本時間）なら

D:20181219084559+09'00'

のような出力を行いますので，この出力文字列をパースすることで「毎秒異なる値」を取得できます。

TeX で擬似乱数を生成する種々の方法

(La)TeX 上で擬似乱数を利用する方法として，次のように様々な手が考えられます。

• lcg パッケージ
• fp パッケージの \FPrandom
• pgf パッケージの rnd / rand / random
• 専用プリミティブ
• expl3 の \fp_eval:n
• LuaTeX の math.random
• LuaTeX でメルセンヌ・ツイスターの Lua 実装ライブラリを呼び出す

このそれぞれを見てゆきましょう。

lcg パッケージ

lcg パッケージはかなり昔からある LaTeX パッケージです。名前の通り LCG（線形合同法）で擬似乱数を生成します。パラメータが $a=7^5=16807$，$m=2^{31}-1=2147483647$，$c=0$ の線形合同法となっています。これは Park–Miller random number generator や MINSTD (minimal standard) と呼ばれる標準的な疑似乱数生成器です。

seed は手動で指定できるほか，自動にした場合は現在の日時・ファイル上の行番号・ページ番号から自動設定されます。

【コード例】

\documentclass{article}
\usepackage{lcg}
\newcount\cnt
\begin{document}

%% seedを固定して指定された範囲の乱数を得る
\reinitrand[counter=random,first=10,last=20,seed=1234]

\cnt=0
\loop
  \rand %%% 次の乱数へ進める
  \arabic{random},
  \advance\cnt by 1
\ifnum\cnt<10\repeat

% seedは現在時刻から自動取得
\reinitrand[counter=random2,first=300,last=500]

\cnt=0
\loop
  \rand %%% 次の乱数へ進める
  \arabic{random2},
  \advance\cnt by 1
\ifnum\cnt<10\repeat

\end{document}

【出力】

毎秒異なる seed を指定したい場合

seed を自動にした場合は現在の日時を使用した seed が設定されるため，先に述べたように，このままでは次の異なるコンパイル結果を得るまでに最大1分待つ必要があります。そこで，現在の秒数を使用した seed を設定してみましょう。例えば pdfTeX / e-pTeX / e-upTeX エンジンが搭載している \pdfcreationdate を使う場合，次のような感じで秒数を取得・使用できます。

\documentclass{article}
\usepackage{lcg}
\newcount\cnt

\makeatletter
\def\seed{\the\numexpr\expandafter\@NOW\pdfcreationdate\relax}
\catcode`D=12
\def\@NOW D:#1#2#3#4#5#6#7#8{#1#2#3#4#5#6#7#8+\@@NOW}
\catcode`D=11
\def\@@NOW#1#2#3#4#5#6+#7\relax{#1#2#3#4#5#6\relax}
\makeatother

\begin{document}
current seed = \seed
% 2018年12月19日08時45分59秒なら 20181219+084559 を計算して 20265778 という文字列を得る

% 秒数を含む現在時刻に基づくseedを指定
\reinitrand[counter=random,first=10,last=20,seed=\seed]

\cnt=0
\loop
  \rand %%% 次の乱数へ進める
  \arabic{random},
  \advance\cnt by 1
\ifnum\cnt<10\repeat
\end{document}

fp パッケージの \FPrandom

次に，これも昔からある fp パッケージです。fp-random.sty で定義された rand 関数を使用することで，[0,1) 上の一様乱数が得られます。乱数生成はやはり $a=7^5=16807$，$m=2^{31}-1=2147483647$，$c=0$ の線形合同法(MINSTD)で，得られた 0, ..., $m-1$ の値を $m$ で割ることで [0,1) の実数値としています。seedは \FPseed というカウンタで指定します。それが 0（デフォルト）の場合は 123456789 という固定値が採用されます。

【コード例】

\documentclass{article}
\newcount\cnt
\usepackage[nomessages]{fp}

\begin{document}
\FPseed=987654 %%% seedを固定値に指定
\loop
  % 新しい乱数値を \rnd に格納
  \FPrandom{\rnd}
  \rnd,
  \advance\cnt by 1
\ifnum\cnt<10\repeat
\end{document}

【出力】

毎秒異なる seed を指定したい場合は，先述の \pdfcreationdate などを使う方法でOKです。

pgf パッケージの rnd / rand / random

pgf パッケージが提供する pgfmathfunctions.random.code.tex の中には，$a=69621$，$m=2^{31}-1=2147483647$，$c=0$ の線形合同法による擬似乱数生成器が実装されています。

• rnd / \pgfmathrnd：0 ～ 1 の実数値一様乱数
• rand/ \pgfmathrand：-1 ～ 1 の実数値一様乱数
• random(x,y) / \pgfmathrandom{x,y}：x ～ y（両端含む）の整数値一様乱数

が利用できます。seed は \pgfmathsetseed で指定できます。デフォルトでは \time および \year をもとにした値が設定されます。 したがって異なる結果を得るには最大1分待つ必要がありますので，毎秒異なる seed を指定したい場合は，先述の \pdfcreationdate などを使う必要があります。

【コード例】

\documentclass{article}
\usepackage{pgffor}
\begin{document}
\pgfmathsetseed{18}%%% seedを固定値に指定
% 10～20の整数値乱数を10個生成
\foreach \x in {1,...,10}{\pgfmathparse{random(10,20)}\pgfmathresult,} 
\end{document}

【出力】

専用プリミティブ

最近の TeX エンジンには，乱数生成のための専用のプリミティブが内蔵されているものが多いです。一様乱数から正規乱数へはボックス＝ミュラー変換で変換できますが，プリミティブの側で変換してくれるのはありがたいです。また，プリミティブであるからこそ完全展開可能であるという点も嬉しい性質です。

pdfTeX / e-pTeX / e-upTeX の場合は

• \pdfuniformdeviate<number>：0以上 number 未満の整数値一様乱数を生成
• \pdfnormaldeviate：平均値 0，標準偏差 65536 の正規分布に従う整数値乱数を生成

LuaTeX / XeTeX の場合は

• \uniformdeviate‌<number>：0以上 number 未満の整数値一様乱数を生成
• \normaldeviate‌：平均値 0，標準偏差 65536 の正規分布に従う整数値乱数を生成

となります（pTeX 系は TeX Live 2017 以降，XeTeX は TeX Live 2019 以降）。

これらのエンジンには現在の時刻を正確に取得する機能が備わっているため，seed はデフォルトで「現在時刻によって異なる値」が指定されます。正確には，現在時刻をマイクロ秒単位で取得して seed に利用されるため，1秒以内に再コンパイルしてもほぼ確実に異なる乱数列が得られます。

seed を手動指定して出力を固定したい場合は \pdfsetrandomseed<number> で指定できます。現在の seed は \pdfrandomseed で取得できます。（LuaTeX / XeTeX の場合はそれぞれ \setrandomseed, \randomseed となります。）

【コード例】

\documentclass{article}
\newcount\cnt
\begin{document}
%% seed を手動指定で固定。自動の場合はシステム時刻から毎回異なる値になる。
\pdfsetrandomseed123456
\loop
  % 0以上100未満の一様乱数を生成して50を足すことで50以上150未満の一様乱数を得る（完全展開可能！）
  \edef\rnd{\the\numexpr\pdfuniformdeviate100+50\relax}
  \rnd,
  \advance\cnt by 1
\ifnum\cnt<10\repeat
\end{document}

【出力】

実装

e-(u)pTeX でのこれらのプリミティブの実装はこちらで見られます。これは pdfTeX の実装から来ており，さらにそれは Knuth による METAFONT の実装 から来ています。

この実装は，線形合同法ではなく，Knuth が The Art of Computer Programming に記載した “portable random number generator” によるものです。The Art of Computer Programming, Vol.2, 2nd ed., §3.6 には次の FORTRAN コードが載っています。

   FUNCTION IRN55(IA)
   DIMENSION IA(1)
C  ASSUMING THAT IA(1), ..., IA(55) HAVE BEEN SET UP PROPERLY,
C  THIS SUBROUTINE RESETS THE IA ARRAY TO THE NEXT 55 NUMBERS
C  OF A PSEUDO-RANDOM SEQUENCE, AND RETURNS THE VALUE 1.
   DO 1 I = 1, 24
     J = IA(I) - IA(I+31)
     IF (J .LT. 0)  J = J + 1000000000
     IA(I) = J
1  CONTINUE
   DO 2 I = 25, 55
     J = IA(I) - IA(I-24)
     IF (J .LT. 0)  J = J + 1000000000
     IA(I) = J
2  CONTINUE
   IRN55=1
   RETURN
   END

※このコードは現在の 3rd edition では削除されているようです。実際，Preface to the Third Edition のところに most significant changes の一つとして “§3.6 (about portable random-number generators)” が挙げられていました。

既知の不具合

2016年に報告されていますように，このプリミティブが作り出す疑似乱数は，ときに強い規則性を持ってしまいます。具体的には，

\loop\ifodd\pdfuniformdeviate 536870912\relax\else\repeat

が無限ループとなります。つまり，\pdfuniformdeviate 536870912 が毎回必ず偶数を返すようになっています。この不具合は未だ解消されていません。

expl3 の \fp_eval:n

expl3 の \fp_eval:n が解釈できる関数として，次が用意されています。

• rand()：[0,1) 上の実数値一様乱数を得る
• randint(x)：1以上x以下の整数値一様乱数を得る
• randint(x,y)：x以上y以下の整数値一様乱数を得る

これらは，上記のエンジン内蔵乱数生成プリミティブを内部で利用しています（そのため TeX Live 2018 では，そのプリミティブを持たない XeTeX では使えません）。また，seed は \sys_gset_rand_seed:n で指定できます。

【コード例】

\documentclass{article}
\usepackage{xparse}
\begin{document}
\ExplSyntaxOn
  % 1以上99以下の整数値乱数を10個生成
  \prg_replicate:nn{10}{\fp_eval:n{randint(1,99)},}
\ExplSyntaxOff
\end{document}

【出力】

LuaTeX の math.random

LuaTeX の場合は，Lua の標準ライブラリに含まれる math.random で乱数を生成できます。

• math.random()：[0,1) 上の実数値一様乱数を得る
• math.random(x)：1以上x以下の整数値一様乱数を得る
• math.random(x,y)：x以上y以下の整数値一様乱数を得る

seed は math.randomseed(1234) のようにして指定できます。

内部実装としては，Luaインタプリタを実装するC言語の乱数生成関数が呼ばれています。

【コード例】

\documentclass{article}
\begin{document}
% 30以上40以下の整数値乱数を10個生成
\directlua{for i = 1, 10 do
  tex.print(math.random(30, 40) .. ", ")
end
}
\end{document}

【出力】

LuaTeX でメルセンヌ・ツイスターの Lua 実装ライブラリを呼び出す

せっかくの LuaTeX ですから，線形合同法ではなく，もっと質の良い乱数を生成するメルセンヌ・ツイスターを使ってみましょう。メルセンヌ・ツイスターは，有限体 $\mathbb{F}_2$ とメルセンヌ素数の性質を高度に利用した高品質な疑似乱数生成アルゴリズムです。 Caldas Lopes 氏による randomlua ライブラリ には，Pure Lua によるメルセンヌ・ツイスターの実装が含まれています。これをダウンロードして

require('randomlua')

して使えば，次のようにしてメルセンヌ・ツイスターが使えます。

m = twister(seed)

とすれば，メルセンヌ・ツイスターの generator object を生成します。引数 seed を与えなければ現在時刻をもとに自動的に設定されます。

この generator に対して，次のメソッドが使えます。

• :random(0)：メルセンヌ・ツイスターが生成した加工前の 31bit 整数値乱数を得る
• :random()：[0,1) 上の実数値一様乱数を得る
• :random(x)：1以上x以下の整数値一様乱数を得る
• :random(x,y)：x以上y以下の整数値一様乱数を得る

【コード例】

メルセンヌ・ツイスターが生成した 31bit 整数値乱数を10個並べます。

\documentclass{article}
\begin{document}
\directlua{
require('randomlua')
m = twister()
for n = 1, 10 do
  tex.print(m:random(0) .. ", ")
end
}
\end{document}

【出力】

乱数を使う：モンテカルロ・シミュレーション

さて，TeX における様々な疑似乱数生成法が分かったところで，それを使って何かをやってみましょう。

まずは，乱数の活用法として定番のモンテカルロ・シミュレーションで，これまたド定番の円周率$\pi$の近似値計算をやってみましょう。円周率の近似値計算法としては収束速度が最悪な部類に属しますが，まあ練習です。可視化した方が面白いので，TikZ で可視化しましょう。それなら pgf にビルトインされた乱数生成器を使うのが素直でしょうかね。

[0,1) の一様乱数を2つ生成して (x,y) としてプロットする試行を繰り返せば，$x^2+y^2<1$ を満たした点の個数の割合が，四分円の面積 $\displaystyle\frac{\pi}{4}$ に近づくはずです。それを4倍することで $\pi$ の近似値を得ようというわけです。

【コード例】

\documentclass[pdftex]{article}%%% 適切なドライバ指定をする
\usepackage{scsnowman}

\newcount\cnt
\newcount\hitcount
\newcount\samplecount
\samplecount=3000 %% プロットするサンプル数

\newsavebox\redsnowman
\newsavebox\bluesnowman
% ヒットした点を表す赤☃
\setbox\redsnowman\hbox{\scsnowman[hat=red,buttons=red,muffler=red,scale=0.8]}
% ヒットしなかった点を表す青☃
\setbox\bluesnowman\hbox{\scsnowman[hat=blue,buttons=blue,muffler=blue,scale=0.8]}

\begin{document}
\centering
\begin{tikzpicture}[x=1pt,y=1pt,scale=200,>=stealth]
\coordinate (O) at (0,0);
\node[anchor=north east] at (O) {0};

\draw[->] (O) -- (1.1,0) node[anchor=west] {$x$};
\draw[->] (O) -- (0,1.1) node[anchor=south] {$y$};

\coordinate (A) at (1,0);
\draw (A) node [anchor=north] {1} arc(0:90:1) node [anchor=east] {1} -| (A);

\loop
  \pgfpointdiff{\pgfpoint{0pt}{0pt}}{\pgfpoint{rnd}{rnd}}
  \pgfgetlastxy{\xdim}{\ydim}
  \pgfmathsetmacro{\dist}{veclen(\xdim,\ydim)}
  \ifdim\dist pt<1pt\relax
    \node at (\xdim,\ydim) {\usebox\redsnowman};
    \advance\hitcount by 1
  \else
    \node at (\xdim,\ydim) {\usebox\bluesnowman};
  \fi
  \advance\cnt by 1
  \ifnum\cnt<\samplecount\repeat
  \global\hitcount=\hitcount
\end{tikzpicture}

\makeatletter
$\pi\simeq\strip@pt\dimexpr\the\hitcount pt*4/\the\samplecount\relax$
\makeatother
\end{document}

【出力】

次に，$\displaystyle\int_0^1\frac{dx}{1+x^2}=\frac{\pi}{4}$ であることを利用して，同様のモンテカルロ・シミュレーションを実行してみましょう。

【コード例】

\documentclass[pdftex]{article}%%% 適切なドライバ指定をする
\usepackage{scsnowman}

\newcount\cnt
\newcount\hitcount
\newcount\samplecount
\samplecount=3000 %% プロットするサンプル数

\newsavebox\redsnowman
\newsavebox\bluesnowman
% ヒットした点を表す赤☃
\setbox\redsnowman\hbox{\scsnowman[hat=red,buttons=red,muffler=red,scale=0.8]}
% ヒットしなかった点を表す青☃
\setbox\bluesnowman\hbox{\scsnowman[hat=blue,buttons=blue,muffler=blue,scale=0.8]}

\begin{document}
\centering
\begin{tikzpicture}[x=1pt,y=1pt,scale=200,>=stealth]
\coordinate (O) at (0,0);
\node[anchor=north east] at (O) {0};
\draw[->] (O) -- (1.1,0) node[anchor=west] {$x$};
\draw[->] (O) -- (0,1.1) node[anchor=south] {$y$};
\node[anchor=west] at (1, 0.5) {$\displaystyle y=\frac{1}{1+x^2}$};

\coordinate (A) at (1,0);
\coordinate (B) at (0,1);
\draw (A) node [anchor=north] {1} |- (B) node [anchor=east] {1};
\draw[thick, smooth,samples=100,domain=0:1] plot (\x, {1/(1+\x*\x)});

\loop
\pgfpointdiff{\pgfpoint{0pt}{0pt}}{\pgfpoint{rnd}{rnd}}
\pgfgetlastxy{\xdim}{\ydim}
\pgfmathparse{1/(1+\xdim*\xdim)}

\ifdim\ydim <\pgfmathresult pt\relax
  \node at (\xdim,\ydim) {\usebox\redsnowman};
  \advance\hitcount by 1
\else
  \node at (\xdim,\ydim) {\usebox\bluesnowman};
\fi
\advance\cnt by 1
\ifnum\cnt<\samplecount\repeat
\global\hitcount=\hitcount
\end{tikzpicture}

\makeatletter
$\pi\simeq\strip@pt\dimexpr\the\hitcount pt*4/\the\samplecount\relax$
\makeatother
\end{document}

【出力】

TeX で $\pi$ の近似値計算をやっても大した精度は期待できませんので，次回は乱数のもっと実用的な用途を考えてみましょう。