.NET System.Random の実装と欠陥について～重箱の隅をつつきたおす～

はじめに

.NET (C#) には、組み込みの擬似乱数生成器 System.Random が用意されています。ここでは、 System.Random の実装と性質・ひいては欠陥について、可能な限り深くまで調べて難癖をつけていきます。

結構いろいろあって内容が増えてしまったので、雰囲気をつかみたい方は目次をみてください。

はじめに
内部実装
外部実装・その他の問題
おわりに

内部実装

ここでは、 System.Random に実装されているメソッドの中身を詳しく見ていきます。

内部実装の参照

System.Random の内部実装は以下サイトで確認することができます。

両者は基本的には同一ですが、 new Random() の動作だけ異なります。詳細については new Random() の項を参照してください。

System.Random の定義だけ並べると以下のようになっています。ひとつひとつ追っていきたいと思います。

public class Random
{
    public Random();
    public Random(int Seed);

    public virtual int Next();
    public virtual int Next(int maxValue);
    public virtual int Next(int minValue, int maxValue);
    public virtual void NextBytes(byte[] buffer);
    public virtual void NextBytes(Span<byte> buffer);   // .NET Standard 2.1~
    public virtual double NextDouble();
}

`new Random()`

ライブラリ側で int 型のシード値を選択し、後述する new Random(int Seed) を呼び出します。このコンストラクタだけは、.NET Framework と .NET Core 以降で動作が異なります。

.NET Framework: 同タイミングでの初期化でシードが重複する問題

.NET Framework では、 Environment.TickCount (システム起動後のミリ秒単位の経過時間) をシードに用います。そのため、同じタイミングで生成された Random はかなり高い確率で同じ乱数列を返すようになってしまいます。例えば、同一フレーム中に別々の場所でインスタンス生成した場合などで、意図せず出力が重複する事故が予想されます。余談ですが、手元の C# Interactive で試してみたところ、1 ミリ秒の間に 2 万個程度の同じ系列を返すインスタンスを生成できました。

.NET Core: シード重複問題の改善

対して .NET Core 以降では、このような問題が起こらないよう対策がとられています。具体的には、スレッドごとに独立な static Random インスタンスの Next() をシードとします。そして、シードを生成するほうは普通の static な Random インスタンスの Next() をシードとしています。そしてそして、そのシードを生成するほうは、extern 関数 SystemNative_GetNonCryptographicallySecureRandomBytes をシードとしています。
要するに、外部取得シード → static Random → ThreadStatic Random → インスタンス Random 、という流れで初期化されていきます。

なお、 SystemNative_GetNonCryptographicallySecureRandomBytes の Unix 環境における実装は、

Arc4Random が利用できればそれを、
できなければ、/dev/urandom/ と lrand48() を xor したものを

返すようです。

ところで、 Next() の性質上、シードに int.MaxValue は絶対に選択されなくなりました。 (.NET Framework では起動から 25 日弱経てば生成される可能性がありました。) 詳細は Next() の項で述べます。

余談: 同じアルゴリズムによる初期化について

乱数のシード設定に同じ擬似乱数生成アルゴリズムを用いるのは果たして大丈夫なのか、という懸念が少しだけあります。シード設定には根本的に構造が異なる擬似乱数生成器を用いるほうがよい、気がするので… xoshiro/xoroshiro の作者のサイトにはこのような記述があります。

as research has shown that initialization must be performed with a generator radically different in nature from the one initialized to avoid correlation on similar seeds.

(肝心の "research" の中身にアクセスできないので厳しいですが、 Mersenne Twister の作者が書かれているようです)

`new Random(int Seed)`

Seed をもとに内部状態を初期化します。同じ Seed を指定すれば、完全に同じ出力列が得られます。

初期化に使われているアルゴリズムは、微妙に根拠が怪しい気もするラグ可変の Fibonacci 法らしきものと、遷移関数と異なるラグを持つ Lagged Fibonacci 法 4 周分です。 (もし根拠をご存知の方がいらっしゃれば教えていただけるとありがたいです。)

なるべく分かりやすいように整理した等価なコードを示します。

// 内部状態:
int[] _seedArray = new int[56];      // [0] は未使用
int _inext;   // I(ndex) Next ?
int _inextp;  // I(ndex) Next P ?


// Seed の絶対値を取る
Seed = (Seed == int.MinValue) ? int.MaxValue : Math.Abs(Seed);

// 定数 161803398 は黄金比が由来とのこと
_seedArray[55] = 161803398 - Seed;


// ラグ可変の Fibonacci 法らしきもの (？) 
{
    int prevState = _seedArray[55];
    int nextState = 1;

    for (int i = 1; i < 55; i++)
    {
        int leapIndex = 21 * i % 55;
        _seedArray[leapIndex] = nextState;

        // (prevState - nextState) mod int.MaxValue のイメージ
        nextState = prevState - nextState;
        if (nextState < 0)
            nextState += int.MaxValue;

        prevState = _seedArray[leapIndex];
    }
}

// ラグの異なる Lagged Fibonacci 法で 4 周する
for (int k = 0; k < 4; k++)
{
    // lagIndex は (1 + (index + 30) % 55) に等しい  
    for (int index = 1, lagIndex = 32; index < 56; index++, lagIndex++)
    {
        // ([index] - [lagIndex]) mod int.MaxValue のイメージ
        _seedArray[index] -= _seedArray[lagIndex];
        if (_seedArray[index] < 0)
            _seedArray[index] += int.MaxValue;

        if (lagIndex == 55)
            lagIndex = 0;
    }
}

// ラグ (インデックス) の初期化
_inext = 0;
_inextp = 21;        // 💀💀💀

なお、初期化後の _seedArray の各要素は $0 \le x \le 2^ {31} - 2$ の範囲になります。 int.MaxValue は含まない非負の値をとるのが重要なポイントです。

絶対値が同じ `Seed` は同じ乱数列を生成する

Seed は絶対値をとってから計算するため、例えば new Random(401) と new Random(-401) は同じ乱数列を生成します。なお、Seed == int.MinValue の場合はなぜか int.MaxValue として扱われます。したがって、全ての int 値の入力 ( $2^ {32}$ パターン) に対して、実質のシード値は

0 : 0 の 1 個
1 <= x < int.MaxValue : x, -x の 2 個
int.MaxValue : int.MaxValue, -int.MaxValue, int.MinValue の 3 個

と微妙に偏ります (？) 。

初期状態のパターン数が少ない

以上から、初期状態は $2^ {31}$ パターンに限定されます。内部状態のとりうるパターン数 $(2^ {31} - 1) ^ {55} \approx 2^ {1075}$ *1 よりかなり小さく、ブルートフォースでもなんとかなる程度の量です。

内部状態 `_seedArray` の無駄

_seedArray は要素が 56 個あるにもかかわらず、実際に使っているのは [1] ~ [55] までで [0] は使用されておらず、残念感があります。55 個しか使わないのが理論上正しいので、 0 基点にして 55 要素で良かったのでは、と思ってしまいます。

余談: 初期化時の制約について

System.Random では上記のアルゴリズムで初期化が行われていますが、別のアルゴリズムで差し替えようと思ったとき (?) 、問題なく初期化するにはどうすればよいでしょうか。擬似乱数生成器の設計によりますが、一般に内部状態を初期化する際は特殊な条件を満たす必要があります。例えば Xorshift では全 0 以外で初期化する必要があり、そうでないと永遠に 0 を吐き出し続けてしまいます。

System.Random では、内部状態にひとつでも非 0 の要素があれば (全 0 でなければ) 大丈夫です。ただしこれは法 $m = 2 ^ {31} - 1$ 上で正しく設計されていた場合の話です。正しく、がどういう意味かは Next() の項で説明しています。

もし法 $m = 2 ^ n$ ならば、ひとつでも奇数の要素があれば大丈夫です。 (つまり、どれかは最下位ビットが 1 になっている必要があります。) 少しだけ条件が厳しくなります。

ただし、これらは最低条件 (最長周期を満たす条件) であって、まともな初動の出力が得られるかはまた別の話です。例えば {1, 0, 0, 0, ..., 0} といった極端な初期状態を設定すれば、しばらくの間 0 とか 1 のようなものが出力され続けるでしょう。現実には、なるべくランダムなバイト列で埋めてから、これらの条件を満たすかどうかチェックする (ダメならやり直す) 、という流れになります。ただ、 System.Random ではその手のチェックは入っていません。要素が 55 個もあるのでどれかは条件を満たすだろう、という楽観的な見込みなのかもしれません (実際問題ないとは思いますが) 。

また _inext や _inextp といったインデックスについては、別種の注意を払って初期化する必要があります。コメントの意味深な 💀💀💀 については、次の Next() の項で解説します。

`int Next()`

0 <= x < int.MaxValue の範囲の擬似乱数を出力します。

Next() は他の出力系メソッドのベースになっています。他のメソッドは、このメソッドの出力を加工して出力しています。 (実際のコード上では InternalSample() というメソッドに実装されていますが、Next() と完全に等価なので、簡単のためこのように表現します。)

Next() 、つまり System.Random の遷移関数は、 MSDN にもあるように Lagged Fibonacci 法の一種である「引き算法」で実装されています。

The current implementation of the Random class is based on a modified version of Donald E. Knuth's subtractive random number generator algorithm.

ちなみに、文献としては TAOCP が参照されていますが、実際のコードは Numerical Recipes の ran3 が原典のようです。

引き算法の遷移関数を漸化式で書くと以下のようになり、名前通り引き算で遷移していることが分かります。

$\begin{aligned} x_{n} = x_{n-j} - x_{n-k} \pmod{m} \end{aligned}$

なお、ここでの $\mathrm{mod}$ は剰余演算子 % とは微妙に異なり、必ず非負の値をとるものとします。例えば、 $-1 \equiv 4 \pmod{5}$ です。

System.Random においては、 $j = 55, k = 34, m = 2^ {31} - 1$ です。

Next() を分かりやすいように整理した等価なコードを示します。

// index を +1 して、配列範囲からあふれたらループさせる (% 55 (+1) のイメージ)
if (++_inext >= 56) 
    _inext = 1;
if (++_inextp >= 56)
    _inextp = 1;

// ([_inext] - [_inextp]) mod int.MaxValue のイメージ
_seedArray[_inext] = _seedArray[_inext] - _seedArray[_inextp];
if (_seedArray[_inext] < 0)
    _seedArray[_inext] += int.MaxValue;

return _seedArray[_inext];

さて、擬似乱数の紹介と言えば「周期」が代表的ですが、System.Random の周期は不明です。正しく実装された場合よりも短くなっていると推測されます。 なぜかというと、誠に遺憾ながら System.Random の遷移関数の実装に重大な誤りが複数含まれているためです。

遷移関数の実装ミス (1) - ラグ k の指定

誤りのひとつは、 Random(int seed) コンストラクタの項で意味深な 💀💀💀 を貼っていた _inextp = 21; で、正しくは 31 を代入すべきでした。この値は 2 つめのラグ (漸化式での $k$ ) を指定するものです。

Lagged Fibonacci 法で最大周期を達成するためには、法 $m = 2^ n$ の場合はラグ $j, k$ に $\mathrm{GF}(2)$ 上の原始多項式の指数を用いる必要があります。 (とても雑に言えば、 $j, k$ は適当に決めてはダメで、特別な場合だけ最大周期になるということです。) 例えば、以下に示す $f(x)$ は $\mathrm{GF}(2)$ 上での原始多項式です:

$\begin{aligned} f(x) = x^ {55} + x^ {24} + 1 \end{aligned}$

この指数である 55 と 24 をラグとして用いることで、最大周期を達成できます。 _inextp の正しい初期値 31 は、 55 - 24 から来ています。

そして、System.Random で使われているラグ $(55, 34)$ を多項式で表した $x^ {55} + x^ {34} + 1$ は、原始多項式にはなりません (既約多項式ですらありません。詳細は後述) 。そのため、 System.Random の周期は、達成可能な最大周期よりも短くなっており、かつシードによって周期長が変動してしまう可能性があります。

余談ですが、皮肉にも new Random(int Seed) における 2 つめの初期化フローにおいては、 _inextp = 31; 相当の "正しい" 実装になっています。

遷移関数の実装ミス (2) - 法 m の指定

また、 ~~コピペ元~~ 実装のベースになったと思われる Numerical Recipes の ran3 の実装と比べると、 $m$ の値も $10^ {9}$ から $2^ {31} - 1$ に変更されています。さらに言えば、原典となる TAOCP の ran_array では $2^ {30}$ でした。MSDN の modified version という記述はこの $m$ の変更を指しているものと思われます。 *2

法 $m = 2^ n$ の場合は、下位 $i$ 番目のビット (最下位は 0 とします) の周期が $2^ {i} (2^ {j} - 1)$ となるため、最大周期は $2^ {n-1} (2^ {j} - 1)$ となります。 *3 System.Random でいえば $2 ^ {30} (2 ^ {55} - 1) \approx 2 ^ {85}$ です。

しかし、 System.Random においては $m = 2^ {31} - 1$ で、これは奇数 (もっと言えば素数) なので、以上の議論は適用できません。私が調べた限りでは、少なくとも「 $\mathrm{GF}(2^ {31} - 1)$ 上の」多項式 $f(x) = x^ {55} + x^ {34} + 1$ が原始多項式である必要がありそうです。しかし、 (そもそもラグが誤っているのを無視しても) $\mathrm{GF}(2)$ 上で原始多項式であることと $\mathrm{GF}(2^ {31} - 1)$ 上で原始多項式であることは全く関係がないのと、 $\mathrm{GF}(2^ {31} - 1)$ 上での原始性判定は法の大きさのために困難であると推測されるのとがあります。原始性判定には $(2^ {31} -1) ^ {55} - 1$ の素因数分解を知っている必要がありますが、巨大なため FactorDB でも完全な分解が載っておらず (297 桁の未分解の数を含みます) 、まだ人類は知らない可能性が高いです。

ちなみに $m = 2 ^ {31} - 1$ の場合は、理論上の最長周期は $(2 ^ {31} - 1) ^ {55} - 1 \approx 2 ^ {1705}$ になります。

余談: m = 2 ^ 31 の場合の最下位ビットの周期について

注記: この項での内容は実際の System.Random には適用できません。参考情報程度にお考え下さい。

オリジナルの $m = 2^ {31} - 1$ の場合は解析が厳しいことを上で述べましたが、 +1 して $m = 2^ {31}$ とした場合なら、法が 2 のべき乗になるので本来の $\mathrm{GF}(2)$ 上での理論が適用可能であり、話が比較的簡単になります。

System.Random の $\mathrm{GF}(2)$ 上での多項式は、 Wolfram Alpha 先生曰く以下のように因数分解されます:

$\begin{aligned} f (x) &= x^ {55} + x^ {34} + 1 = f _ 1(x) \cdot f _ 2(x) \cdot f _ 3(x) \\ f _ 1(x) &= x^ {9} + x^ {8} + x^ {7} + x^ {6} + x^ {5} + x^ {3} + 1 \\ f _ 2(x) &= x^ {19} + x^ {18} + x^ {17} + x^ {16} + x^ {14} + x^ {11} + x^ {7} + x^ {6} + x^ {5} + x^ {4} + 1 \\ f _ 3(x) &= x^ {27} + x^ {25} + x^ {23} + x^ {21} + x^ {20} + x^ {17} + x^ {15} + x^ {13} + x^ {11} + x^ {8} + x^ {7} + x^ {6} + x^ {4} + x^ {3} + 1 \\ \end{aligned}$

なお、 $f _ 1, f _ 2, f _ 3$ は原始多項式となることを確認しています。

ここで、内部状態の最下位ビットだけを考えることにすると、遷移式 $x _ i = x _ {i - 55} - x _ {i - 34} \pmod m$ の減算を $x _ i = x _ {i - 55} \oplus x _ {i - 34} \pmod m$ のように XOR に置き換えても意味は変わりません (参考) 。 XOR に置き換えた式から、最下位ビットは LFSR (Xorshift や Mersenne Twister の仲間) と同じ働きをすることが分かります。

LFSR では、 $s$ 回の遷移を $x ^ s \pmod{f}$ と多項式で表すことができます。ここで、 $x = 0$ にすると $s$ にかかわらず = 何度遷移しても 0 になるので、周期が 1 になってしまいます。 (そのため、LFSR ベースの擬似乱数生成器では「全 0 以外で初期化せよ」とよく書かれています。) 今回は法とする多項式が原始ではないので、計算が少しだけ複雑になります。

$f$ を因数分解した $f _ 1, f _ 2, f _ 3$ について、 $x _ 0 \not \equiv 0 \pmod{f _ i}$ を満たす $f _ i$ の位数 ( $x ^ k \equiv 1 \pmod{f _ i}$ を満たす最小の自然数 $k$ ) の最小公倍数が、初期値 $x _ 0$ から始まる系列の周期となります。 ( $f _ i$ は $\mathrm{GF}(2)$ 上の原始多項式なので、位数はそれぞれ $2^ {9} - 1, 2^ {19} - 1, 2^ {27} - 1$ となります。) 逆に考えると、 $x _ 0 \equiv 0 \pmod{f _ i}$ となる $f _ i$ が存在すれば、その分だけ周期が短くなる、ということです。

以上から、全てのとりうる初期状態 $2 ^ {55}$ パターン内での、最下位ビットの周期の内訳を計算したものを下表に示します。最初の 3 列が $x _ 0 \pmod{f _ i}$ が 0 かそれ以外かを示します。「ストリームの本数」は初期値の個数を周期で割ったもので、互いに交わらない乱数列の数を示します。 (ストリーム A で出現する数値は他のストリームでは一切出現せず、その逆も同様です。初期値を変更しない限りストリームが変わることはありません。) なお、最初の 0, 0, 0 の列は 0 で初期化した場合で、一応含めましたが実質ノーカンです。

f1	f2	f3	属する初期値の個数	周期	ストリームの本数
0	0	0	1	1	1
0	0	any	134217727	134217727	1
0	any	0	524287	524287	1
0	any	any	70368609435649	70368609435649	1
any	0	0	511	511	1
any	0	any	68585258497	134217727	511
any	any	0	267910657	267910657	1
any	any	any	35958359421616639	70368609435649	511

割合にすると 99.9998 % ぐらいは周期 $70368609435649 \approx 2^ {46}$ になりますが、それでも実現可能な最大周期 $2^ {55} - 1$ よりかなり小さくなります。また、とても運が悪いとたった 511 回でループする出力列が割り当てられる可能性があります。

具体例を挙げてみます。最下位ビットだけを集めて構成した MinimalBottomSystemRandom において、初期状態 0x6c98d1cf は表の 2 行目に分類されます。そのため、内部状態は 55 ビット分あるにもかかわらず周期として 134217727 が出力されます。

class MinimalBottomSystemRandom
{
    // State の最下位ビット (0x1) が x[0] の最下位ビット, 次のビット (0x2) が x[1] の最下位ビット, ... を表します
    public ulong State { get; private set; }
    const ulong Polynomial = 1 ^ (1ul << 34) ^ (1ul << 55);

    public MinimalBottomSystemRandom(ulong state) => State = state;

    public ulong Next()
    {
        // 55 bit LFSR として構成します
        State <<= 1;
        if ((State & (1ul << 55)) != 0)
            State ^= Polynomial;

        return State;
    }
}

static void Test()
{
    ulong initialState = 0x6c98d1cf;
    var rng = new MinimalBottomSystemRandom(initialState);

    {
        ulong i = 0;
        do
        {
            rng.Next();
            i++;
        } while (rng.State != initialState);

        Console.WriteLine($"loops at {i}");
        // `loops at 134217727`
    }
}

ちなみに、初期値 0x192b762df1 は 3 行目に分類され、周期 524287 となります。

本項の内容は、文献 *4 を参考にしました。ただし、記載している値に誤りがあるような気がします。Fig. 30 の累乗の計算に誤りがあったり、 System.Random の多項式の係数を (55, 21) としていたりします。

遷移関数の性質の要約 - 周期は不明

要約すると、 System.Random には周期を保証できる要素が存在しないという、とても悲しい状態にあります。しかし、これらを修正すると生成する乱数列が変わってしまうため、シードを固定して同じ動作を再現することを期待しているプログラムを壊してしまいます。標準ライブラリのつらさを感じます。

一応 .NET (Framework) のメジャーバージョンをまたいだ場合には、乱数列の再現性は保証しないと明記はされているので、変更することは理論上は可能だとは思いますが、 .NET 5 でも変わらなかったところを見ると…という感じです。

この項は dotnet/runtime に立てられた Issue をもとにしています。別の問題についての言及もあるので、リンクをたどってみると面白いです。

`Next()` は `int.MaxValue` を含まない

以上では、 Next() の遷移関数としての性質について述べました。ここからは、出力関数としての特性を観察します。

MSDN にもあるように、このメソッドは 0 <= x < int.MaxValue (== 2147483647) の範囲の乱数を返します。

A 32-bit signed integer that is greater than or equal to 0 and less than MaxValue.

int.MaxValue は範囲に含まれない のが重要なポイントです。 (この制限は、内部状態が同様に 0 <= x < int.MaxValue の範囲を取っており、それをそのまま出力しているためです。遷移関数の説明での $m$ によります。)

大変面倒なことに、すべてのベースとなるメソッドの値域が 32 bit 全域を埋めないどころかそもそも 2 のべき乗ではありません。そのため、ビット演算・構造を活用した効率化ができないことが問題になります。例えば、4 バイトのランダムな byte[] を得るために BitConverter.GetBytes(random.Next()) などと書いてしまった場合、[3] の最上位ビットが常に 0 だったり、[0] の確率が偏ったりします。

`Next()` の出力パターン数は奇素数である - 範囲加工における偏り

また、出力されうるパターンの数 $2^ {31} - 1 = 2147483647$ は素数です。したがって、 $1, 2^ {31} - 1$ 以外のすべての x に対して 2147483647 % x > 0 を満たします。これが何を意味するかというと、出力をある範囲に加工する際、乗算や剰余を使って安直に変換すると必ず余りが出る = 確率に偏りが生まれることを示します。

例えば、 random.Next() % 401 によって 0 から 400 までの値を得たいとします。 *5 この時、 2147483647 % 401 == 327, 2147483647 / 401 == 5355320 となります。したがって、0 ～ 326 までの範囲は 5355321 / 2147483647 、327 ～ 400 までの範囲は 5355320 / 2147483647 の確率で得られます。このように微々たる差ではありますが、確実に差が生まれます。最も相対的な差が大きくなるのは $x \ge 2 ^ {30}$ の場合で、多いほうと少ないほうで 2 倍確率が異なる (2 / 2147483647 と 1 / 2147483647) ことになります。

一番影響がありそうな事例としては、 random.Next() & 0xff といったコードでしょうか。 $m = 2^ n$ の生成器なら偏りなく各ビットを得ることができますが、そうではないので 0xff が 8388607 / 2147483647 、それ以外が 8388608 / 2147483647 の確率で出現します。

最も、偏りについては「変換手法を工夫すれば」回避可能です。ただ、そんな工夫は System.Random には存在せず、むしろ上述した剰余を使う手法よりも偏りが発生する手法を採用しています。これについては Next(int maxValue) の項で説明します。

`_inextp` は必要か

既に述べた _seedArray の無駄に近い話ですが、これが使用されているのは Next() のほうなのでこちらに書きます。

_inextp は常に (_inext + 20 % 55) + 1 となるので削除可能ではあります。ただ、都度剰余などの計算をしなくてもよくなるので、一概に無駄とは言えないかもしれません。

実際に確かめてみましょう。変更前後で BenchmarkDotNet を使ってベンチマークを取ってみたところ、速度は全くと言っていいレベルで変わりませんでした (変更前 11.72 ± 0.560 ns 、変更後 11.65 ± 0.303 ns 。± は標準偏差) 。 int 1 つ分とはいえメモリを占有することを考えると、やっぱり不要かもしれません。

// 変更後バージョン。 Next() の先頭だけ抜粋: 
if (++_inext >= 56) 
    _inext = 1;
int _inextp = _inext + 21;      // note: ローカル変数なので保存不要
if (_inextp >= 56)
    _inextp -= 55;

余談: 下位ビットの品質について

線形合同法ベースの擬似乱数生成器では「下位ビットは品質が良くない (周期が短い) のでなるべく上位ビットを使え」「範囲に変換するときに剰余 % maxValue を使うな」とよく言われますが、 System.Random ではどうでしょうか。

ふわっとした説明しかできないので申し訳ないのですが、 System.Random においては問題なさそうです (その他もろもろの問題は見なかったことにするなら…) 。というのも、法が $m = 2 ^ n$ ではないので、最上位ビットの情報が剰余を通じて最下位ビットに伝播するため、すべてのビットの周期が最上位ビットと同様になります。

対して $m = 2 ^ n$ の場合は、ビット間の情報はキャリー (繰り上がり・下がり) のみによって伝播するため下位ビットから上位ビットへの一方通行になります。その場合は前述したように最下位から $i$ ビット目の周期は $2^ {i} (2^ {j} - 1)$ となって下位ほど周期が短くなるので、上位ビットを使ったほうが良いことになります。

このように、この問題は線形合同法固有の話ではありません。逆に、線形合同法でも法を奇素数などに設定すればすべてのビットの周期が同じになります。最も、その場合は System.Random と同じようにバイト・ワード全域への均等な出力 (例えば Next() & 0xff) が困難になるので一長一短です。

余談ですが、Xorshift や Mersenne Twister などの LFSR ベースの擬似乱数生成器ではこの問題は生じません。 $n$ ビットの状態を持つとき、すべてのビットの周期が $2 ^ n - 1$ になります。

`int Next(int maxValue)`

0 <= x < maxValue を満たすランダムな x を返します。

実は Next(0) は有効で、 Next(1) と同様に常に 0 を返します。 Next(-1) ではきちんと ArgumentOutOfRangeException がスローされます。

変換関数は、一旦 0.0 <= x < 1.0 の double 値に変換してから maxValue を掛ける、という手法で実装されています。等価な実装を示します:

return (int)(Next() * (1.0 / int.MaxValue) * maxValue);

浮動小数点数の丸めによって確率が偏る

double 上で計算しているので、浮動小数点数の丸めによる計算誤差が生じ、結果として確率が偏る場合があります。

Next(int.MaxValue) の場合が分かりやすいので、例に挙げて説明します。出力範囲としては Next() と同じ 0 <= x < int.MaxValue になりますが、Next() と Next(int.MaxValue) の結果は同一の入力に対して等価になりません。 Next(int.MaxValue) は均等な出力を生成せず、絶対に出力されない値とより多く出力される値が存在します。

この問題は、 Next() * (1.0 / int.MaxValue) * maxValue (double 値) の小数点以下のビットがすべて 1 になっていた場合に発生します。具体例を挙げると、Next() == 4194305 の場合の計算結果は内部表現 0x415000003fffffff (10進 4194304.9999999991) をとり、これを丸めた結果 4194304 と元々の Next() と違う値が返ります。なお、 Next() が前後の数値 4194304, 4194306 を返した場合の計算結果は Next() と同じになります。つまり、Next(int.MaxValue) から 4194305 を得る確率は 0 となり、対して 4194304 を得る確率は 2 / 2147483647 と通常の 2 倍になります。

Next(int.MaxValue) において、入力 ( Next() ) と出力が異なる要素は 9437184 個存在し、全体の約 0.44 % にのぼります。また、出力が異なった結果重複して振り分けられる数も存在することから、絶対出力されない要素が 9437184 個、通常の 2 倍の確率で出力される要素も 9437184 個存在することになります。

このことは目に見える問題となって現れます。より詳細な解析を行っている文献によれば、 Next(int.MaxValue) の返り値が奇数となる確率は約 50.34 % (厳密には 1081081855 / 2147483647) となり、明らかに想定される確率 (50 %) から離れています。

`int Next(int minValue, int maxValue)`

minValue <= x < maxValue の範囲の擬似乱数を出力します。

こちらも、 minValue == maxValue を許容して常に minValue を返します。もちろん minValue > maxValue ならば ArgumentOutOfRangeException をスローします。

範囲が広い場合に確率が偏る

範囲が狭い maxValue - minValue <= int.MaxValue のときは NextDouble() * (maxValue - minValue) + minValue を返すので、 Next(int maxValue) と同様です。

特筆すべきは maxValue - minValue > int.MaxValue の場合で、計算手法が異なります。前提として、元々の範囲より広い範囲に安直に変換すると絶対に出ない値が生じるなど、出力が偏ってしまいます。 (前述した浮動小数点数の丸め誤差とは別です。) 一番単純な例を挙げると、 0 <= x < 2L * int.MaxValue にしようとして x * 2 とすると全て偶数になってしまう、という感じです。余談ですが、この問題は実際に PHP で発生しているようです (雑談の項を参照) 。 System.Random では、この問題を回避するために乱数を継ぎ足す処理が入っています。

maxValue - minValue > int.MaxValue の場合の等価なコードを示します。

int result = Next();            // 0 <= x < int.MaxValue

if (Next() % 2 == 0)            // ※均等ではない
    result = -result;           // ※ 0 が出やすくなる。 -int.MaxValue < x < int.MaxValue

double d = result;
d += int.MaxValue - 1;          // 0.0 <= x < 2.0 * int.MaxValue - 1
d /= 2u * int.MaxValue - 1;     // 0.0 <= x < 1.0

return (int)((long)(d * ((long)maxValue - minValue)) + minValue);

上から順に難癖をつけていきます。まず Next() % 2 == 0 ですが、前述したように Next() のパターン数は $2^ {31} - 1$ 個なので偏りが生じます。 1073741824 / 2147483647 (50.000000023283064 %) の確率で true 、 1073741823 / 2147483647 (49.999999976716936 %) の確率で false になります。

次に result = -result ですが、 0 を入れると値が変わらないことに気づけます。したがって、他の数値の 2 倍の確率で 0 が出現することになります。 (0 が概ね $2^ {-31}$ 、他は $2^ {-32}$ 程度になります。)

残りの d の計算については、 Next(int maxValue) の処理を 0.0 <= x < 2u * int.MaxValue - 1 の範囲で行うイメージです。ただし、ここでは逆数の乗算 d *= 1.0 / constant; ではなく除算 d /= constant; を使用しています。この 2 つの演算はパッと見同じに見えますが、結果は異なります。具体的には、 Next(int.MinValue, int.MaxValue) で呼び出した際、 Next(int maxValue) の項で挙げた浮動小数点数の丸めによる誤差 (result != d * (maxValue - minValue) となる場合) が生じないようです。その代わり、除算なので多少遅くなりそうです。

`double NextDouble()`

0.0 <= x < 1.0 の double 値を返します。

等価な実装を示します。 (実際の実装は Sample() メソッドにあります。) Next(int maxValue) の計算で使われているものと同じです。

return Next() * (1.0 / int.MaxValue);

出力される最大値の誤記

0 より大きい最小の出力値は 4.6566128752457969E-10 ( 内部表現 0x3e000000_00200000 ) 、最大値は 0.99999999953433871 ( 0x3fefffff_ffc00000 ) になります。 MSDN には、取りうる最大値は 0.99999999999999978 ( 0x3fefffff_fffffffe ) であると書いてありますが、Next() が返す最大値を考えれば誤りであることは明らかです。

精度が低い

さて、 Next() の取りうる値は $2^ {31} - 1$ 種類ですから、 NextDouble() が取りうる値も同様に $2^ {31} - 1$ 種類となります。実際に double の取りうる値の数と比べてみましょう。

double において、 0.0 <= x < 1.0 の範囲には $1023 \times 2 ^ {52} \approx 2 ^ {62}$ 個の数値が存在します。 ( 0.0 の内部表現は 0x00000000_00000000 、 1.0 の内部表現は 0x3ff00000_00000000 であることからわかります。) したがって、本来取りうる値の $2 ^ {-31}$ 程度しか埋められないことになります。だいぶ疎です。

最も、全領域をカバーするのは計算コスト的に割に合わないので、実用的には 52 ～ 53 bit 精度で妥協する場合が多いです。この 52 は、 double 型の仮数部の桁数から来ています。 (+1 はケチ表現の分です。) 例えば、以下の手法では 53 bit 精度の double 値が得られます。

// ここでの Next() は ulong 型 (64bit) 全域のランダムな値を返すものとします
return (Next() >> 11) * (1.0 / (1ul << 53));

何故素直に Next() * (0.5 / (1ul << 63)) としないのかというと、乗算の前に Next() が double にキャストされて結局 53 bit 精度になってしまうためです。

範囲加工における注意 - 最大値を含む可能性

Next() には範囲を指定するオーバーロードがありますが、 NextDouble() にはありません。 MSDN を見ると、以下の式で手動変換するとよい、とあります。

NextDouble() * (maxValue - minValue) + minValue

このコードには注意すべき点があります。 NextDouble() は 0 <= x < 1 なので、この式も minValue <= x < maxValue となりそうですが、 minValue の絶対値が大きく maxValue - minValue が相対的にかなり小さい場合などに x == maxValue となる場合があります。具体例としては (987654321, 987654444) が挙げられます。余談ですが、Java の java.util.Random#doubles(double, double) では、このような場合に Math.BitDecrement() 相当のメソッドを使って 1 つだけ小さくすることで回避しているようです (実装)。

通常はあまり気にしなくても良いかと思いますが、 x == maxValue になるとゼロ除算や予期しない ∞ ・ NaN が発生しうる場合 (Math.Log() など) には注意が必要です。

`NextBytes(byte[] buffer), NextBytes(Span<byte> buffer)`

byte[] または Span<byte> 全域をランダムな値 ( 0x00 <= x <= 0xff ) で埋めます。 Span<byte> のほうは .NET Standard 2.1 から使えます。

等価な実装を示します。処理的にはどちらも同様に Next() の下位 8 ビットを使用して埋めていきます。

for (int i = 0; i < buffer.Length; i++)
    buffer[i] = (byte)(Next() & 0xff);

`0xff` が生成される確率が低い

さて、何度か述べたように Next() は $2^ {31} - 1$ 通りなので、 Next() & 0xff (言い換えれば Next() % 256) は均等になりません。 0xff が 8388607 / 2147483647 、それ以外はそれぞれ 8388608 / 2147483647 の確率となり、 0xff が出現する確率がわずかに低くなります。

`Next()` 1 回あたり 1 byte しか利用しない非効率な実装

本来 4 byte に近い情報量のある Next() 1 回あたり 1 byte 分しか生成しない、とても効率の悪い実装となっています。 Next() の値域問題もあるので 4 byte フルに使うわけにはいきませんが、少なくとも 3 byte 分は使えるので 3 倍近く速くできそうな気がします。

外部実装・その他の問題

以上で System.Random に実装されているメソッドの紹介は終わりです。あとは MSDN などを見ながら、さらなる重箱の隅をつついてきます。

スレッドセーフではない

System.Random はスレッドセーフではありません。乱数生成のたびに内部状態を更新する必要があるので、それはそうです。 MSDN のマルチスレッド処理のサンプルを見ると、乱数を取得する箇所を lock ステートメントで囲って対応しています。

lock で囲む以外の対応策としては、スレッドごとに System.Random インスタンスを生成する手法があります。無理にロックしなくていいのでこちらのほうが速くなりそうな気がします。 (試せてはいませんが……)

ただ、前述したように初期状態のパターンが $2 ^ {31}$ に限定される問題 (生成する乱数列のバリエーションが少ない、最悪衝突する可能性) があったり、.NET Framework ならより深刻なシードが同一になる問題があったりするので、難しいところがあります。 xoroshiro128+ などの jump() や、 java.util.SplittableRandom の split() のようなものをサポートしていればこの辺りが解決できるのですが、残念ながら System.Random では非対応です。

再現が難しい / シリアライズできない

System.Random では、同じ乱数列を再現したい場合 new Random(int Seed) を使用することになります。しかし難点がいくつかあります:

内部状態の情報量に対して、初期状態がしかない
- シードだけを保存する場合、とりうる状態 (→出力される乱数列) がかなり制限される
途中から再開することができない
- 例えば、乱数を $2^ {30}$ 個消費した実験の続きをしたい場合、もう一度 $2^ {30}$ 個消費しなければならない

内部状態の保存 (シリアライズ) ができれば簡単に途中から再開できるようになりますが、 System.Random はそのままではシリアライズできません。内部状態は private フィールドで、属性も特についていないためです。 protected なら継承すれば合法的に触れたかもしれませんが…

リフレクションを使えばできないことはないですが、ちょっと面倒です。保存すべきフィールドは int[] _seedArray, int _inext, int _inexp の 3 つです。以下に状態を読み書きするサンプルを示します。テストのために内部状態を適当に書き換えていますが、実際に保存用途に使う場合はスキップしてお好みの入出力コードを書いてください。ただ private フィールドなので外部アクセスは本来想定外の操作であり、将来的な動作は一切保証されません。

var rng = new Random();

var bindingFlag = BindingFlags.NonPublic | BindingFlags.Instance;
var stateField = typeof(Random).GetField("_seedArray", bindingFlag);
var inextField = typeof(Random).GetField("_inext", bindingFlag);
var inextpField = typeof(Random).GetField("_inextp", bindingFlag);

// 内部状態からの読み込み
var state = stateField.GetValue(rng) as int[];      // int[56]
var inext = (int)inextField.GetValue(rng);
var inextp = (int)inextpField.GetValue(rng);

// 書き込みテスト用(次の出力で 0 が出るようにする + ラグ指定修正)
state[2] = 12345;
state[33] = 12345;
inext = 1;
inextp = 32;

// 内部状態への書き込み
stateField.SetValue(rng, state);
inextField.SetValue(rng, inext);
inextpField.SetValue(rng, inextp);

Console.WriteLine(rng.Next());      // `0`

ところで、内部状態全体の保存が面倒だからと言って、シード値だけを保存して乱数のセーブポイントを通るたびに (頻繁に) new Random(seed) しなおすような実装はくれぐれも行わないようにしてください。繰り返しになりますが初期状態は $2 ^ {31}$ パターンしかないので、出力列が著しく限定され、最悪の場合意図せず同一の乱数列を出力してしまう可能性があります。

拡張性

継承による実装の拡張

MSDN を見ると、 System.Random を継承して別の擬似乱数生成器を実装したり、NextBytes(byte[] bytes, int min, int max) などの追加メソッドを生やしたりする手法が紹介されています。

しかし、継承した場合もれなく System.Random の内部状態 280 byte がついてきます (実測値。実装によって変わるかもしれません) 。もし別の擬似乱数生成器を実装するなら無駄になってしまいます。また、Sample() メソッドを上書きすれば大体 OK 、といったことが書かれていますが、 double ベースの計算になるので上述した問題が生じます。といってもすべての public メソッドが virtual なので、すべて自前実装で置き換えてしまえば問題はありません。

アルゴリズムを差し替えたり出力関数を柔軟に追加したりすることを考えると、 interface IRandom を定義しておいて、出力関数などは拡張メソッドにする (今なら interface のデフォルト実装でしょうか？) 手法が思いついたりしますが、 System.Random の実装当時 (.NET Framework 1.1 → C# 1.2) には拡張メソッドが存在しないのでダメです。歴史の重みです。

いくつかの擬似乱数生成器を実装している Math.NET Numerics を覗いてみると、擬似乱数生成器側は class RandomSource を基底として遷移関数と基礎的な出力関数を実装し、複雑な変換関数 (正規分布など) の実装は class Normal といった別クラスに分離して new Normal(new Xoshiro256StarStar()).Sample() といった感じで利用できるようにしているようです。これが綺麗な気がします。

各種分布への変換

正規分布への変換

上の項で分布の話が出たので、一様分布以外では多分最もメジャーな正規分布への変換についても触れておきます。

Java には java.util.Random#NextGaussian がありますが、我らが System.Random にはもちろん実装されていないので、自力で対応する必要があります。

ざっくり Polar Method で実装すると以下のようになります。

/// <summary>
/// 平均 0, 標準偏差 1 の正規分布に従う独立な乱数のペアを取得します
/// </summary>
public (double x, double y) NextGaussian()
{
    double u, v, s;
    do
    {
        u = random.NextDouble() * 2 - 1;
        v = random.NextDouble() * 2 - 1;
        s = u * u + v * v;
    } while (s >= 1 || s == 0);

    s = Math.Sqrt(-2 * Math.Log(s) / s);
    return (u * s, v * s);
}

x, y は独立なので、どちらを使っても・ペアで使っても大丈夫です。1 つだけ要る場合は単に捨てればよいです。もし平均と標準偏差を変えたければ returnValue * stdev + mean としてください。

正規分布への変換には Box-Muller Transform もありますが、Polar Method のほうが若干高速です。雑にベンチマークを取ったところ、手元の環境では Polar Method が 44.90 ± 1.041 ns 、Box-Muller Transform が 56.32 ± 0.702 ns (± は標準偏差) でした。 Box-Muller Transform の微妙な利点としては、工夫すると消費する乱数が 2 個固定になる点が挙げられます。また、より高速な実装手法として Ziggurat algorithm がありますが、テーブルが必要だったり実装が大変だったりするのでここでは取り上げません。大量に正規分布乱数が必要であれば検討するとよいと思います。

その他の分布への変換

Math.NET Numerics を活用するとよいと思います。 (丸投げ)

数学的な確率分布以外では、円|球の中|外周にランダムに配置する計算が (特に Unity で) よく使われます。こちらの説明が分かりやすいです。

シャッフル

MSDN のサンプル実装は非効率

コレクションの要素の順序をランダムに入れ替える「シャッフル」について検討します。 System.Random には組み込みのシャッフル関数は存在しません。ただ、 MSDN にはカードのシャッフルのサンプル実装があります。

var deck = new Card[52] { ... };

// シャッフル用インデックス
var order = new double[deck.Length];
for (int i = 0; i < order.Length; i++)
    order = random.NextDouble();

// order をキーとして deck をソート
Array.Sort(order, deck);

この手法には何点か問題があります:

NextDouble() の出力値が同じ要素が存在した (衝突した) 場合の動作が不定です。
- MSDN の Array.Sort を見た限りではイントロソート (不安定) のようですので、衝突した場合どちらに偏るかは場合によります。
  - ちなみに安定ソート ( Enumerable.OrderBy ) なら確実に偏ります。元々先頭に近いほうが必ず前に来ます。
- 直感的には奇跡でも起こらない限りあり得ないのでは、と思いそうですが、誕生日のパラドックスより 54562 個程度の要素に対して実行すると概ね 1/2 の確率でどれかが衝突します。
時間的に非効率です。一般にソートの計算量は $O(n \log n)$ ですが、より効率的に $O(n)$ でシャッフルする手法が存在します。
空間的に非効率です。同じ要素数の double[] を確保する必要があります。

素直に Fisher-Yates シャッフルを実装しましょう。例では直書きしましたが、 IList<T> の拡張メソッドにしたりすると便利です。

var deck = new Card[52] { ... };

for (int i = deck.Length - 1; i > 0; i--)
{
    int r = rand.Next(i + 1);
    (deck[i], deck[r]) = (deck[r], deck[i]);
}

副作用をなくして IEnumerable<T> に対応する場合は、取り出し版のアルゴリズムで 0 ～ n - 1 の連番をシャッフルした配列を作成して、 source.ElementAt(randomArray[i]) といった感じで参照する手法が考えられます。

余談: ランダム値のソートによるシャッフルの注意点

話はそれますが、ランダムな値をキーとしてソートする場合、横着して以下のように書くと例外 (ArgumentException など) がスローされる場合があります。

// DANGER!
Array.Sort(deck, (_, _) => random.Next(-1000, 1000));

そもそもソートは値を一定の順番に並べるためのものなので、比較するたびに大小が変わるような一貫性のない操作を行うと破綻してしまうためです。しかも性質が悪いのが「確実に例外がスローされるわけではない」点で、初回は成功したからといって油断していると後で刺されたりします。

爆発するかどうかはソートの形式によります。 (a, b) => Next() のように 2 引数を取って比較するもの (Comparison<T> や IComparer<T> をとるものなど) は危険で、ランダムに例外が投げられます。対して、 a => Next() のように 1 引数をキーとして取るもの ( Enumerable.OrderBy(Func<TSource, TKey>) など) は、問題ないこともあります。内部実装として、メソッドを 1 回評価して配列などに格納 → デフォルトの比較関数でソート、というようになっていれば、比較自体は一貫して行えるため問題ありません。

乱数性テスト

出力列の品質、言い換えれば「出力列が乱数らしい振る舞いをしているか」を調べてみます。遷移・出力関数の各所に問題があるのは今まで述べたとおりですが、テストを通しても検出できることを示す、という位置づけにします。

「テスト」は、出力された乱数列に統計的検定や構造解析を行うことで、何らかの乱数らしからぬ構造が見出せるか、を調べるシステムをいいます。「良い乱数であることを示す」よりは、「特定条件下で明らかに悪いことを示す」ようなイメージです。ただ、たくさんのテストを実施して問題が検出されなければ、品質が良い (真の乱数列に近い性質を持つ) とみなせそうです。 *6

PractRand のテストで速やかに失敗する

強力で汎用的な乱数性テストツールとしては PractRand が有力です。自動で複数種のテストをかけることができ、検出力が高く誤検知 (真の乱数を入力してもテストに失敗すること) もほぼありません。私は PractRand でデフォルトの最大値である 32 TB の出力を検定し、 FAIL となるテストが 1 つも存在しないことを前提として擬似乱数の設計・選定を行っています。

今回は、より多くのテストを行う設定の -te 1 -tf 2 -tlmaxonly を指定しました。出力には、 Next() では 32 bit 全域の出力ができないため NextBytes() 相当の実装を使用しました。大抵のテストツールでは基本的に 32 or 64 bit 全域を出力するメソッドがあることを前提としているので、そのままではテストにかけられないこと自体がとても厳しいです。というのも、この変換がテスト結果に影響を及ぼす可能性があるので、妥当かどうか微妙なところがあるのです…

結果を以下に示します。疑わしいテスト結果だけが表示されており、成功したテストは ...and から始まる行にまとめて書いてあります。結果は Evaluation のところに注目してください。雰囲気としては以下のような感じです。

usual : まれによくある
- これはあまり問題になりません。別のシードでは起きなかったりする微妙なエラーです
suspicious : だいぶ怪しい
- 大抵の場合次の length で FAIL します
FAIL : 確実にアウト
- 明らかに何らかの欠陥があります

テスト結果 (長いので折り畳み - クリックすると開きます)

RNG_test using PractRand version 0.94
RNG = .netrandom, seed = 0xb2540eef
test set = expanded, folding = extra

rng=.netrandom, seed=0xb2540eef
length= 32 megabytes (2^25 bytes), time= 2.7 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/8]BCFN(0+1,13-5,T)          R= +20.1  p =  4.8e-8   very suspicious
  ...and 936 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 64 megabytes (2^26 bytes), time= 14.9 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/8]BCFN(0+1,13-5,T)          R= +41.1  p =  2.9e-16    FAIL !
  [Low1/8]DC6-6x2Bytes-1            R=  +6.5  p =  4.7e-4   unusual
  ...and 1006 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 128 megabytes (2^27 bytes), time= 30.8 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/8]BCFN(0+1,13-4,T)          R= +79.5  p =  1.0e-34    FAIL !!!
  [Low1/8]DC6-6x2Bytes-1            R=  +8.7  p =  3.4e-5   mildly suspicious
  [Low1/8]DC6-5x4Bytes-1            R=  +8.9  p =  1.5e-5   mildly suspicious
  ...and 1078 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 256 megabytes (2^28 bytes), time= 48.0 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/8]BCFN(0+0,13-3,T)          R= +44.0  p =  1.2e-20    FAIL !!
  [Low1/8]BCFN(0+1,13-3,T)          R=+129.4  p =  4.8e-61    FAIL !!!!
  [Low1/8]DC6-6x2Bytes-1            R= +14.5  p =  1.1e-7   very suspicious
  [Low1/8]DC6-5x4Bytes-1            R= +20.1  p =  5.5e-13    FAIL
  [Low1/32]FPF-14+6/16:cross        R=  +5.7  p =  4.8e-5   unusual
  ...and 1146 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 512 megabytes (2^29 bytes), time= 68.1 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/8]BCFN(0+0,13-3,T)          R= +91.4  p =  4.8e-43    FAIL !!!
  [Low1/8]BCFN(0+1,13-3,T)          R=+257.6  p =  1.0e-121   FAIL !!!!!
  [Low1/8]DC6-6x2Bytes-1            R= +30.8  p =  4.1e-20    FAIL !!
  [Low1/8]DC6-5x4Bytes-1            R= +36.5  p =  2.3e-24    FAIL !!
  ...and 1216 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 1 gigabyte (2^30 bytes), time= 96.9 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(0+3,13-1,T)                  R= +11.4  p =  1.4e-5   mildly suspicious
  BCFN(0+4,13-1,T)                  R= +13.0  p =  1.7e-6   suspicious
  [Low1/8]BCFN(0+0,13-2,T)          R=+388.0  p =  6.9e-198   FAIL !!!!!!
  [Low1/8]BCFN(0+1,13-2,T)          R=+465.4  p =  1.9e-237   FAIL !!!!!!
  [Low1/8]BCFN(0+2,13-3,T)          R= +18.9  p =  9.2e-9   very suspicious
  [Low1/8]DC6-6x2Bytes-1            R= +57.2  p =  5.1e-32    FAIL !!!
  [Low1/8]DC6-5x4Bytes-1            R= +62.9  p =  1.5e-35    FAIL !!!
  [Low4/16]BCFN(0+2,13-2,T)         R= +26.9  p =  2.5e-13    FAIL
  ...and 1286 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 2 gigabytes (2^31 bytes), time= 146 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(0+3,13-0,T)                  R= +27.2  p =  4.7e-14    FAIL
  BCFN(0+4,13-1,T)                  R= +28.3  p =  1.0e-14    FAIL
  [Low1/8]BCFN(0+0,13-1,T)          R= +1303  p =  2.3e-701   FAIL !!!!!!!
  [Low1/8]BCFN(0+1,13-1,T)          R=+819.6  p =  8.3e-441   FAIL !!!!!!!
  [Low1/8]BCFN(0+2,13-2,T)          R= +47.7  p =  5.7e-24    FAIL !!
  [Low1/8]DC6-6x2Bytes-1            R=+110.9  p =  1.1e-52    FAIL !!!!
  [Low1/8]DC6-5x4Bytes-1            R=+117.3  p =  2.9e-74    FAIL !!!!
  [Low4/16]BCFN(0+2,13-1,T)         R= +42.0  p =  4.8e-22    FAIL !!
  [Low4/16]BCFN(0+3,13-1,T)         R= +13.9  p =  5.8e-7   suspicious
  ...and 1359 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 4 gigabytes (2^32 bytes), time= 241 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(0+3,13-0,T)                  R= +57.5  p =  2.8e-30    FAIL !!!
  BCFN(0+4,13-0,T)                  R= +50.5  p =  1.6e-26    FAIL !!
  [Low1/8]BCFN(0+0,13-1,T)          R= +2589  p =  1e-1393    FAIL !!!!!!!!
  [Low1/8]BCFN(0+1,13-1,T)          R= +1659  p =  9.2e-893   FAIL !!!!!!!
  [Low1/8]BCFN(0+2,13-1,T)          R= +82.8  p =  4.8e-44    FAIL !!!
  [Low1/8]BCFN(0+4,13-2,T)          R= +10.7  p =  4.6e-5   unusual
  [Low1/8]DC6-6x2Bytes-1            R=+204.2  p =  6.8e-109   FAIL !!!!!
  [Low1/8]DC6-5x4Bytes-1            R=+219.3  p =  1.7e-100   FAIL !!!!!
  [Low1/32]BCFN(0+6,13-4,T)         R= +12.2  p =  2.9e-5   unusual
  [Low4/16]BCFN(0+2,13-1,T)         R= +89.6  p =  1.0e-47    FAIL !!!
  [Low4/16]BCFN(0+3,13-1,T)         R= +24.5  p =  1.1e-12    FAIL
  ...and 1433 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 8 gigabytes (2^33 bytes), time= 497 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(0+3,13-0,T)                  R=+115.8  p =  2.0e-61    FAIL !!!!
  BCFN(0+4,13-0,T)                  R= +88.9  p =  4.5e-47    FAIL !!!
  [Low1/8]BCFN(0+0,13-0,T)          R= +5800  p =  9e-3101    FAIL !!!!!!!!
  [Low1/8]BCFN(0+1,13-0,T)          R= +2740  p =  8e-1465    FAIL !!!!!!!!
  [Low1/8]BCFN(0+2,13-1,T)          R=+169.1  p =  1.7e-90    FAIL !!!!!
  [Low1/8]BCFN(0+4,13-1,T)          R= +15.6  p =  7.9e-8   very suspicious
  [Low1/8]DC6-6x2Bytes-1            R=+386.9  p =  1.0e-179   FAIL !!!!!!
  [Low1/8]DC6-5x4Bytes-1            R=+415.2  p =  1.1e-253   FAIL !!!!!!
  [Low4/16]BCFN(0+2,13-0,T)         R=+156.1  p =  5.4e-83    FAIL !!!!
  [Low4/16]BCFN(0+3,13-0,T)         R= +34.5  p =  5.3e-18    FAIL !
  ...and 1530 test result(s) without anomalies

