Случайный float в C с использованием getrandom

Ответ 1

У OP есть 2 проблемы:

• Как начать последовательность очень случайным образом.

• Как сгенерировать double в диапазоне [0... 1].

Обычный метод заключается в том, чтобы взять очень случайный источник, например /dev/urandom, или результат из syscall() или, может быть, даже seed = time() ^ process_id; и семя через srand(). Затем вызовите rand() по мере необходимости.

Ниже приведено быстро преобразованный метод для создания однородного [0.0 to 1.0) (линейного распределения). Но, как и все случайные генерирующие функции, действительно хорошие основы являются основой обширного изучения. Это просто несколько раз вызывает rand() на основе DBL_MANT_DIG и RAND_MAX,

[Edit] Оригинал 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;
}

Ответ 2

Если вам нужно генерировать удвоения, может потребоваться следующий алгоритм:

CPython генерирует случайные числа, используя следующий алгоритм (я изменил имя функции, typedefs и возвращаемые значения, но алгоритм остается тем же):

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

Источником этого алгоритма является генератор случайных чисел Mersenne Twister 19937 от Takuji Nishimura и Makoto Matsumoto. К сожалению, исходная ссылка, упомянутая в источнике, больше не доступна для скачивания.

Комментарий к этой функции в CPython отмечает следующее:

[эта функция] - это функция genrand_res53 в исходном коде; генерирует случайное число на [0,1) с 53-битным разрешением; Обратите внимание, что 9007199254740992 == 2**53; Я предполагаю, что они написали "/2**53" как умножить на взаимность в (вероятно, напрасно) надеемся, что компилятор будет оптимизируйте разделение во время компиляции. 67108864 - 2**26. В эффект, a содержит 27 случайных битов, сдвинутых влево 26, и b заполняет более низкие 26 бит 53-битного числителя.

Оригинальный код, зачисленный Исаку Вада за этот алгоритм, 2002/01/09

Упрощение из этого кода, если вы хотите быстро создать 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 и перейти оттуда.