C++ の std::uniform_real_distribution はあまり信用できない。

これは何？

C++11 で導入された uniform_real_distribution。
仕様と実装で不思議なことがあったので書いておく。

仕様

C++ の仕様書のようなもの を見ると

produces random numbers x, a ≤ x < b (引用者註: a, b は 2 引数コンストラクタの引数)

とあるのに

Preconditions: a ≤ b and b − a ≤ numeric_limits<RealType>::max().

となっている。

Preconditions では a=b でも良さそうなことが書いてあるが、作られる値を見ると a=b だと困る。
仕様書のバグかな？ それとも何か見落としている？

色々試す

上端が出る例

まあソースコード。

c++11

#include <cstdio>
#include <random>

int main(int argc, char *argv[])
{
  std::mt19937_64 engine(0);
  auto delta = 1.0/(1ull<<50);
  auto maxval = 1+delta;
  auto minval = 1-delta;
  std::uniform_real_distribution<> dist(minval,maxval);

  constexpr unsigned long long TRIAL = 1ull<<28;
  for (unsigned long long i = 0; i<TRIAL; i++)
  {
    auto v = dist(engine);
    if (v<minval || maxval<=v ){
      std::printf("v=%.17f  range is [%.17f, %.17f)", v, minval, maxval );
      return 0;
    }
  }
  std::puts( "No unexpected values were found." );
  return 0;
}

仕様通りなら No unexpected values were found. と出るはずだが、手元の以下のコンパイラは全部 v=1.00000000000000089 range is [0.99999999999999911, 1.00000000000000089) と出力して終了した。

• Apple clang version 13.0.0 (clang-1300.0.29.3) on macOS
• g++-11 (Homebrew GCC 11.2.0) 11.2.0 on macOS
• g++-7 (Homebrew GCC 7.5.0_4) 7.5.0 on macOS
• clang version 7.0.1-8+rpi3+deb10u2 on Raspberry Pi OS
• g++-8 (Raspbian 8.3.0-6+rpi1) 8.3.0 on Raspberry Pi OS

上端を超える例

先程は上端ぴったりが出る例だったけど、上端を超える例もある。

これ

c++11

#include <iostream>
#include <random>
#include <cmath>
#include <cstdio>
#include <utility>

int main(int argc, char *argv[])
{
  using real = float;
  std::mt19937_64 rng(0);
  constexpr real d0 = 1e-36;
  constexpr real d1 = d0 * (1ull << 22);
  real minval = -d1 * 16;
  real maxval = d1 + d0;
  std::cout << "[" << minval << ", " << maxval << ")\n";
  std::uniform_real_distribution<real> dist(minval, maxval);
  constexpr unsigned long long TRIAL = 1ull << 28;
  for (unsigned long long i = 0; i < TRIAL; i++)
  {
    auto v = dist(rng);
    if (v < minval || maxval < v)
    {
      std::printf("v=%.10e  range is [%.10e, %.10e)\n", v, minval, maxval);
      return 0;
    }
  }
  std::puts("No unexpected values were found.");
  return 0;
}

を Apple clang on macOS で実行すると

[-6.71089e-29, 4.19431e-30)
v=4.1943082885e-30  range is [-6.7108866412e-29, 4.1943052792e-30)

などと表示される。試した範囲内では、 Apple clang on macOS 以外は No unexpected values were found. となる。

なにが起こっているのか

g++ も clang も、ソースコードを見ると、シンプルに

rand_val() * (b-a) + a

というような計算をしている。rand_val は、0≤rand_val()<1 の一様乱数だと思うけど、詳細はちょっとわからない。

この計算だと

• b-a+a != b の場合に上端があやふやになる
• a 〜 b の間で取りうる値の数が rand_val() がとりうる値の数よりもずっと少ない場合、上端があやふやになる

という2つの問題が起こる。

clang on macOS とそれ以外で何が違うのかはよくわからなかった。

そもそも

そもそも、浮動小数点数の一様乱数は難しい。

たとえば [a,b) の値を出す倍精度の一様乱数を考える。これが a<0<b なら 0 を出すかもしれないことを保証するのは大変だと思う。
あるいは [-1e30,1) の値を出す倍精度の一様乱数。これが正の値を出すかもしれないことを保証するのは大変だと思う。

理想的には、$a≤x<b$ の浮動小数点数（-0.0 は含めないものとする）を小さい順に整列したものを $v_0, v_1, ..., v_{n-1}$ とし、$v_n = b$ としたときに。
$v_n$ が出る確率が $(v_{n+1}-v_n)÷(b-a)$ であるような乱数がいい。

でもどうやったら安く実現できるだろうか。

対策

現実としては a≤x<b の乱数がほしいのに b=x だったり b<x だったりするのは本当に困る場合がある。

対策としては、 uniform_real_distribution を使わない、しかないと思う。
整数の一様乱数から浮動小数点数の乱数を自力で生成するのが良い。

素朴に計算すると丸められてすぐに b≤x の値を生成してしまうので、そこは慎重に実装する必要がある。
[-1e30, 1) の単精度の値が必要で、たまには正の値も出てほしい、みたいなのを倍精度を使わずに実装する必要があるならわりとたくさん書く必要があると思うし、そんなに簡単ではないと思う。