C++とC#で各種乱数発生器を実装してみる（勉強用）

はじめに

筆者の先の記事

において，C++ の各種乱数発生器を原理も分からずにとりあえず使ってみましたが，疑似乱数の勉強がてら，それらと同等なものを C# と C++ で実装してみることにしました。

お品書き

今回挑戦する乱数発生器は以下の６種（７個）です。生成する乱数の型は C# のものです。C++ の場合は int → int32_t, uint → uint32_t，ulong → uint64_t と読み換えて下さい。

アルゴリズム	型	範囲	備考
線形合同法	uint	$1 \sim 2^{31} - 2$	C++ の std::minstd_rand
減算乱数ジェネレータ	int	$0 \sim 2^{31} - 2$	.NET Framework の Random
リオーダーアルゴリズム	uint	$1 \sim 2^{31} - 2$	C++ の std:knuth_b
メルセンヌツイスター	uint	$0 \sim 2^{32} - 1$	C++ の std:mt19937
メルセンヌツイスター64bit版	ulong	$0 \sim 2^{64} - 1$	C++ の std:mt19937_64
XorShift 64bit版	ulong	$1 \sim 2^{64} - 1$

線形合同法

一般的に線形合同法とは，乱数の種 $S$ とおくと初期値 $x_0 = S$ として

x_{n + 1} = (A \cdot x_n + C) \bmod M

という漸化式で与えられます。特に $C = 0$ としたものを乗法合同法ともいい，C++ の標準ライブラリに採用されています。std::minstd_rand0 は $A = 16807$，std::minstd_rand は $A = 48271$ とし，いずれも$M = 2^{31} - 1 = 2147483647$ としています。

なぜ $M = 2^{31} - 1$ なのかというと，$2^{31} - 1$ が符号付き 32bit 整数で表せる最大の数であり，素数だからと思われます。$C = 0$ かつ $M$ が素数であれば $0 < x_n < M$ のとき必ず $0 < x_{n + 1} < M$ になることを簡単に証明できるからです。

素数 $M$ より小さい自然数 $a$，$b$ の積 $a \times b$ は $M$ では割り切れない。

ただし $M$ は符号付き 32bit 整数で表せる最大の数なので $A \cdot x_n + C$ は符号付き 32bit 整数の範囲を超え得ることに注意する必要があります。もしもオーバーフローして下位 32bit のみを用いた場合，それは $2^{32}$ での剰余を求めたことになるからです。すなわちこれらの演算を正しく行うためには通常 64bit ÷ 32bit の除算（剰余演算）を用いるため，ソースコード上はとてもシンプルに見えて実は少しだけ重い処理になっています。

線形合同法（C# 版）

class LinearCongruential {
	private	uint	x = 1;
	public	LinearCongruential(){}
	public	LinearCongruential(uint s) {
		if(s < 1 || s >= 2147483647) s = 1;
		x = s;
	}
	public	uint	Next() {
		x = (uint)(((ulong)x * 48271) % 2147483647);
		return x;
	}
}

線形合同法（C++ 版）

class LinearCongruential {
private:
	uint32_t	x = 1;
public:
	LinearCongruential(){}
	LinearCongruential(uint32_t s) {
		if(s < 1 || s >= 2147483647) s = 1;
		x = s;
	}
	uint32_t	operator()() {
		x = (uint32_t)(((uint64_t)x * 48271) % 2147483647);
		return x;
	}
};

＜参考文献＞

• linear_congruential_engine クラス - microsoft.com

減算乱数ジェネレーターアルゴリズム（Knuth-A）

.NET の Random クラスに採用されているもので，D.E.Knuth の著書「The Art of Computer Programming Volume 2 Seminumerical Algorithms」では「加算式乱数生成器」アルゴリズム A として紹介されています。ただし，Knuth-A と呼ぶ人はあまり聞いたことがないですね。

