乱数生成の備忘録２ - エフアンダーバー

他人が書いた乱数生成のコードを読んでいたところ、前回自分が書いたものと違う実装をしていたので改めていろいろと調査。

問題のコード

異なる実装となっていたのは整数の乱数から[0,1)の範囲の実数の乱数を生成するコード。

前回自分が書いたコードでは次のように定数を掛けて変換していた。

static Double ToDouble(UInt64 x)
{
    return (x * ((1.0 - 1.0 / (1ul << 52)) / UInt64.MaxValue));
}

今回見つけたコードは次のような実装となっていた*1。

static Double ToDouble(UInt64 x)
{
    const UInt64 ExponentBits = 0x3ff0_0000_0000_0000;
    const UInt64 MantissaMask = 0x000f_ffff_ffff_ffff;
    UInt64 bits = (x & MantissaMask) | ExponentBits;
    Double v = BitConverter.Int64BitsToDouble(unchecked((Int64)bits));
    return (v - 1.0);
}

double型（倍精度浮動小数点数）はIEEE754規格に従っているので、 1bitの符号部と11bitの指数部および52bitの仮数部から成る。

この実装では符号部を正(0)とし、指数部が0(3ff)となるようにして、残りの仮数部を乱数の値で埋めている。このようにすると、数値は単純に仮数部の先頭に1を加えた1.XXX...XXXで表される[1,2)の範囲の値となる。最後にその値から1を引けば値は[0,1)の範囲に収まることになる。

メジャーなプロジェクトとしてV8のコードを見てみたところ、こちらも後者で実装されていた。

何が違うか

前者のコードでは一様分布な整数の乱数の範囲を[0,1)に調整した後で値が丸められる。異なる実数であっても、丸めによって同じ値となってしまうことがあるため、結果の浮動小数点数自体の分布は一様ではなくなる。浮動小数点数が表現できる値には偏りがあり、0に近いほど高い精度で表現できるため、比較的精度の低い1付近では0付近に比べて丸めによる同じ値の発生が起きやすい。ただし、表現している実数自体は一様に分布している。

一方、後者のコードは計算機イプシロンの間隔で並ぶ[0,1)の範囲の浮動小数点数の部分集合のみを利用して、それらに一様の分布を与えている。つまり、[0,1)の範囲で最も精度の悪い部分に合わせることにより、浮動小数点数（の部分集合）自体の分布をも一様としている。欠点は[0,1)の範囲を表現する浮動小数点数62bit分の内の52bit分(1/1024)しか使用できないこと。しかし、削られた数値からの距離が計算機イプシロンの半分以下である値が必ず含まれるので気にするほどではないかもしれない。

どちらがいいか

いくらか調べてはみたもののどちらかが決定的にダメだと思えるような例は見当たらなかった。ただし、後者を強く支持する人がいる一方で、前者はこれも悪いとは言えないよねくらいの感じだったので後者の方が無難かもしれない。

しかし、後者は実装する言語によっては少し面倒くさい。例えば、C#ではBitConverter.Int64BitsToDoubleに対応するfloat(単精度浮動小数点数）用のメソッドが見当たらない。そのため次のように少し回りくどいことをしなければならない。

byte配列を介してBitConverterで変換
- ヒープメモリの確保が必要
- オブジェクトを使いまわすならスレッドセーフでなくなる（あるいはThreadStaticAttributeが必要）
unsafeでreinterpret_cast的な変換
StructLayoutAttributeで共用体的な変換

・・・とか考えていたところ、最初に見つけたコード(JavaScript)が次のようなコードに更新されていた（実際には定数部分は最適化されている）。

function toDouble(upper, lower) {
    return upper * Math.pow(2, -32) + (lower >> 12) * Math.pow(2, -52);
}

C#で書くならこんな感じ。

static Double ToDouble(UInt64 x)
{
    return (x >> 12) * (1.0 / (1ul << 52));
}

これならfloat(Single)でも同じように書ける。

static Single ToSingle(UInt32 x)
{
    return (x >> 9) * (1.0f / (1u << 23));
}

乗算が含まれるとはいえ、回りくどい変換をするよりはこちらの方が速度も良さそうな気がする。

参考

おわりに

調べてみて改めて思ったのは乱数とか浮動小数点数とかすごくメンドクサイな、と。正しいか正しくないかの検証が楽ならいいのだけど、どうなれば正しいのかがよくわからないのでいまいち結果に自信を持てない。

書いている途中で気づいたのだけど、前回の記事でdouble(64bit)の乱数を生成するのにuint(32bit)の乱数を元にしているという。分布がどうのこうの言うより前にまともなコードを書かねば・・・。

*1:下記は同様の実装となるようにC#で書き換えたもの