rng=.netrandom, seed=0xb2540eef
length= 16 gigabytes (2^34 bytes), time= 976 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(0+3,13-0,T)                  R=+234.3  p =  8.5e-125   FAIL !!!!!
  BCFN(0+4,13-0,T)                  R=+172.1  p =  1.4e-91    FAIL !!!!!
  BCFN(0+6,13-0,T)                  R= +12.4  p =  3.6e-6   unusual
  [Low1/8]BCFN(0+0,13-0,T)          R=+11611  p =  1e-6207    FAIL !!!!!!!!
  [Low1/8]BCFN(0+1,13-0,T)          R= +5513  p =  2e-2947    FAIL !!!!!!!!
  [Low1/8]BCFN(0+2,13-0,T)          R=+273.9  p =  5.6e-146   FAIL !!!!!
  [Low1/8]BCFN(0+3,13-0,T)          R=  +9.6  p =  1.1e-4   unusual
  [Low1/8]BCFN(0+4,13-1,T)          R= +22.3  p =  1.7e-11   VERY SUSPICIOUS
  [Low1/8]DC6-6x2Bytes-1            R=+735.4  p =  1.2e-265   FAIL !!!!!!
  [Low1/8]DC6-5x4Bytes-1            R=+780.1  p =  6.6e-429   FAIL !!!!!!!
  [Low4/16]BCFN(0+2,13-0,T)         R=+342.4  p =  1.2e-182   FAIL !!!!!!
  [Low4/16]BCFN(0+3,13-0,T)         R= +70.4  p =  3.7e-37    FAIL !!!
  ...and 1625 test result(s) without anomalies