Knuth 本では二つの乱数の「加算」を行っていたものが Micorosft の実装では「減算」になっていたり，55 個のバッファの掃引方向が逆になっていたりなど実装の細部は異なっています。Microsoft 自身も

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

と述べています。もっとも Microsoft の実装は必要な 55 個のバッファに対して配列を 56 個で宣言し，かつ配列の先頭要素（0 番目）を使用していないので，おそらくは 1 オリジンの言語（FORTRANとか？）からコードを流用したきたものではないかと推察しています。以下は，確保する配列を 55 個で済むように筆者が作り替えたものです。

なお，個々の乱数生成時の演算には加減算しか使用しないことから Knuth 本では高速アルゴリズムとして紹介されています。乗算が高コストだった時代には正しかったのかもしれませんが，現在においてはむしろ乗算を用いる線形合同法よりも遅いのは悲しいところです。

減算乱数ジェネレーターアルゴリズム（C# 版）

class Knuth_A {
	private	int[]	a = new int[55];
	private	int		p = 0;
	private	int		q = 21;
	public	Knuth_A() : this(0) {}
	public	Knuth_A(int s) {
		/**/ if(s == -2147483648) s = 2147483647;
		else if(s < 0) s = -s;
		int	j = 161803398 - s;
		int	k = 1;
		a[54] = j;
		for(int i = 0; i < 54; i++) {
			int	n = k;
			a[(21 * i + 20) % 55] = n;
			if((k = j - k) < 0) k += 2147483647;
			j = n;
		}
		for(int h = 0; h < 4; h++) {
			for(int i = 0; i < 55; i++) {
				int	n = a[i] - a[(i + 31) % 55];
				if(n < 0) n += 2147483647;
				a[i] = n;
			}
		}
	}
	public	int	Next() {
		int	x = a[p] - a[q];
		/**/ if(x == 2147483647) x--;
		else if(x < 0) x += 2147483647;
		a[p] = x;
		if(++p >= 55) p = 0;
		if(++q >= 55) q = 0;
		return x;
	}
}

減算乱数ジェネレーターアルゴリズム（C++ 版）

class Knuth_A {
private:
	int	a[55];
	int	p = 0;
	int	q = 21;
public:
	Knuth_A() : Knuth_A(0) {}
	Knuth_A(int s) {
		/**/ if(s == -2147483648) s = 2147483647;
		else if(s < 0) s = -s;
		int	j = 161803398 - s;
		int	k = 1;
		a[54] = j;
		for(int i = 0; i < 54; i++) {
			int	n = k;
			a[(21 * i + 20) % 55] = n;
			if((k = j - k) < 0) k += 2147483647;
			j = n;
		}
		for(int h = 0; h < 4; h++) {
			for(int i = 0; i < 55; i++) {
				int	n = a[i] - a[(i + 31) % 55];
				if(n < 0) n += 2147483647;
				a[i] = n;
			}
		}
	}
	int	operator()() {
		int	x = a[p] - a[q];
		/**/ if(x == 2147483647) x--;
		else if(x < 0) x += 2147483647;
		a[p] = x;
		if(++p >= 55) p = 0;
		if(++q >= 55) q = 0;
		return x;
	}
};

＜参考文献＞

• Random クラス - microsoft.com
• The Art of Computer Programming Volume 2 Seminumerical Algorithms Third Edition 日本語版
• .NET System.Random の実装と欠陥について ～ 重箱の隅をつつきたおす ～ 屋根裏工房改 - HatenaBlog

リオーダーアルゴリズム（Knuth-B）

線形合同法である std::minstd_rand0 を内部で利用しています。256 個のリオーダーバッファを持ち，線形合同法によって得られた乱数の順番を入れ替えているだけにも関わらず，乱数の品質は向上するという不思議な乱数発生器です。D.E.Knuth の著書「The Art of Computer Programming Volume 2 Seminumerical Algorithms」にアルゴリズム B として紹介されていることから通称 Knuth-B とも呼ばれているようです。

