基于随机比特流生成随机浮点值

Question

给定一个随机源（随机比特流的生成器），如何在给定范围内生成一个均匀分布的随机浮点值？

假设我的随机源看起来像：

unsigned int GetRandomBits(char* pBuf, int nLen);

我想实现

double GetRandomVal(double fMin, double fMax);

注意事项：

• 我不希望结果精度受到限制（例如只有 5 位）。
• 必须严格均匀分布
• 我不是要求对现有库的引用。我想知道如何从头开始实施它。
• 对于伪代码/代码，C++ 将不胜感激

Answer 1

这个问题是不恰当的。浮点数上的均匀分布是什么意思？

根据discrepancy 的提示，解决您的问题的一种方法是定义您想要最小化以下值的分布：

其中x 是您使用GetRandomVal(double fMin, double fMax) 函数采样的random variable， 表示随机x 小于或等于t 的概率。

现在您可以继续尝试评估例如a dabbler's answer。 （提示所有未能使用整个精度并坚持例如 52 位的答案将不符合此最小化标准。）

但是，如果您只是希望能够以相同的可能性生成落在您指定范围内的所有浮点位模式，即使这意味着例如请求GetRandomVal(0,1000) 将在 0 和 1.5 之间创建比在 1.5 之间更多的值和 1000，这很容易：当解释为位模式时，IEEE 浮点数的任何间隔都可以轻松映射到 unsigned int64 的极少数间隔。参见例如question。在任何给定间隔内生成均匀分布的unsigned int64 随机值很容易。

Answer 2

令我惊讶的是，对于这个古老的问题，没有人提供最佳答案的实际代码。 User515430's answer 做对了——您可以利用 IEEE-754 双精度格式直接将 52 位放入双精度中，而无需任何数学运算。但他没有给出代码。所以在这里，来自我的公共领域ojrandlib：

double ojr_next_double(ojr_generator *g) {
    uint64_t r = (OJR_NEXT64(g) & 0xFFFFFFFFFFFFFull) | 0x3FF0000000000000ull;
    return *(double *)(&r) - 1.0;
}

NEXT64() 获取一个 64 位随机数。如果您有更有效的方法来仅获取 52 位，请改用它。

Answer 3

我想我永远不会相信你真的需要这个，但写起来很有趣。

#include <stdint.h>

#include <cmath>
#include <cstdio>

FILE* devurandom;

bool geometric(int x) {
  // returns true with probability min(2^-x, 1)
  if (x <= 0) return true;
  while (1) {
    uint8_t r;
    fread(&r, sizeof r, 1, devurandom);
    if (x < 8) {
      return (r & ((1 << x) - 1)) == 0;
    } else if (r != 0) {
      return false;
    }
    x -= 8;
  }
}

double uniform(double a, double b) {
  // requires IEEE doubles and 0.0 < a < b < inf and a normal
  // implicitly computes a uniform random real y in [a, b)
  // and returns the greatest double x such that x <= y
  union {
    double f;
    uint64_t u;
  } convert;
  convert.f = a;
  uint64_t a_bits = convert.u;
  convert.f = b;
  uint64_t b_bits = convert.u;
  uint64_t mask = b_bits - a_bits;
  mask |= mask >> 1;
  mask |= mask >> 2;
  mask |= mask >> 4;
  mask |= mask >> 8;
  mask |= mask >> 16;
  mask |= mask >> 32;
  int b_exp;
  frexp(b, &b_exp);
  while (1) {
    // sample uniform x_bits in [a_bits, b_bits)
    uint64_t x_bits;
    fread(&x_bits, sizeof x_bits, 1, devurandom);
    x_bits &= mask;
    x_bits += a_bits;
    if (x_bits >= b_bits) continue;
    double x;
    convert.u = x_bits;
    x = convert.f;
    // accept x with probability proportional to 2^x_exp
    int x_exp;
    frexp(x, &x_exp);
    if (geometric(b_exp - x_exp)) return x;
  }
}

int main() {
  devurandom = fopen("/dev/urandom", "r");
  for (int i = 0; i < 100000; ++i) {
    printf("%.17g\n", uniform(1.0 - 1e-15, 1.0 + 1e-15));
  }
}

Answer 4

这是一种方法。

IEEE Std 754 双精度格式如下：

[s][     e     ][                          f                         ]

其中 s 是符号位（1 位），e 是偏置指数（11 位），f 是小数（52 位）。

注意内存中的布局在 little-endian 机器上会有所不同。

对于0

(-1)**(s)   *  2**(e – 1023)  *  (1.f)

通过将 s 设置为 0，将 e 设置为 1023，将 f 设置为比特流中的 52 个随机位，您将在区间 [1.0, 2.0) 中获得一个随机双精度。这个区间的独特之处在于它包含 2 ** 52 个双精度数，并且这些双精度数是等距的。如果然后从构造的双精度数中减去 1.0，则会在区间 [0.0, 1.0) 中得到一个随机双精度数。此外，关于等距的性质是保留的。 从那里您应该能够根据需要进行缩放和翻译。

Answer 5

这可能不是你想要的答案，而是这里的规范：

http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2010/n3225.pdf

在 [rand.util.canonical] 和 [rand.dist.uni.real] 部分中，包含足够的信息来实现您想要的，尽管语法略有不同。这并不容易，但这是可能的。我从个人经历说。一年前，我对随机数一无所知，但我能够做到。虽然花了我一段时间... :-)

Answer 6

要在 [0..1[ 中获取随机值，您可以执行以下操作：

double value = 0;
for (int i=0;i<53;i++)
   value = 0.5 * (value + random_bit());  // Insert 1 random bit
   // or value = ldexp(value+random_bit(),-1);
   // or group several bits into one single ldexp
return value;

Answer 7

我可能误解了这个问题，但是是什么阻止了您从随机比特流中采样接下来的 n 位并将其转换为范围为 0 到 2^n - 1 的以 10 为基数的数字。

Answer 8

这很容易，只要您有一个精度与double 一样多的整数类型。例如，一个 IEEE 双精度数有 53 位精度，所以 64 位整数类型就足够了：

#include <limits.h>
double GetRandomVal(double fMin, double fMax) {
  unsigned long long n ;
  GetRandomBits ((char*)&n, sizeof(n)) ;
  return fMin + (n * (fMax - fMin))/ULLONG_MAX ;
}