64 MB まで検定した時点で BCFN テストに、 256 MB から DC6 テストに FAIL しました。ほぼ瞬殺です。無事失敗したのと時間がかかるのとで 16 GB で中断しました。

参考として他の擬似乱数生成器はどうなのかというと、有名どころとしては線形合同法 (LCG) ・ Mersenne Twister ・ Xorshift ・ xoroshiro128+ は FAIL していて、PCG ・ xoshiro256** ・splitmix (≒ java.util.SplittableRandom) は 32 TB でもパスしています。詳しくはこちらの "Failed Test" 列を参照してください。 hwd 以外は PractRand のテストです。

Hamming-weight dependencies

次に、xoroshiro/xoshiro の作者が作った Hamming-weight dependencies テスト (略称 hwd) *7 を試してみます。これは TinyMT において BRank 系以外のテストで初めて問題を検出できた実績のあるテストです。こちらは 1 PB (ペタバイト！) テストして p < 1e-20 にならなければパスしたといえます。

こちらでは試しに Next() 相当の実装を使用しました。また、HWD_BITS=32 (出力のワードサイズが 32 bit) としました。

結果の一部を以下に示します。本項執筆時点では 25 TB まで検定して p = 0.761 なので、意外なことに問題は検出されていなさそうです。最後までやると 2-3 週間ぐらいかかるのでこの辺りで許してください…

