C语言中使用getrandom生成随机浮点数

Question

我想生成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()的示例。

- user3030010

---

您可以尝试使用fopen('/dev/urandom', 'rb')打开文件并读取4个字节，而不是使用syscall()函数。然后将这4个字节传递给srand()函数。这样做可以实现同样的随机数生成效果。 - chux - Reinstate Monica

一个更可移植的方法是打开 /dev/urandom 并从中读取 read(2)。 - Nate Eldredge

@user3030010，您认为系统调用比从urandom设备读取更不容易失败的原因是什么？ - John Bollinger

如果你使用随机位填充一个双精度或单精度浮点数，可能会产生NaN或无穷大，但是没有人建议这样做。如果你生成一个填满随机位的无符号整数，将其转换为double，然后除以最大可能值，那么你可靠地得到一个在闭区间[0,1]内的double。 - John Bollinger

@tofro 他/她说的不是重复播种，而是程序重复运行，在启动时程序只播种一次。 - user253751

Answer 1

OP有两个问题：

1. 如何随机开始序列。

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 种不同的 double，而不是 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;
}

- chux - Reinstate Monica

---

我在思考使用比结果具有更多尾数位的随机位开始会对精度产生什么影响。为什么不在一个适当的整数中生成DBL_MANTISSA_BITS个随机位，将其转换为double，并使用ldexp()将其缩放到范围[0, 1)中呢？ - John Bollinger

@John Bollinger 添加了第二种方法（针对long double），不使用额外的随机位。 - chux - Reinstate Monica

不错。如果可以的话，我会再给你一个+1。 - John Bollinger

又一篇有用的帖子。关于float.h和limits.h有什么反对意见吗？ - David C. Rankin

@David C. Rankin float.h和limits.h来晚了，它们并不是要被“排除”的。 - chux - Reinstate Monica

Answer 2

如果您需要生成双精度浮点数，以下算法可能会有用：

CPython使用以下算法生成随机数（我更改了函数名称、typedef和返回值，但算法保持不变）：

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的函数；生成一个具有53位分辨率的[0,1)之间的随机数；请注意，9007199254740992 == 2 ** 53；我假设他们将“/2**53”拼写成乘以倒数的形式，以便编译器可以在编译时优化掉除法。67108864等于2 ** 26。实际上，a包含了27个左移26位的随机位，而b填充53位分子的较低的26位。

原始代码将Isaku Wada归功于这个算法，时间是2002年1月9日。

---

简化代码后，如果你想快速创建一个float，你应该使用uint32_t的位掩码和(1 << FLT_MANT_DIG) - 1除以(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，然后从那里继续。

- Antti Haapala -- Слава Україні