Upgrade to Pro — share decks privately, control downloads, hide ads and more …

[0.0, 1.0) の乱数を得るための“本当の”方法

hole
September 03, 2023

[0.0, 1.0) の乱数を得るための“本当の”方法

レイトレ合宿9(*)のセミナー発表スライドです。
* https://sites.google.com/view/rtcamp9/home

- 2023/09/08 “除算法2”追記。(@Reputelessさんありがとうございました)

hole

September 03, 2023
Tweet

More Decks by hole

Other Decks in Programming

Transcript

  1. 考察する手法 • 今回は以下の手法について考察する。 • std::generate_canonical • 除算法 • mantissa(仮数部)法 •

    除算法2 • exponent+mantissa法 • Abseilの最適化 • Downeyの方法 • 名称は便宜的に付けたものも含まれる。
  2. std::generate_canonical • まずは標準C++ライブラリから始める。標準C++ライブラリには[0.0, 1.0) 区間で一様乱数を得る関数が用意されている。 • std::uniform_real_distribution の実装に使われていること が多い。 実数区間[0.0,

    1.0)に展開(事実上正規化)された一様分布乱数を得る ための関数テンプレート。 テンプレート引数に与える RealType 型について、 bits だけの分解能を持つ仮数部を URBG gを必要な回数だけ繰り返し呼 び出して生成する。 https://cpprefjp.github.io/reference/random/generate_canonical.html
  3. std::generate_canonical • 使い方 • std::generate_canonical<float, 23>(random_generator); • std::generate_canonical<double, 52>(random_generator); •

    返却する型(float,double等)と、仮数部の分解能の最低要求(float, doubleのフルスペックなら23bit, 52bit)を指定する。 • random_generator の生成ビット数と、分解能要求に応じて、複数 回の乱数生成が行われる。 • 32bit乱数生成器、52bit分解能の場合、二回乱数生成が行われる。 • 今回はfloatに限定しているので、32bit乱数生成器、23bit分解能 で良い。
  4. std::generate_canonicalの問題 • これを使えばいいのか? → ダメ! • 幾つかの問題が存在している。 • 問題1:一部のライブラリ実装において、[0.0, 1.0)が保証されない。

    • 問題2:一部のライブラリ実装において、非常に低速。 • 問題3:一部のライブラリ実装において、[0.0, 1.0)の範囲における、表現可 能な全てのfloat値を得ることが出来ない。
  5. 問題1:一部のライブラリ実装において、[0.0, 1.0) が保証されない • ライブラリ実装についてエッジケースを観察する。 • 0xffffffffを乱数として与える。1.0fが出力されなければOK。 • 実験方法 •

    返却型として、floatを使用。 • 乱数生成器として、32bitを使用。 • 以下のライブラリコードを抽出、共通のコンパイラ上でコンパイル・実行。ターゲット はx64。 • 対象ライブラリ • msvc new (>=19.32) • msvc old (<= 19.31) • gcc libstdc++ (2023/08/22 時点でのrepo最新) • clang libc++ (2023/08/22 時点でのrepo最新) • Abseil (2023/08/22 時点でのrepo最新) • Abseil については同様の機能を持つ別の関数、GenerateRealFromBitsを使用。
  6. 結果 • msvc new • 32bit rng/float: 1.0000000000000000e+00 • msvc

    old • 32bit rng/float: 1.0000000000000000e+00 • gcc libstdc++ • 32bit rng/float: 9.9999994039535522e-01 • clang libc++ • 32bit rng/float: 1.0000000000000000e+00 • Abseil • 64bit rng/float: 9.9999994039535522e-01
  7. • msvc new • 32bit rng/float: 1.0000000000000000e+00 • msvc old

    • 32bit rng/float: 1.0000000000000000e+00 • gcc libstdc++ • 32bit rng/float: 9.9999994039535522e-01 • clang libc++ • 32bit rng/float: 1.0000000000000000e+00 • Abseil • 64bit rng/float: 9.9999994039535522e-01 結果 多くのライブラリで、1.0fを出力してしまう!
  8. 原因:generate_canonicalの実装 • いずれのライブラリにおいても、本質的には以下の様になっている(Abseil 除く) • (return_type)random() / ((return_type)(RandomMax)+(return_type)1) • たとえばreturn_type=float、32bit乱数生成器の場合、丸め誤差に

    より以下のような計算結果になる! • RandomMax = 0xffffffff = 4294967295 • 分子: (float)4294967295 → 4294967296.0f • 分母: ((float)(4294967295) + (float)1) → 4294967296.0f • よって、random()の結果がRandomMaxだと最終結果が1.0fとなる。
  9. 原因:generate_canonicalの実装 • gccにおける対策:next_after • if (__ret >= _RealType(1)) • {

    • __ret = std::nextafter(_RealType(1), _RealType(0)); • } • 1より小さい、最大の浮動小数値を返す。 • 実用上は問題にならないと考えられるが、該当の値の生成確率がバイアスさ れる。 指定方向への次の表現可能な値を取得する。この関数は、パラメータxの値を パラメータyの値の方向に対して、その環境で表現可能な最小の値だけ進める。 https://cpprefjp.github.io/reference/cmath/nextafter.html
  10. 問題2:一部のライブラリ実装において、非常に低速 • msvc old (<= 19.31) における実装は非常に低速。 • 1億回の実行に要する平均時間(Intel®Xeon®Gold 6140

    CPU @ 2.30GHz) • msvc new: 0.23 [sec] • msvc old: 4.18 [sec] • gcc libstdc++: 0.36 [sec] • clang libc++: 0.21 [sec] • Abseil 0.41 [sec] • 呼び出しごとに、log2(bits)を計算しており遅い。が、bitsはコンパイル時 に決定できるため、無駄。 • 他の実装はコンパイル時に計算するための工夫が行われている。 • msvc 19.32からコンパイル時に計算するようになり、高速化された! • かなり最近。
  11. 除算法 • generate_canonicalのアルゴリズムをそのまま取り出す。 • float r = (float)random() / ((float)(RandomMax)+(float)1);

    • 先述したような問題がある。 • [0.0, 1.0) が保証されない。 • 丸め誤差によって1.0fが出てしまう。これはすでに述べた通り。 • 解決策は、エッジが出たら再生成 or クランプ、など。 • [0.0, 1.0) の範囲における、表現可能な全てのfloat値を得ることが出来ない。 • この問題について考察する。
  12. 除算法:浮動小数の構造について • ビット上の構造 • fractionはmantissaとも呼ぶ 0.0 1.0 0.5 0.25 指数部(exponent)によって範囲を決定

    仮数部(mantissa, fraction)によって位置を決定 −1 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 ⋅ 2𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒−127 ⋅ (1. 𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓𝑓)
  13. 除算法:浮動小数の構造について • ビット上の構造 0.0 1.0 0.5 0.25 指数部(exponent)によって範囲を決定 この場合exponent=125 仮数部(mantissa,

    fraction)によって位置を決定 0.0 1.0 0.5 0.25 指数部(exponent)によって範囲を決定 この場合exponent=126 仮数部(mantissa, fraction)によって位置を決定
  14. 除算法 • 除算法は以下のようなイメージとなる。 0.0 1.0 0.5 0.25 乱数 0.0 1.0

    [0.0, 1.0) 区間を等分し、いずれかを乱数で選び出す 浮動小数表現にマッピングする(丸める) float
  15. 表現可能な値 • 除算法によって得られる、0より大きい最小の値 • float f = (float)1 / ((float)0xffffffff

    + (float)1); • これは 2.32830644e-10 となる。 • 一方、floatで表現可能な、0より大きい最小の値は • 正規化数の範囲なら 1.17549435e-38 • 非正規化数も含むと 1.40129846e-45 • floatの指数部は正規化数の範囲で2^-126(=1.17549435e-38)まで 小さくなる。しかし、除算法で得られる、2.32830644e-10 というのは指数 でいえば2^-29から2^-30の間くらいの精度しかない。
  16. 表現可能な値に関する問題 • exponential distribution • 乱数ξを使って、u=-logn(ξ)によって得られる分布。ボリュームレンダリングと かでも良く使う。 • ξは(0.0, 1.0]の乱数が使われる。

    • ξとして理想的なfloatの範囲をとるならば、最小値は(既に述べたよ うに)正規化数で1.17549435e-38となり、uの最大値は87.34となる。 しかし、除算法ではuの最大値は22.18になる。 • 本来、uの最大値は∞なので、いずれも真の分布になっていないという本質的な 問題は別途存在する。
  17. mantissa(仮数部)法 • [0.0, 1.0) の一様乱数を得る方法として、mantissa(仮数部)を 直接埋める方法がある。 0.0 1.0 0.5 0.25

    指数部(exponent)によって範囲を決定 仮数部(mantissa, fraction)によって位置を決定
  18. mantissa(仮数部)法 • ちゃんと開区間を表現できる。 • 除算法などと違って、単なるビット演算とfloatへの変換のみなので、あ る種のハードウェアでは高速に実行されることが期待される。 • とはいっても、今時のCPUでは、そこまでパフォーマンスは有利ではない。 • 1億回の実行に、0.20

    [sec]程度。 • 得られる、0より大きい最小値は1.19209290e-07となる。 • 表現可能な値に比べると、かなり大きめ。 • 大きな問題として、最終結果の乱数のmantissaにおける最下位ビットが 常にゼロになるという事実がある。 • 表現可能な値が少ない。
  19. mantissa(仮数部)法 • 最終結果の乱数の最下位ビットが常にゼロになる。 • 0.348863: 00111110101100101001111000110100 • 0.971680: 00111111011110001100000000000110 •

    0.060423: 00111101011101110111111001100000 • 0.093260: 00111101101111101111111011100000 • 0.683346: 00111111001011101110111111000110 • 0.770189: 00111111010001010010101100100000 • 0.236677: 00111110011100100101101101101000 • 0.465761: 00111110111011100111100000111000 • 0.824375: 00111111010100110000101000111010 • 0.243600: 00111110011110010111001001101000 • 0.395521: 00111110110010101000000111000100 • 0.182249: 00111110001110101001111110010000 • 0.198098: 00111110010010101101101000100000 • 0.712343: 00111111001101100101110000010100 • 0.576150: 00111111000100110111111010001100 • 0.880710: 00111111011000010111011000111000
  20. mantissa(仮数部)法まとめ 実装が自明、かつ単純。 [0.0, 1.0) が保証される。 doubleにも自然に拡張可能。 [0.0, 1.0) の範囲における、表現可能な全てのfloat値を得ること が出来ない。

    しかも、mantissaの最下位ビットが常にゼロになる。 除算法と比べても、出力される最小値が大きい。 除算法:2.32830644e-10 mantissa法:1.19209290e-07 〇 × × 〇 〇
  21. 除算法2 • 除算法を修正する • float r = (float)(random() >> 8)

    / 16777216.0f; • 改善点 • [0.0, 1.0) が保証される! • 乱数として0xffffffffが出た時、9.99999940e-01となる。 • mantissa法より分解能が高い! • 2^24の分解能となる。mantissa法は2^23。 • 分解能について • 表現可能なfloat値全てを得ることはできない。 • 出力される0より大きい最小の値は、5.9604644775e-08 • 除算法:2.32830644e-10 • mantissa法:1.19209290e-07
  22. exponent+mantissa法 • アルゴリズムの概要 • 1bit乱数を得る。 • 1なら、[0.5, 1.0)区間のexponentに決定。 • 0なら、再び1bit乱数を得る。

    • 1なら、[0.25, 0.5)区間のexponentに決定。 • 0なら、再び1bit乱数を得る。以上を繰り返す。
  23. Abseilの最適化 • Abseilは、exponent+mantissa法の簡略版。 • 単一の64bit乱数で処理を完結させている。 • exponent決定時、反復回数を63回に制限(=乱数のbit数に制限)。 • 何度も乱数を引かない、ということ。exponentの範囲が狭くなり、出力される 最小値が大きくなる。

    • といっても、出力される最小値は5.42101086e-20なので、そこまで悪くない。 • mantissa決定時、exponent決定時に余ったbitを流用。 • exponent次第で、mantissaが常に23bit使い切らない。 • といっても、非常に小さい値の範囲での話なので、そこまで問題にならない。
  24. 結論 • なんだかんだ普通に除算法するのでもあまり問題にならない。 • パフォーマンスは良い。表現可能なfloat値の全ては得られない。 • [0.0, 1.0)を保証するかどうかなどについても考えておく必要がある。 • [0.0,

    1.0)を保証しつつ、表現可能なfloat値の全てを得たいなら、 exponent+mantissa法やDowneyの方法もあり得る。 • ただし、乱数を何度も引く可能性などについて考えておく必要がある。 • Abseilの最適化が良いバランスになることもあり得る。 • 目的やパフォーマンス、得られる乱数のbit数(と乱数生成速度)など に応じて、最適な選択がきまる、かもしれない。
  25. おまけ:std::uniform_real_distributionの問題 • std::uniform_real_distributionを使うことで、任意の範囲 [a, b) の一様分布を得られる。 • が、与える値次第で容易に区間に関する条件を満たさなくなる。 • つまり、bが返ってしまうことがある。

    • [0.0, 1.0)の乱数を用いても、丸めによってbになりえる。 • 一般的に解決するにはおそらく難しく、再生成する、クランプする、などしかない かもしれない。 • 世間の実装も概ねそのように対処しているか、そもそも対処していないか。 • ただし、得たい値がfloatならdouble除算法の結果を最後にキャストするよう にすればある程度精度向上を図れる(今どきなら、パフォーマンスもあまり変わ らない)