mix3 extreme = 1.51920 (sig = 00000100) weight 1 (16), p-value = 0.89
mix3 extreme = 2.67035 (sig = 00010010) weight 2 (112), p-value = 0.573
mix3 extreme = 3.41468 (sig = 00101200) weight 3 (448), p-value = 0.249
mix3 extreme = 3.23865 (sig = 12002010) weight 4 (1120), p-value = 0.74
mix3 extreme = 3.64975 (sig = 12222100) weight >=5 (4864), p-value = 0.721
bits per word = 32 (analyzing bits); min category p-value = 0.249

processed 2.5e+13 bytes in 4.14e+04 seconds (0.604 GB/s, 2.174 TB/h). Sun Jan  3 00:05:51 2021

p = 0.761
------

ジャンプ

ジャンプ機能がない

擬似乱数生成器によっては、「ジャンプ」機能を備えているものがあります。具体例を挙げると、 xoroshiro128+ では $2 ^ {64}$ 回の next() 呼び出しと等価な操作を速やかに行うことができる jump() 関数が用意されています。この機能は並列実行などで複数の乱数列が必要な場合に便利で、 new Random() ではなく var child = original.Jump(); のようにしてインスタンスを生成することで、先の例でいえば $2 ^ {64}$ 回 next() を呼び出すまでは重複しない保証のある乱数列を得ることができます。