position() は生成した乱数の値を用いてリオーダーバッファの位置を求めていますが，その際，乱数が [1, 2147483646] の範囲で一様分布することを想定しています。乱数の値が 0 にならないことは線形合同法の章で説明した通りです。

線形合同法を内部で呼び出していることからスピードは線形合同法より遅くなることは確実ですが，C++ の標準ライブラリのように 2～3 倍も遅くなるというのはちょっと信じられないですね。筆者が実装した C# 版では線形合同法に対してほんの僅か劣る程度ですから，C++ の標準ライブラリの実装があまり良くないのかもしれません。

リオーダーアルゴリズム（C# 版）

class Knuth_B {
	private	uint[]	a = new uint[256];
	private	uint	x;
	private	int		p;
	private	uint	update(uint y) {
		return (uint)(((ulong)16807 * y) % 2147483647);
	}
	private	uint	position(uint y) {
		return (uint)((ulong)256 * (y - 1) / (2147483647 - 1));
	}
	public	Knuth_B() : this(1) {}
	public	Knuth_B(uint s){
		if(s < 1 || s >= 2147483647) s = 1;
		x = s;
		for(int i = 0; i < a.Length; i++)
			a[i] = x = update(x);
		x = update(x);
		p = position(x);
	}
	public	uint	Next() {
		uint	y = a[p];
		a[p] = x = update(x);
		p = position(y);
		return y;
	}
}

リオーダーアルゴリズム（C++ 版）

class Knuth_B {
private:
	uint32_t	a[256];
	uint32_t	x;
	uint32_t	p;
	uint32_t	update(uint32_t y) {
		return (uint32_t)(((uint64_t)16807 * y) % 2147483647);
	}
	uint32_t	position(uint32_t y) {
		return (uint32_t)((uint64_t)256 * (y - 1) / (2147483647 - 1));
	}
public:
	Knuth_B() : Knuth_B(1) {}
	Knuth_B(uint32_t s) {
		if(s < 1 || s >= 2147483647) s = 1;
		x = s;
		for(int i = 0; i < 256; i++)
			a[i] = x = update(x);
		x = update(x);
		p = position(x);
	}
	uint32_t	operator()() {
		uint32_t	y = a[p];
		a[p] = x = update(x);
		p = position(y);
		return y;
	}
};

＜参考文献＞

• shuffle_order_engine クラス - microsoft.com
• The Art of Computer Programming Volume 2 Seminumerical Algorithms Third Edition 日本語版

メルセンヌツイスター（MT19937）

C++ のデフォルト乱数 default_random_engine にもなっている日本製の乱数発生器です。理論的にはさっぱり理解できていませんが，個々の乱数生成時の演算にはビットシフトと論理演算のみ使用しており，おそらくこれが高速性の一因でしょう。なお，筆者の勝手な想像ですが，後発のより高速な乱数発生器である XorShift にも影響を与えたのではないのでしょうか？

なお，C# 版では x の偶奇に応じて加算される値を配列 mag[] として実装していますが，C++ 版では三項演算子で実装しています。これは，それぞれその方が速かったからです。C++ のコンパイラ（clang-cl）では三項演算子が cmov 命令のような効率の良い命令に変換されるのに対し，C# のコンパイラはまだそこまで最適化できていないからだと考えています。

メルセンヌツイスター（C# 版）

class MersenneTwister {
	const	uint	N   = 624;
	const	uint	M   = 397;
	const	uint	MSK = 0x7FFFFFFF;
	private	uint[]	mag = new uint[]{ 0, 0x9908B0DF };
	private	uint[]	mt  = new uint[N];
	private	uint	mti;
	public MersenneTwister() : this(5489) {}
	public MersenneTwister(uint s) {
		mt[0] = s;
		for(mti = 1; mti < N; mti++)
			mt[mti] = 1812433253 * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti;
	}
	public	uint	Next() {
		uint	x;
		if(mti >= N) {
			uint	i, j;
			for(i = 0, j = M; i < N - M; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ mag[x & 1];
			}
			for(j = 0; j < M - 1; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ mag[x & 1];
			}
			x = (mt[i] & ~MSK) | (mt[0] & MSK);
			mt[i] = mt[j] ^ (x >> 1) ^ mag[x & 1];
			mti = 0;
		}
		x = mt[mti++];
		x ^= (x >> 11);
		x ^= (x <<  7) & 0x9D2C5680;
		x ^= (x << 15) & 0xEFC60000;
		x ^= (x >> 18);
		return x;
	}
}

