Как создать равномерно распределенный случайный дубль в C, используя libsodium в диапазоне [-a, a]? - PullRequest
Новая
       14

Как создать равномерно распределенный случайный дубль в C, используя libsodium в диапазоне [-a, a]?

0 голосов
/ 25 апреля 2019

В библиотеке libsodium есть функция

uint32_t randombytes_uniform(const uint32_t upper_bound);

, но, очевидно, это возвращает целое число без знака.Могу ли я каким-то образом использовать это для генерации равномерно распределенного случайного числа double в диапазоне [-a,a], где a также является двойным, заданным пользователем?Я особенно сосредоточен на том, чтобы результат был равномерно распределен / несмещен, поэтому я хотел бы использовать библиотеку libsodium.

Ответы [ 3 ]

1 голос
/ 25 апреля 2019

Позвольте мне попытаться сделать это шаг за шагом.

Во-первых, вам, очевидно, нужно объединить два вызова, чтобы получить до 64 бит случайности для одного выхода с двойным значением.

Во-вторых, вы конвертируете его в интервал [0 ... 1].Есть несколько способов сделать это, все они в каком-то смысле хороши, я предпочитаю единообразные случайные двоичные рациональные числа в форме n * 2 -53 , подробности см. здесь ,Вы можете попробовать и другие методы, перечисленные выше.Примечание: методы в ссылке дают результаты в диапазоне [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;
}
тривиально
1 голос
/ 25 апреля 2019 
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);
0 голосов
/ 26 апреля 2019

Сначала признается, что нахождение случайного числа [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.

Некоторые проблемы округления применяются, что этот ответ замаскирован, но распределение остается однородным - за исключением потенциально на краях.

...