さて System.Random ではどうなのかというと、少なくとも用意はされていません。設計するのもなかなか難しそうな気がします。 $\mathrm{GF}(2 ^ {31} - 1)$ 上での多項式計算を頑張ればできるかもしれませんが、 Next() の項で述べたように原始か分からなかったり、そもそも計算がよくわからなかったりするので私は何も分かりません。

余談: 「正しく」実装されていた場合のジャンプの設計

注記: 例によって System.Random には直接適用できません。

また机上の空論になってしまいますが、もしラグ指定が正しく法も 2 のべき乗である $k = 24, m = 2 ^ {31}$ だったのならば、限定的ながら可能です。

先に述べたように、 $m = 2 ^ n$ ならば最下位ビットだけを見ると LFSR として扱うことができます。 LFSR でのジャンプ手法を用いることで、「少なくとも $s$ ステップ以上の距離がある状態」を作り出すことが可能です。「少なくとも」というのは最下位ビット以外については制御できないためで、実際の距離は $(2 ^ {55} - 1)i + s \; (i \ge 0)$ という感じになると思います、たぶん…

LFSR のジャンプは、 $j(x) = x ^ s \pmod{f(x)}$ を事前に計算しておいて ( $f(x)$ は遷移関数の特性多項式です) 、実行時に $j(x)$ の $x$ に現在の状態を代入することで行います。例えば、 $2 ^ {54}$ ステップのジャンプ多項式 $j _ {2 ^ {54}}(x)$ は、 xoroshiro128+ の JUMP[] と同じフォーマットなら 0x0010010010200200 として表せます。あとは xoroshiro128+ の jump() と同じ雰囲気で実装すれば少なくとも $2 ^ {54}$ ステップ以上の距離がある状態を作り出すことができます。