メルセンヌツイスター（C++ 版）

class MersenneTwister {
	static	const	uint32_t	N   = 624;
	static	const	uint32_t	M   = 397;
	static	const	uint32_t	MAG = 0x9908B0DF;
	static	const	uint32_t	MSK = 0x7FFFFFFF;
private:
	uint32_t	mt[N];
	uint32_t	mti;
public:
	MersenneTwister() : MersenneTwister(5489) {}
	MersenneTwister(uint32_t s) {
		mt[0] = s;
		for(mti = 1; mti < N; mti++)
			mt[mti] = 1812433253 * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti;
	}
	uint32_t	operator()() {
		uint32_t	x = 0;
		if(mti >= N) {
			uint32_t	i, j;
			for(i = 0, j = M; i < N - M; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ ((x & 1) == 0 ? 0 : MAG);
			}
			for(j = 0; j < M - 1; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ ((x & 1) == 0 ? 0 : MAG);
			}
			x = (mt[i] & ~MSK) | (mt[0] & MSK);
			mt[i] = mt[j] ^ (x >> 1) ^ ((x & 1) == 0 ? 0 : MAG);
			mti = 0;
		}
		x = mt[mti++];
		x ^= (x >> 11);
		x ^= (x <<  7) & 0x9D2C5680;
		x ^= (x << 15) & 0xEFC60000;
		x ^= (x >> 18);
		return x;
	}
};

＜参考文献＞

• 2002年版(mt19937ar.c)の 説明・Ｃコード集 - Mersenne Twister Home Page

• A 623-dimensionally equidistributed uniform pseudorandom number generator - ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, January pp.3-30 (1998)

メルセンヌツイスター64bit版（MT19937_64）

メルセンヌツイスターの 64bit 版です。筆者はこれを 32bit 版とまったく同じ乱数列を生成するもので，かつ 64bit プロッサ向けに高速なアルゴリズムを適用したものだと勘違いしていました。正しくは 64bit 版のメルセンヌツイスターは 64bit 整数の乱数列を生成するものです。なお，高速版は別途用意されており，SFMT（SIMD-oriented Fast Mersenne Twister）というものです。

なお，メルセンヌツイスター 32bit 版とコードを比較すれば分かりますが演算量自体はほぼ変わりません。配列 mt[] も 64bit 幅に増えた一方，要素数 N は半減しているからです。これは 623 次元の乱数生成に必要なテーブル個数が半分になったということでしょう。一見，計算量は半減したように思えますが更新周期も半分になったので，長期的に稼働させれば計算量は変わらないはずです。

メルセンヌツイスター64bit版（C# 版）

class MersenneTwister64 {
	const	uint	N   = 312;
	const	uint	M   = 156;
	const	ulong	MSK = 0x000000007FFFFFFFL;
	private	ulong[]	mag = new ulong[]{ 0, 0xB5026F5AA96619E9L };
	private	ulong[]	mt  = new ulong[N];
	private	uint	mti;
	public MersenneTwister64() : this(5489) {}
	public MersenneTwister64(ulong s) {
		mt[0] = s;
		for(mti = 1; mti < N; mti++)
			mt[mti] = 6364136223846793005L * (mt[mti - 1] ^ (mt[mti - 1] >> 62)) + mti;
	}
	public	ulong	Next() {
		ulong	x;
		if(mti >= N) {
			uint	i, j;
			for(i = 0, j = M; i < N - M; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ mag[x & 1];
			}
			for(j = 0; j < M - 1; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ mag[x & 1];
			}
			x = (mt[i] & ~MSK) | (mt[0] & MSK);
			mt[i] = mt[j] ^ (x >> 1) ^ mag[x & 1];
			mti = 0;
		}
		x = mt[mti++];
		x ^= (x >> 29) & 0x5555555555555555L;
		x ^= (x << 17) & 0x71D67FFFEDA60000L;
		x ^= (x << 37) & 0xFFF7EEE000000000L;
		x ^= (x >> 43);
		return x;
	}
}

