メニエスの名のもとに

プログラミング関係を中心としたぐだぐだブログ

一様乱数, それとは別に単精度浮動小数点数乱数

コメントに答えるには、コメント欄が狭いかなと思ったが、書いてみたらそれほど長くもなかったような。

一様乱数

C++ のstd::uniform_int_distributionが遅いという話。ある範囲の数を一様に出力する乱数が欲しい。これは簡単なようでいてそれほど簡単ではない。
メルセンヌ・ツイスタのページのQandAコーナーよく聞かれる質問

サイコロみたいに、[0..N-1]の一様乱数を生成したい
N があまり大きくなくて(だいたい < 2^16) かつ、あまり精度を必要としなければ、 [0,1)区間版(倍精度)を使用して、出力をN倍して切り捨てれば十分です。 (注意:必ず半開区間バージョン[0,1) を使うこと。閉区間バージョンでこれをやると、約42億回に1回の割合で 1.0(16進整数で言うと0xFFFFFFFF) が出現し、値として N が返されるというたちの悪いバグが起こります。 和田維作さんの指摘。2002/Jan.8th added)
WARNING: この方法は、Nが2の累乗でない限り、分布に誤差が生じます。 誤差のオーダーはN×2^-32になります。 より優れた、より安全な方法は、 整数乱数を生成し、N<=2nとなる最小のnをとって、二進で上位n桁 を見て、N以上の乱数は捨てることです。

関係ないけど、上付き文字を指定しても下付き文字に見えるんだよな。このページで下付き文字に見えるとしても上付き文字だと思ってください。
C++ のstd::uniform_int_distributionはここでいうところの「より優れた、より安全な方法」を使っていると思われる。
この生成法ならば、(可能な限り)正確な一様乱数を生成することができるが、ifによる条件判断が入る。それどころか、運が悪ければ、1個生成するのに、複数回条件判断で乱数を捨てなければならない。
しかし、サイコロ(N=6)のように疑似乱数生成器の出力ビット数に対して、Nが十分に小さければ、誤差は小さいので最初の方法で生成しても良いかもしれない。そこはユーザーの判断となるだろう。つまり、Nが小さく速度が欲しければstd::uniform_int_distributionは使わないという選択もあるということである。
また上記のQandAで「誤差のオーダーはN×2^-32」と言っているのは、32ビットMTから倍精度に変換したものから作っているから、53ビット精度のdouble を使用すれば、誤差のオーダーはN×2^-53になる。
分布の誤差を小さくする方法としては、64ビット整数から80ビット拡張倍精度浮動小数点数にして、簡易法を使えばN×2^-64の誤差に出来る。あるいは、4倍精度を使う。
どれが速いかは、例によって全部やって比較しないと分からない。

単精度、半精度浮動小数点数

32bit整数から作ればよいのだが、dSFMTのように浮動小数点数形式で作ればより速いと考えられる。しかし、ここで一つ問題がある。それは擬似乱数生成器の統計的検定をするために、TestU01というパッケージを使っていることである。いや、このパッケージの中のBigCrushというテストが非常に多くの統計的検定をやってくれるので大変便利なのである。だが、BigCrushは基本的に32ビット用である。64bitの疑似乱数を検定する時は、上位32ビットだけの出力、下位32ビットだけの出力、上位下位交互の出力というように分けて検定している。
しかし、単精度をダイレクトに作ると、23ビットの生成、半精度だと10ビットの生成である。これをそのままBigCrushに掛けると、まず検定を通らない。BigCrushを使わずに個々の統計的検定関数を呼べばよいのだが、パラメータが多くて相当に大変(だと予想される)。さらに、適当なパラメータを設定して検定を通ったとしても、そのパラメータが正しいのか、意味のある検定をしているのかという疑問が出てくる。まあ、論文を投稿すれば査読者に指摘されるだろう。
(実をいえば、BigCrushで使っているパラメータは本当に適切なのかという疑問もある。だが、エライ教授が監修して作ったパッケージでみんなが使っているので、信頼されている)
32ビット整数から作った場合に比べて本当に速いのかというのも、作ってみないと分からない。
この辺は、学部生なり院生なりが卒論、修論として取り組むと面白いかも知れない。