ただ、変わるのは最下位ビットだけなので、この手法だけでは出力自体はしばらく元の出力列とほとんど同じになることが予想されます。難儀です。とっても雑に見積もると、キャリーが 1 bit 上に伝播するのに 1 周かかる (= 55 回) x 31 bit なので 1705 回ぐらい乱数を読み飛ばすことで目に見える相関はなくなりそうな気がします。 (気がするだけで確たる根拠はないです…) もちろん重い処理なので手軽ではないですが…

「ガチャシミュレーション」における異常な振る舞い

本項の検討は、 dotnet/runtime の Issue #40490 およびオリジナルの StackOverflow への投稿に基づきます。ここまでの問題は誤差レベルのものが多かったですが、これは実用的なユースケースで目に見える問題が生じます。

System.Random を用いて、単発で確率 $p$ で当たるガチャを当たるまで引く (!) とき、最初に当たるまでに引いた回数 $n$ はどんな分布になるかのシミュレーションをしてみましょう。数学的には幾何分布というらしいです。

以下のように実装してみます。ガチャで当たるまでの回数の分布を、100 万回ぐらいシミュレーションしてみます。

var rand = new Random();

double p = 0.05;        // 当たる確率
var distribution = new int[256];        // (i + 1) 回目で初めて当たりを引いた回数

