如何使用 [-a,a] 范围内的 libsodium 在 C 中生成均匀分布的随机双精度数？答案

【问题标题】：How to generate a uniformly distributed random double in C using libsodium in range [-a,a]?如何使用 [-a,a] 范围内的 libsodium 在 C 中生成均匀分布的随机双精度数？
【发布时间】：2019-04-25 13:49:17
【问题描述】：

libsodium 库有一个函数

uint32_t randombytes_uniform(const uint32_t upper_bound);

但显然这会返回一个无符号整数。我可以以某种方式使用它在[-a,a] 范围内生成均匀分布的随机double，其中a 也是用户给出的双精度数？我特别关注均匀分布/无偏见的结果，这就是我想使用libsodium 库的原因。

【问题讨论】：

rand() 有什么问题？
我担心它的均匀分布的保证，重要的是生成的结果尽可能不带偏见。
如果你搜索“random double”，你会发现很多关于这个的其他问题。但是你的问题定义不明确。浮点格式在零附近比远离零更好。如果您从一个区间内的均匀分布中选择一个实数，并将其四舍五入为可在double 中表示的最接近的值，则某些值会比其他值更频繁地出现，因为它们离邻居更远，而实数的跨度对它们的舍入比double 接近零的值要长。
此外，端点在区间内的频率将是其邻居频率的一半。有时在选择“随机”浮点值时，人们希望使用均匀间隔的值群，忽略更精细值。这取决于输入数字的目的。你用随机数做什么？
@EricPostpischil 我看不出一个数学上简单的任务，比如“选择一个范围内均匀分布的随机实数”是如何定义不明确的。我不精通数字精度和实数表示等主题，所以这项任务不可能吗？另外，我想如果这意味着我可以选择“粒度”并且结果将是0.001 的倍数，我们可以选择使用“均匀分布的值的群体，忽略更精细的值”均匀分布在 [-1.245, 3.523] 范围内？

标签： c random double libsodium

【解决方案1】：

const uint32_t mybound = 1000000000; // Example
const uint32_t x = randombytes_uniform(mybound);
const double a = 3.5; // Example
const double variate = a * ( (2.0 * x / mybound) - 1);

【讨论】：

我也希望绑定为双精度，例如结果在 [-3.5,+3.5] 范围内，这可能吗？
是的，我给你举了几个数字作为例子。请考虑点击接受作为答案，谢谢。
变量mybound的值有什么意义，比如可以设置为UINT_MAX吗？
这是假的。要使用 53 位尾数使 64 位加倍，您只需使用 32 个随机位。你怎么能保证精确地在 [-a...a] 范围内并且覆盖均匀？

【解决方案2】：

让我试着一步一步来。

首先，您显然需要结合两个调用来为一个双精度值输出获得高达 64 位的随机性。

其次，将其转换为 [0...1] 区间。有几种方法可以做到这一点，所有这些在某种意义上都很好，我更喜欢 n*2^-53 形式的统一随机二元有理数，详情请参阅here。您也可以尝试上面列出的其他方法。注意：链接中的方法在 [0...1) 范围内产生结果，我尝试接受/拒绝以关闭 [0...1] 范围。

最后，我将结果缩放到所需的范围内。

抱歉，仅限 C++，但转换为 C 很简单

#include <stdint.h>
#include <math.h>
#include <iostream>
#include <random>

// emulate libsodium RNG, valid for full 32bits result only!
static uint32_t randombytes_uniform(const uint32_t upper_bound) {
    static std::mt19937 mt{9876713};
    return mt();
}

// get 64bits from two 32bit numbers
static inline uint64_t rng() {
    return (uint64_t)randombytes_uniform(UINT32_MAX) << 32 | randombytes_uniform(UINT32_MAX);
}

const  int32_t bits_in_mantissa = 53;
const uint64_t max  = (1ULL << bits_in_mantissa);
const uint64_t mask = (1ULL << (bits_in_mantissa+1)) - 1;

static double rnd(double a, double b) {
    uint64_t r;
    do {
        r = rng() & mask; // get 54 random bits, need 53 or max
    } while (r > max);
    double v = ldexp( (double)r, -bits_in_mantissa ); // http://xoshiro.di.unimi.it/random_real.c
    return a + (b-a)*v;
}

int main() {
    double a = -3.5;
    double b = 3.5;

    for(int k = 0; k != 100; ++k)
        std::cout << rnd(a, b) << '\n';

    return 0;
}

【讨论】：

【解决方案3】：

首先认识到找到一个随机数[0...a] 就足够了，然后是+/- 的硬币翻转。

第 2 步。找到expo 使得a < 2**expo 或ceil(log2(a))。

int sign;
do {
  int exp;
  frexp(a, &exp);

步骤 3. 形成一个完整的 63 位随机数 [0...0x7FFF_FFFF_FFFF_FFFF] 和随机符号。 63 应该至少与 double 的精度一样宽 - 通常是 53 位。此时r 肯定是统一。

  unit64_t r = randombytes_uniform(0xFFFFFFFF);
  r <<= 32;
  r |= randombytes_uniform(0xFFFFFFFF);
  // peel off one bit for sign
  sign = r & 1;
  r >>= 1;

第 4 步。缩放并测试是否在范围内。根据需要重复。

  double candidate = ldexp(r/pow(2 63), expo);
} while (candidate > a);

第 5 步。应用标志。

if (sign) {
  candidate = -candidate;
}
return candidate;

避免使用(2.0 * x / a) - 1，因为计算不是关于 0.0 对称的。

代码将受益于处理DBL_MAX 附近的a 的改进。

此答案掩盖了一些舍入问题，但分布保持均匀 - 可能在边缘除外。

【讨论】：