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

Question

libsodium库有一个功能

uint32_t randombytes_uniform(const uint32_t upper_bound);

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

Answer 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);

Answer 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;
}

Answer 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不对称。

---

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

一些舍入问题适用于这个答案掩盖，但分布仍然是均匀的 - 除了潜在的边缘。