for (int i = 0; i < 1 << 20; i++)
{
    int done = 0;       // 引いた回数(-1)
    while (rand.NextDouble() >= p)      // 当たるまで回す!!
        done++;

    distribution[Math.Min(done, distribution.Length - 1)]++;
}

// 出力
for (int i = 0; i < distribution.Length; i++)
    Console.WriteLine($"{i + 1,3}: {distribution[i]}");

結果をグラフに書くと以下のようになります。

f:id:andantesoft:20201230224030p:plain

55 (回目に引いたときに初めて当たった回数) に不自然なピークがあるのが分かります。拡大して見てみましょう。ここで、追加の「理論値」は幾何分布の確率質量関数から求めています。

f:id:andantesoft:20201230224036p:plain

あるべき数の大体 1/2 ぐらいになっていそうです。なお、この現象は確率 $p$ に依存せず発生するようです。

原因を考えてみましょう。遷移関数についての今までの説明から、55 回目の遷移・出力は $x _ {55} = x _ {0} - x _ {34} \pmod{ m}$ と表せます。確率 $p$ は簡単のため整数で定義することとし、 $1 \le p \lt m$ として $x _ i \lt p$ のときに「成功」とします。

まず、式中の $x _ 0, x _ {34}$ がとりうる値の範囲を考えます。 $x _ {34}$ については、「55 回目に初めて成功した」すなわち「54 回目まではすべて失敗した」前提条件から、 $x _ {34} \ge p$ です。 $x _ 0$ は、現在の 55 回前ですから「直前の試行の最後の出力」を指します。 (一番最初の試行でない限り) 最後は成功で終わっているので、 $x _ {0} \lt p$ です。

以上の条件下で $x _ {55} \lt p$ を満たすには、 $x _ {0} - x _ {34} + m \lt p$ である必要があります。ここで $m$ を足したのは mod を打ち消すためです ( $x _ {0}, x _ {34}$ の値域から、 $x _ {0} - x _ {34}$ は必ず負値になります) 。この 3 つの関係式を満たす $x _ {0}, x _ {34}$ ペアの個数は $p(p-1) / 2$ 個となります。 $x _ 0 \lt p, x _ {34} \ge p$ の範囲条件を満たすペアの総数 (上記関係式は満たさなくてもよい場合) は $p(m - p)$ 個ですから、この範囲内で成功する確率は $(p(p-1)/2)/(p(m-p)) = (p-1) / 2(m - p)$ となります。本来の確率は $p/m$ ですから、 $p$ がある程度小さければ本来の確率の大体 1/2 ぐらいになり、実験と一致しそうです。

イメージしやすいように図に描いてみましょう。 $x _ {0}, x _ {34}$ 全域 (グラフ全体) に対して、 $x _ {0} \lt p, x _ {34} \ge p$ を満たすのは赤い部分です。色が濃い部分が成功、薄い部分は失敗を表します。 (色が濃い部分 / グラフ全体) に対して、 (赤色の濃い部分 / 赤色全体) の割合が半分ぐらいになっているのが分かるかと思います。

f:id:andantesoft:20201230224120p:plain

さて、原因は分かったので対策を考えてみましょう。直前の試行の最後の出力を参照して $x _ 0 \lt p$ となっている、すなわち直前の試行の影響を受けているのがよくないので、例えば試行が終わった直後にランダムな個数だけ乱数を消費してみるのはどうでしょうか。あまり行儀は良くないですがサンプルということで許してください…

// 抜粋
for (int i = 0; i < 1 << 20; i++)
{
    int done = 0;       // 引いた回数(-1)
    while (rand.NextDouble() >= p)      // 当たるまで回す!!
        done++;

    distribution[Math.Min(done, distribution.Length - 1)]++;

    // 追加コード: ランダムに 0 - 255 個消費する
    int skipAmount = Guid.NewGuid().ToByteArray()[15];
    for (int k = 0; k < skipAmount; k++)
        rand.Next();
}

結果をグラフに書くと以下のようになり、本項の現象がみられなくなっています。

f:id:andantesoft:20201230224033p:plain

一応補足ですが、乱数の消費回数を決定するのに rand.NextInt(256) を使わなかったのには理由があります。自分自身の乱数を使用すると乱数列が重複する可能性が上がってしまうためです。説明を簡単にするために 1 byte のカウンタを擬似乱数生成器と言い張ることにしましょう:

int rand() => (x = (x + 1) % 256);

ここで、 x == 128 のときは rand() => 129 となるので 129 回消費して x == 2 になります。また、 x == 0 のときは rand() => 1 となり 1 回消費して x == 2 になります。これらより、別の内部状態からランダムスキップによって同じ内部状態になってしまい、以降の出力列が重複してしまうことが分かります。 rand() とは関連性のない別系列の乱数 (ここでは Guid.NewGuid() の最後のバイト) を用いることで、この問題を回避している、というわけです。

話を戻してまとめに入ります。遷移・出力関数の実装がシンプルなため、「内部での状態更新処理」と「外部での乱数の利用処理」が嚙み合ってしまった場合このように問題が表出します。とはいえ「内部処理を把握したうえで注意して書け」というのは現実的ではありません。強いて対策を挙げるとすれば、別の擬似乱数生成器を採用するなどでしょうか……？

乱数列の過去と未来を知る

未来の予測 (内部状態復元攻撃)

出力から内部状態を復元して、将来の乱数の予想を試みます。

一番簡単な Next() について検討します。 Next() では内部状態をそのまま返しているので、連続する 55 個の出力列を記録できれば、それがそのまま内部状態 _seedArray になります。あとは _inext = 0; _inextp = 21; としてインスタンスに設定すれば出力列が再現できます。

// コピー元
var rand = new Random();

// 出力の記録 (内部状態に合わせるため [0] は飛ばす)
int[] output = new int[56];
for (int i = 1; i < output.Length; i++)
    output[i] = rand.Next();

// クローンの生成・設定
var clone = new Random();
SetInternalState(clone, output, 0, 21);

// 検証
for (int i = 0; i < 256; i++)
    Debug.Assert(rand.Next() == clone.Next());

// ---

// リフレクションで内部状態を設定する
static void SetInternalState(Random rng, int[] state, int inext, int inextp)
{
    var bindingFlag = BindingFlags.NonPublic | BindingFlags.Instance;
    var stateField = typeof(Random).GetField("_seedArray", bindingFlag);
    var inextField = typeof(Random).GetField("_inext", bindingFlag);
    var inextpField = typeof(Random).GetField("_inextp", bindingFlag);

    stateField.SetValue(rng, state);
    inextField.SetValue(rng, inext);
    inextpField.SetValue(rng, inextp);
}

また、もし内部状態全域が分からなくとも、遷移式 $x _ {i} = x _ {i-55} - x _ {i-34} \pmod{2^ {31} - 1}$ より 55 回前と 34 回前の出力を得ることによって、次に出力される値を求めることができます。

NextDouble() の場合は、 (int)Math.Round(value * int.MaxValue) などで整数に逆変換すれば、 Next() と同じ手法が適用できます。

その他の範囲に加工するもの (Next(int maxValue), NextBytes など) を情報源として内部状態の完全な復元を行う手法については、パッとは思いつきません。Lagged Fibonacci 法に対して適用できるアルゴリズムがあれば、それが使えるかもしれません。

別の手法としては、初期状態が高々 $2^ {31}$ パターンしかないことから、同じ乱数列を出力する系列を $2^ {31}$ パターンの中から力ずくで探索する手法が挙げられます。初期化してからそれほど使用されていない (Next() 呼び出し回数が少ない) と推測される場合に特に有効です。また、 .NET Framework では初期化は起動時間ベースなので、全空間を探索しなくともある程度絞り込むことができます。例えば、初期化タイミングが 1 分オーダーでも分かっているなら、 60000 (ms) 個分試せばどれかは当たるでしょう。

過去の状態・出力列の逆算

Next() の逆演算 (直前の状態に戻す演算) を考えてみましょう。これが分かると、ある時点での内部状態が判明した際、過去の状態 (乱数列) の逆算ができます。難しいことは何もなく、減算を加算にしてインデックスを戻す処理を書くだけです:

uint prev = (uint)_seedArray[_inext] + (uint)_seedArray[_inextp];
if (prev >= int.MaxValue)
    prev -= int.MaxValue;

_seedArray[_inext] = (int)prev;

if (--_inext <= 0) _inext += 55;
if (--_inextp <= 0) _inextp += 55;

均等分布次元

System.Random の「均等分布次元」について考えてみましょう。よく Mersenne Twister の説明で「623 次元に均等分布する」と言われているアレです。

前提: 均等分布次元とは

そもそも均等分布次元とは何かというと、雑に説明すると「要素数 n の配列を連続する乱数で埋めたとき、全てのパターンが出現しうる最大の n」です。 623 次元に均等分布する Mersenne Twister では、 int[623] の配列を連続する Next() で埋めるとき、 int[623] で表現できるありとあらゆるパターンが出現しうることが保証されています。

この性質は当たり前に存在するわけではありません。例えば線形合同法は 1 次元均等分布なので、 int[2] ですらごく限られたパターンしか生成できません。より具体的には、以下の線形合同法では [0] = 0 なら [1] = 2531011 以外にはあり得ません。それ以外の 0 から始まるベクトル ({0, 1} など) は絶対に生成されません。線形合同法でランダムな 3 次元の点をプロットすると平面上に並ぶ、というのは有名な話です。

uint lcg(uint x) => x * 214013 + 2531011;

System.Random の均等分布次元

話を戻して、 System.Random ではどうでしょうか。そもそも遷移関数の実装ミスによって周期決定すらおぼつかない状態なので、実際の均等分布次元を求めることは困難です。はい……

