使用 getrandom 在 C 中随机浮点数答案

【问题标题】：Random float in C using getrandom使用 getrandom 在 C 中随机浮点数
【发布时间】：2016-06-25 22:29:15
【问题描述】：

我正在尝试生成一个介于 0 和 1 之间的随机浮点数（无论它是在 [0,1] 还是 [0,1）对我来说都无关紧要）。网上关于此的每个问题似乎都涉及rand() 调用，以time(NULL) 为种子，但我希望能够每秒多次调用我的程序并每次都获得不同的随机数。这将我引向 Linux 中的 getrandom 系统调用，它来自 /dev/urandom。我想出了这个：

#include <stdio.h>
#include <sys/syscall.h>
#include <unistd.h>
#include <stdint.h>

int main() {
  uint32_t r = 0;
  for (int i = 0; i < 20; i++) {
    syscall(SYS_getrandom, &r, sizeof(uint32_t), 0);
    printf("%f\n", ((double)r)/UINT32_MAX);
  }
  return 0;
}

我的问题只是我是否正确执行此操作。它似乎有效，但我担心我在滥用某些东西，并且几乎没有在线使用 getrandom() 的示例。

【问题讨论】：

而不是syscall()，你不能打开fopen('/dev/urandom', 'rb') 并读取4 个字节吗？或者然后将其放入srand()？
一种更便携的方法是从中打开/dev/urandom和read(2)。
@user3030010，是什么让您认为系统调用比从 urandom 设备读取更不容易失败？
如果你用随机位填充一个双精度或浮点数，你可以产生一个 NaN 或一个无穷大，但没有人建议这样做。如果您生成一个用随机位填充的无符号整数，转换为 double，然后除以可能的最大值，那么您在闭区间 [0,1] 中可靠地得到 double。
@tofro 她/他说的不是重复播种，而是重复运行程序，程序在启动时播种一次。

标签： c random floating-point

【解决方案1】：

OP 有 2 个问题：

如何随机启动序列。
如何在 [0...1) 范围内生成double。

通常的方法是获取一个非常随机的源，例如 /dev/urandom 或来自 syscall() 或什至可能是 seed = time() ^ process_id; 的结果，并通过 srand() 种子。然后根据需要拨打rand()。

下面包含一个快速转换的方法来生成均匀的[0.0 to 1.0)（线性分布）。但是像所有随机生成函数一样，真正好的随机生成函数是基于广泛的研究。这个简单地基于DBL_MANT_DIG和RAND_MAX调用rand()几次，

[编辑] 原始double rand_01(void) 有一个弱点，它只生成 2^52 个不同的doubles 而不是 2^53。它已被修改。替代方法：double 版本的rand_01_ld(void) 远低于此。

#include <assert.h>
#include <float.h>
#include <limits.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

double rand_01(void) {
  assert(FLT_RADIX == 2); // needed for DBL_MANT_DIG
  unsigned long long limit = (1ull << DBL_MANT_DIG) - 1;
  double r = 0.0;
  do {
    r += rand();
    // Assume RAND_MAX is a power-of-2 - 1
    r /= (RAND_MAX/2 + 1)*2.0;
    limit = limit / (RAND_MAX/2 + 1) / 2;
  } while (limit);

  // Use only DBL_MANT_DIG (53) bits of precision.
  if (r < 0.5) {
    volatile double sum = 0.5 + r;
    r = sum - 0.5;
  }
  return r;
}

int main(void) {
  FILE *istream = fopen("/dev/urandom", "rb");
  assert(istream);
  unsigned long seed = 0;
  for (unsigned i = 0; i < sizeof seed; i++) {
    seed *= (UCHAR_MAX + 1);
    int ch = fgetc(istream);
    assert(ch != EOF);
    seed += (unsigned) ch;
  }
  fclose(istream);
  srand(seed);

  for (int i=0; i<20; i++) {
    printf("%f\n", rand_01());
  }

  return 0;
}

如果想扩展到更宽的 FP，无符号宽整数类型可能不够用。下面是一种没有这种限制的可移植方法。

long double rand_01_ld(void) {
  // These should be calculated once rather than each function call
  // Leave that as a separate implementation problem
  // Assume RAND_MAX is power-of-2 - 1
  assert((RAND_MAX & (RAND_MAX + 1U)) == 0);
  double rand_max_p1 = (RAND_MAX/2 + 1)*2.0;
  unsigned BitsPerRand = (unsigned) round(log2(rand_max_p1));
  assert(FLT_RADIX != 10);
  unsigned BitsPerFP = (unsigned) round(log2(FLT_RADIX)*LDBL_MANT_DIG);

  long double r = 0.0;
  unsigned i;
  for (i = BitsPerFP; i >= BitsPerRand; i -= BitsPerRand) {
    r += rand();
    r /= rand_max_p1;
  }
  if (i) {
    r += rand() % (1 << i);
    r /= 1 << i;
  }
  return r;
}

【讨论】：

我正在考虑以比结果具有尾数位更多的随机位开始对精度的影响。为什么不直接以合适的整数生成DBL_MANTISSA_BITS 随机位，将其转换为double，然后使用ldexp() 将其缩放到[0, 1) 的范围内？
@John Bollinger 面临的挑战是以便携方式对RAND_MAX和DBL_MANT_DIG进行“生成DBL_MANTISSA_BITS随机位”。当然希望代码在编译时知道log2(RAND_MAX + 1) 的整数常量，而不会出现任何范围错误或其他问题。即使RAND_MAX +1u 也会溢出。任何想法表示赞赏。
@John Bollinger 添加了第二个方法（用于long double），不使用额外的随机位。
不错。如果可以的话，我会再给你一个 +1。
另一个内容丰富的帖子。有什么反对float.h 和limits.h 的吗？

【解决方案2】：

如果您需要生成双打，可以使用以下算法：

CPython generates random numbers 使用以下算法（我更改了函数名称、类型定义和返回值，但算法保持不变）：

double get_random_double() {
    uint32_t a = get_random_uint32_t() >> 5;
    uint32_t b = get_random_uint32_t() >> 6;
    return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
}

该算法的来源是 Takuji Nishimura 和 Makoto Matsumoto 的 Mersenne Twister 19937 随机数生成器。不幸的是，源中提到的原始链接不再可供下载。

CPython 中对该函数的注释如下：

【这个函数】就是原代码中名为genrand_res53的函数；在 [0,1) 上生成一个 53 位分辨率的随机数；注意 9007199254740992 == 2**53;我假设他们将“/2**53”拼写为乘以倒数（可能是徒劳的）希望编译器将在编译时优化除法。 67108864 是 2**26。在效果，a 包含 27 个左移 26 的随机位，b 填充 53 位分子的低 26 位。

该算法的原始代码归功于 Isaku Wada，2002/01/09

从该代码简化，如果您想快速创建float，您应该用(1 << FLT_MANT_DIG) - 1 屏蔽uint32_t 的位并除以(1 << FLT_MANT_DIG) 以获得正确的[0, 1) 间隔：

#include <stdio.h>
#include <sys/syscall.h>
#include <unistd.h>
#include <stdint.h>
#include <float.h>

int main() {
    uint32_t r = 0;
    float result;
    for (int i = 0; i < 20; i++) {
        syscall(SYS_getrandom, &r, sizeof(uint32_t), 0);
        result = (float)(r & ((1 << FLT_MANT_DIG) - 1)) / (1 << FLT_MANT_DIG);
        printf("%f\n", result);
    }
    return 0;
}

由于可以假设您的 Linux 具有 C99 编译器，我们可以使用 ldexpf 来代替那个除法：

#include <math.h>

result = ldexpf(r & ((1 << FLT_MANT_DIG) - 1), -FLT_MANT_DIG);

要获得闭区间[0, 1]，可以做效率稍低的操作

result = ldexpf(r % (1 << FLT_MANT_DIG), -FLT_MANT_DIG);

为了快速生成大量高质量的随机数，我只需使用系统调用来获取足够的数据来播种 PRNG 或 CPRNG，然后从那里继续。

【讨论】：