Генерация случайных значений с плавающей точкой на основе случайного битового потока

Question

Учитывая случайный источник (генератор случайного потока битов), как мне сгенерировать равномерно распределенное случайное значение с плавающей точкой в ​​заданном диапазоне?

Предположим, что мой случайный источник выглядит примерно так:

unsigned int GetRandomBits(char* pBuf, int nLen);

И я хочу реализовать

double GetRandomVal(double fMin, double fMax);

Примечания:

• Я не хочу, чтобы точность результата была ограничена (например, только 5 цифрами).
• Необходимо строгое равномерное распределение
• Я не прошу ссылку на существующую библиотеку. Я хочу, чтобы знал, как реализовать это с нуля.
• Для псевдокода / кода C ++ был бы наиболее признателен

Answer 1

Я не думаю, что когда-либо буду убежден, что вам это действительно нужно, но писать было весело.

#include <stdint.h>

#include <cmath>
#include <cstdio>

FILE* devurandom;

bool geometric(int x) {
  // returns true with probability min(2^-x, 1)
  if (x <= 0) return true;
  while (1) {
    uint8_t r;
    fread(&r, sizeof r, 1, devurandom);
    if (x < 8) {
      return (r & ((1 << x) - 1)) == 0;
    } else if (r != 0) {
      return false;
    }
    x -= 8;
  }
}

double uniform(double a, double b) {
  // requires IEEE doubles and 0.0 < a < b < inf and a normal
  // implicitly computes a uniform random real y in [a, b)
  // and returns the greatest double x such that x <= y
  union {
    double f;
    uint64_t u;
  } convert;
  convert.f = a;
  uint64_t a_bits = convert.u;
  convert.f = b;
  uint64_t b_bits = convert.u;
  uint64_t mask = b_bits - a_bits;
  mask |= mask >> 1;
  mask |= mask >> 2;
  mask |= mask >> 4;
  mask |= mask >> 8;
  mask |= mask >> 16;
  mask |= mask >> 32;
  int b_exp;
  frexp(b, &b_exp);
  while (1) {
    // sample uniform x_bits in [a_bits, b_bits)
    uint64_t x_bits;
    fread(&x_bits, sizeof x_bits, 1, devurandom);
    x_bits &= mask;
    x_bits += a_bits;
    if (x_bits >= b_bits) continue;
    double x;
    convert.u = x_bits;
    x = convert.f;
    // accept x with probability proportional to 2^x_exp
    int x_exp;
    frexp(x, &x_exp);
    if (geometric(b_exp - x_exp)) return x;
  }
}

int main() {
  devurandom = fopen("/dev/urandom", "r");
  for (int i = 0; i < 100000; ++i) {
    printf("%.17g\n", uniform(1.0 - 1e-15, 1.0 + 1e-15));
  }
}

Answer 2

Вот один из способов сделать это.

Двойной формат IEEE Std 754 выглядит следующим образом:

[s][     e     ][                          f                         ]

, где s - бит знака (1 бит), e - битсмещенный показатель (11 бит), а f - дробь (52 бита).

Помните, что расположение в памяти на машинах с прямым порядком байтов будет другим.

Для 0

(-1)**(s)   *  2**(e – 1023)  *  (1.f)

Установив s в 0, e в 1023 и f в 52 случайных бита из вашего битового потока, вы получите случайныйдвойной в интервале [1.0, 2.0).Этот интервал уникален тем, что он содержит 2 ** 52 двойных числа, и эти двойные числа равноудалены.Если затем вычесть 1.0 из построенного двойника, вы получите случайный дубль в интервале [0.0, 1.0).Более того, свойство быть равноудаленным - это сохранить.Оттуда вы сможете масштабировать и переводить по мере необходимости.

Answer 3

Я удивлен, что на этот старый вопрос никто не имел действительного кода для лучшего ответа. Ответ User515430 понял все правильно - вы можете воспользоваться двойным форматом IEEE-754, чтобы напрямую сложить 52 бита в двойной без какой-либо математики.Но он не дал код.Так вот, из моего публичного домена ojrandlib :

double ojr_next_double(ojr_generator *g) {
    uint64_t r = (OJR_NEXT64(g) & 0xFFFFFFFFFFFFFull) | 0x3FF0000000000000ull;
    return *(double *)(&r) - 1.0;
}

NEXT64 () получает 64-битное случайное число.Если у вас есть более эффективный способ получить только 52 бита, используйте его вместо этого.

Answer 4

Это легко, если у вас целочисленный тип с таким количеством битов точности, как double.Например, число двойной точности IEEE имеет точность 53 бита, поэтому достаточно 64-разрядного целочисленного типа:

#include <limits.h>
double GetRandomVal(double fMin, double fMax) {
  unsigned long long n ;
  GetRandomBits ((char*)&n, sizeof(n)) ;
  return fMin + (n * (fMax - fMin))/ULLONG_MAX ;
}

Answer 5

Вероятно, это не тот ответ, который вам нужен, а спецификация здесь:

http://www.open -std.org / ОТК1 / SC22 / wg21 / документы / документы / 2010 / n3225.pdf

в разделах [rand.util.canonical] и [rand.dist.uni.real] содержит достаточно информации для реализации того, что вы хотите, хотя и с немного другим синтаксисом. Это не легко, но возможно. Я говорю из личного опыта. Год назад я ничего не знал о случайных числах и смог это сделать. Хотя это заняло у меня некоторое время ...: -)

Answer 6

Вопрос некорректен.Что же означает равномерное распределение по числам с плавающей запятой?

Если взять реплику из несоответствие , то одним из способов реализовать ваш вопрос является определение того, что вы хотите распределение, которое минимизирует следующее значение:

Где x - это случайная величина , в которой вы производите выборку с помощью функции GetRandomVal(double fMin, double fMax), а означает вероятность того, что случайное значение xменьше или равно t.

И теперь вы можете продолжить и попытаться оценить, например, ответ любителя .(Подсказка на все ответы, которые не используют всю точность и придерживаются, например, 52 битов, не выполнит этот критерий минимизации.)

Однако, если вы просто хотите иметь возможность генерировать все шаблоны битов с плавающей запятой, которые попадают в вашуказанный диапазон с равной вероятностью, даже если это означает, что, например, запрос GetRandomVal(0,1000) создаст больше значений от 0 до 1,5, чем от 1,5 до 1000, это легко: любой интервал чисел с плавающей запятой IEEE, если интерпретировать его как битовые комбинации, легко отображается наочень небольшое количество интервалов unsigned int64.Смотрите, например, этот вопрос .Генерация равномерно распределенных случайных значений unsigned int64 в любом заданном интервале очень проста.

Answer 7

Чтобы получить случайное значение в [0..1 [вы можете сделать что-то вроде:

double value = 0;
for (int i=0;i<53;i++)
   value = 0.5 * (value + random_bit());  // Insert 1 random bit
   // or value = ldexp(value+random_bit(),-1);
   // or group several bits into one single ldexp
return value;

Answer 8

Возможно, я неправильно понимаю вопрос, но что мешает вам просто выбрать следующие n битов из случайного потока битов и преобразовать их в число из 10 основных чисел в диапазоне от 0 до 2 ^ n - 1.