ただ、もしも $\mathrm{GF}(2 ^ {31} - 1)$ 上で $x ^ {55} + x ^ {34} + 1$ が原始多項式であり、最長周期を実現できていたなら 55 次元に均等分布します。 (各要素の値域は変わらず $0 \le x \lt 2 ^ {31} - 1$ であることに注意は必要です。)

証明といえるほど厳密ではないですが説明しましょう。最長周期 $(2 ^ {31} - 1) ^ {55} - 1$ を実現できているなら、内部状態は全 0 を除いて全ての表現可能な状態を取ることができます。 (表現可能な状態のパターン数は最長周期と同じ数になります。) なお、全 0 は擬似乱数生成器の構造的にどうしようもないので見なかったことにします。 Mersenne Twister の論文 *8 の k 次元均等分布の定義 (Definition 1.1.) においても "except for the all-zero combination that occurs once less often." と述べられており、例外と考えてよいと思います。

そして、要素数 55 の内部状態と連続する 55 個の出力は一対一に対応します (全単射) 。言い換えれば、内部状態→出力の逆となる出力→内部状態の変換を行えます。具体的に以下のコードを用いて、連続する 55 個の出力から対応する内部状態を求められます。

// note: output.Length == 55, 0 <= output[i] < int.MaxValue
static int[] OutputToState(int[] output)
{
    if (output.Length != 55)
        throw new ArgumentException(nameof(output));

    const int Lag = 34;     // 55 - inext の初期値(21)

    var state = new int[56];

    int index, laggedIndex;       // == inext, inextp

    for (index = 55, laggedIndex = index - Lag;
        index > Lag;
        index--, laggedIndex--)
    {
        state[index] = output[index - 1] + output[laggedIndex - 1];
        if (state[index] < 0)
            state[index] += int.MaxValue;
    }
    for (laggedIndex = 55;
        index > 0;
        index--, laggedIndex--)
    {
        state[index] = output[index - 1] + state[laggedIndex];
        if (state[index] < 0)
            state[index] += int.MaxValue;
    }
    return state;
}

内部状態がすべての可能な値を取ることと、内部状態と連続する出力は一対一対応することから、連続する 55 個の出力もすべての可能な値を取るといえます。したがって、55 次元に均等分布するといえます。

余談: 好きな「乱数」列を出力させる

さて、上の項で内部状態⇔出力列の計算が可能であることを示しました。ということは、欲しい出力列に対応する内部状態を設定することで、55 個までの範囲で好きな値を出力させることができます。

上の OutputToState() メソッドを利用して、文字列を仕込んでみました。実行例を以下に示します。

var rand = new Random(1);
var state = new int[] { 0, 292, 414, 406, 405, 358, 194, 456, 327, 414, 354, 401, 419, 347, 183, 304, 247, 303, 251, 236, 334, 307, 257, 330, 302, 300, 243, 162, 347, 226, 299, 239, 304, 316, 246, 151, 185, 150, 188, 219, 137, 220, 206, 160, 568, 602, 619, 558, 247, 530, 441, 537, 380, 440, 553, 443 };

// リフレクションで _seedArray, _inext, _inextp を設定する
SetInternalState(rand, state, 0, 21);

var randomBytes = new byte[128];
rand.NextBytes(randomBytes);

// randomBytes を ASCII とみなして文字列化
Console.WriteLine(new string(randomBytes.Select(i => (char)i).ToArray()));

// > output:
// > OEEO@aAe?IU?#This message was created by Andante.Shioi-chan kawaii#AsS?oE2uo

# で囲まれているところが仕込んだ文字列です。 (# も含めて計 55 文字ジャストです。) 前述した Prev() を併用すれば、一定時間後に吐き出すような仕込みも可能です。

パフォーマンス

速度

パフォーマンスといえばまずは速度です。とりあえずベンチマークを取ってみましょう。今回は BenchmarkDotNet を使用し、.NET 5 環境で測定しました。

BenchmarkDotNet=v0.12.1, OS=Windows 10.0.19042
Intel Core i7-7700HQ CPU 2.80GHz (Kaby Lake), 1 CPU, 8 logical and 4 physical cores
.NET Core SDK=5.0.101
  [Host]     : .NET Core 5.0.1 (CoreCLR 5.0.120.57516, CoreFX 5.0.120.57516), X64 RyuJIT
  DefaultJob : .NET Core 5.0.1 (CoreCLR 5.0.120.57516, CoreFX 5.0.120.57516), X64 RyuJIT

せっかくなので、比較対象として拙作 Seiran128 擬似乱数生成器を使います。PractRand 32 TB をパスする品質を持ち、けっこう高速です。また、new() では RNGCryptoServiceProvider ベースの初期化を行ったり、 NextInt(int max), NextInt(int min, int max) では偏りがなくなるように操作を行ったりなど、速度を多少犠牲にしてでも品質を維持するように実装しています (つまり、速度だけ見てチューンするなら改善の余地がある、ということです) 。 Seiran 側の実装はこちらです。

各メソッドについて、処理時間は以下のようになりました。± は標準偏差です。

Method	System.Random [ns]	Seiran [ns]	x Speed
`new()`	1741.159 ± 23.02	120.738 ± 1.079	14.421
`new(401)`	470.529 ± 6.842	8.185 ± 0.087	57.487
`Next()`	9.093 ± 0.113	1.330 ± 0.054	6.837
`Next(401)`	11.748 ± 0.169	3.272 ± 0.039	3.590
`Next(168, 401)`	11.473 ± 0.138	4.567 ± 0.063	2.512
`NextDouble()`	9.734 ± 0.125	2.376 ± 0.034	4.097
`NextBytes(Span<byte>[256])`	2254.480 ± 27.31	90.324 ± 1.177	24.960

雑にグラフにすると以下のようになりました。 (グラフガチ勢に刺されそうです…)

f:id:andantesoft:20210109131436p:plain

Next***() がつぶれているので拡大するとこのようになります。

f:id:andantesoft:20210109131439p:plain

上表を見ると、new Random() が予想以上に遅いです。シード指定のある new Random(401) の 3.7 倍も遅いところを見ると、前述したシードが重複しないようにする処理が思いのほか重いのでしょうか。暗号論的擬似乱数生成器による初期化より遥かに遅いとはいったい… .NET Framework 環境ではほぼ new Random(401) と同じになるとは思いますが、ビルド環境がなかったので試せていません。

Seiran と比べると、出力系メソッドでは 2 ～ 4 倍 (NextBytes は 25 倍ぐらい) 、new では 14 ～ 57 倍 (!) ほど遅いことになります。 Seiran.NextBytes が速いのは、 MemoryMarshal.Cast() を使って直接 ulong の値を書き込んでいる (一度に 8 byte 分書き込める) ためです。C 言語でいうと ((uint64_t *)bytes)[i] = Next(); の処理をしているイメージです。対して Random.NextBytes は 1 byte / 回なので、かなりの差が付きます。

メモリ

消費メモリの観点からも見てみましょう。ゲーム分野などではヒープアロケーションに敏感になっている場合があるので、知っておいて損はないでしょう。

System.Random は、.NET 5 環境において 1 インスタンスあたり 280 byte 消費します。これは Visual Studio の診断ツールで確認しました。加えて、一度でも new Random() を呼べば追加で 560 byte 消費します。これはちょうど 2 インスタンス分で、new Random() の項で述べた ThreadStatic なインスタンスと static なインスタンスによるものと思われます。

ちなみに、Seiran では 1 インスタンスあたり 32 byte で済みます。 (内部状態を配列で持つと 64 byte になります。) 本質的には 16 byte 分なのでワンチャン構造体にできるかもしれませんが、今回はそこまでしませんでした。

おわりに

調べれば調べるほど問題や変な性質が見つかるので、何か楽しくなって調べていたら 1 ヵ月が経っていました。誤差レベルのものもあれば、結構影響のあるものもあります。

「だから System.Random は使うな！！！！」といった強い言葉を使うつもりはないですが、ゲームの内部計算・科学シミュレーションといった繊細な用途に使用する場合は、問題が表出する可能性があるため注意が必要だと思います。特に、再現性が必要な場合は流石にやめておいたほうが良いかと思います (前言撤回) 。問題が起きても後で交換できなくなりますので……

対して、品質を気にする必要がなく、かつ再現性が不要な場合なら雑に使ってもよいかとも思います。ゲームで言えばパーティクルの飛ばし方などでしょうか。いくら品質が低いといえども、全力で偏ったりはしないのである程度の使用なら大丈夫です。また、他では真似できない最強の利点が「標準ライブラリなのでいつでもどこでも簡単に使えること」です。適切に使っていきましょう。

さて、 dis るだけ dis って何もしないのもアレなので、ほぼ上位互換の代替実装を作ってみました。 (NextInt(0) が例外を投げるなど、厳密に互換性はないです。) MIT License で公開します。テストは結構しましたが、何かバグがあったらごめんなさい。

Seiran128 sample implementation in C# (.NET 5)

最後の最後に免責事項です。それなりに調査したつもりではありますが、私は数学徒ではないので数学理論周りに誤りやずれがある可能性はあります。もし何かありましたら優しく指摘していただけると大変助かります。よろしくお願いします。

*1:全 0 などの周期が短い状態も含むため、実用的にはこれより小さくなります

*2:ラグの指定ミスは流石に意図したものではない、と信じているので…

*3: Brent, R. P. (1994). On the periods of generalized Fibonacci recurrences. Mathematics of Computation, 63(207), 389-401.

*4:Aviv Sinai. (2011) Pseudo Random Number Generators in Programming Languages. The Interdisciplinary Center, Herzlia Efi Arazi School of Computer Science.

*5:実際には Next(int minValue, int maxValue) があるためこのような操作はしないと思いますが

*6:厳密には「誤検知」や「テスト自身の妥当性」など若干難しいところはありますが、ここでは大目に見てください。

*7:Blackman, D., & Vigna, S. (2018). Scrambled linear pseudorandom number generators. arXiv preprint arXiv:1805.01407.

*8:Matsumoto, M., & Nishimura, T. (1998). Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Transactions on Modeling and Computer Simulation (TOMACS), 8(1), 3-30.