メルセンヌツイスター64bit版（C++ 版）

class MersenneTwister64 {
	static	const	uint32_t	N   = 312;
	static	const	uint32_t	M   = 156;
	static	const	uint64_t	MAG = 0xB5026F5AA96619E9LL;
	static	const	uint64_t	MSK = 0x000000007FFFFFFFLL;
private:
	uint64_t	mt[N];
	uint32_t	mti;
public:
	MersenneTwister64() : MersenneTwister64(5489) {}
	MersenneTwister64(uint64_t s) {
		mt[0] = s;
		for(mti = 1; mti < N; mti++)
			mt[mti] = 6364136223846793005LL * (mt[mti - 1] ^ (mt[mti - 1] >> 62)) + mti;
	}
	uint64_t	operator()() {
		uint64_t	x;
		if(mti >= N) {
			uint32_t	i, j;
			for(i = 0, j = M; i < N - M; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ ((x & 1) == 0 ? 0 : MAG);
			}
			for(j = 0; j < M - 1; i++, j++) {
				x = (mt[i] & ~MSK) | (mt[i + 1] & MSK);
				mt[i] = mt[j] ^ (x >> 1) ^ ((x & 1) == 0 ? 0 : MAG);
			}
			x = (mt[i] & ~MSK) | (mt[0] & MSK);
			mt[i] = mt[j] ^ (x >> 1) ^ ((x & 1) == 0 ? 0 : MAG);
			mti = 0;
		}
		x = mt[mti++];
		x ^= (x >> 29) & 0x5555555555555555LL;
		x ^= (x << 17) & 0x71D67FFFEDA60000LL;
		x ^= (x << 37) & 0xFFF7EEE000000000LL;
		x ^= (x >> 43);
		return x;
	}
};

＜参考文献＞

• 64bit版 Mersenne Twister - Mersenne Twister Home Page

XorShift

2003年に発表されたアルゴリズムで何と言っても実装のシンプルさと高速性が売りでしょう。なお，先ほどメルセンヌツイスターが後発の XorShift に影響を与えたかもしれないと言いましたが，違うような気がしてきました。考案者である Marsaglia 氏の論文を読んでもメルセンヌツイスターの論文を引用していないからです。というか Marsaglia 氏は自著論文しか引用していないのです・・・

XorShift 64bit版（C# 版）

public class XorShift {
	private	ulong	x = 1;
	public	XorShift() {}
	public	XorShift(ulong s) {x = s;}
	public	ulong	Next() {
		x ^= x << 13;
		x ^= x >>  7;
		x ^= x << 17;
		return x;
	}
}

XorShift 64bit版（C++ 版）

class XorShift {
private:
	uint64_t	x = 1;
public:
	XorShift() {}
	XorShift(uint64_t s) {x = s;}
	uint64_t	operator()() {
		x ^= x << 13;
		x ^= x >>  7;
		x ^= x << 17;
		return x;
	}
};

＜参考文献＞

• XorShift - wikipedia
• Xorshift RNGs - Journal of Statistical Software

パフォーマンス比較

まず独自実装した乱数発生器は，.NET Framework のものや C++ 標準ライブラリと比較して，同じ種からスタートすれば同じ乱数列が得られることを確認しました。
※全数検査は無理なのでサンプル検査した範囲に限ります。

テスト題材は拙作の記事「プリキュアカレーのシール全30種コンプするには何箱必要か？」です。計算時間を表１に示します。計測環境ですが，プロセッサは Core-i7 13700，OS は Windows 10 22H2（64bit 版）です。C# コンパイラは Visual Studio 2022 Commnunity Eition のもの，C++ コンパイラは clang-cl 17.0.3（64bit 版）を使いました。

表１　計算時間（単位は秒）

乱数生成器	C# 版		C++ 版		備考
	独自	.NET	独自	標準
線形合同法	397.94		296.40	296.44	C++ の std::minstd_rand
Knuth_A	667.00	717.19	226.29		.NET Framework の Random
Knuth_B	420.32		284.93	645.35	C++ の std::knuth_b
MT19937	372.78		251.84	203.28	C++ の std::mt19937
MT19937_64	409.29		286.33	203.23	C++ の std::mt19937_64
XorShift 64bit版	194.36		180.75

線形合同法

最適化の余地があまりないせいか，C++ の独自実装版と標準ライブラリで差異が殆どありません。一方，C# は三割ほどパフォーマンスが落ちます。整数除算が重いはずなので，C++ のコンパイラはこの演算遅延をうまく隠蔽できているのではないかと思います。

Knuth_A

C# の独自実装版は .NET Framework の Random クラスよりも僅かに速くなりましたが，C++ の独自実装版はさらに三倍近く速くなっています。テーブル参照（メモリアクセス）が足を引っ張るはずですが，C++ のコンパイラは最適化によりメモリレイテンシをうまく隠蔽できているのではないかと思います。

Knuth_B

唯一，C++ の標準ライブラリのパフォーマンスが奮いませんでした。なぜこんなに遅いのか原因は分かりません。

MT19937 と MT19937_64

独自実装版は C++ 版も C# 版のいずれも MT19937_64 のほうが MT19937 よりも一割ほど遅いという結果になりましたが，何故か C++ の標準ライブラリはほぼ同じという結果となりました。メルセンヌ・ツイスターでは内部で大量のデータを必要としますが，MT19937 は 32bit 整数を 624 個使い，MT19937_64 も 64bit 整数を 312 個使うので必要なデータの総量（2,496 bytes）も同じです。一次データキャッシュに入り切りそうな絶妙なサイズ感ですね。ただし，MT19937 では乱数を 624 個生成する毎に 624 個（2,496 bytes）のデータを更新しますが，MT19937_64 では 312 個毎に 312 個（2,496 bytes）のデータを更新します。すなわち乱数一個あたり更新する平均データサイズは MT19937_64 のほうが二倍になりますが，平均データ個数は同じ一個ずつで変わりません。そもそも 64bit 演算（およびメモリアクセス）を 32bit 演算と同じ速度で行えたら，MT19937_64 も MT19937 と同じスピードで実行できるはずで，実際 C++ の標準ライブラリでは等速を実現できています。しかし，独自実装版ではそのレベルに達していません。

XorShift

最速なのは変わらないですね。

結論

• C++ ユーザーの場合，スピードと乱数の品質のバランスから考えて標準ライブラリの MT19937 または MT19937_64 をそのまま使えば良いと思います。

• C# ユーザの場合，標準ライブラリの Random クラスの出来（スピードおよび乱数の品質）があまり良くないので，乱数の品質を要求される方は MT19937 または M19937_64，スピードおよび実装の簡便さを求める方は XorShift を使いましょう。筆者のように独自実装に拘らなくとも，どこかに落ちているでしょう。

• C# コンパイラの最適化性能はまだまだ C++ コンパイラに及ばないようです。

備考

DOTNET のソースコードが置いてあるサイトを覗くと，XorShift の進化版である Xoshiro が用意されているように見えます。

• Xoshiro128 - github.com/dotnet/runtime

• Xoshiro256 - github.com/dotnet/runtime