Предположим, у вас есть генератор обычных чисел, который способен генерировать равномерно распределенные случайные 32-битные числа. И предположим, что вы ищете способ для генерации псевдослучайной последовательности битов, в которой единицы (то есть установленные биты) появляются в последовательности с заранее определенной вероятностью.
Наивным способом создания такой последовательности будет бегущее число генератор на уровне битов, но он ужасно неэффективен для малых вероятностей, таких как 0,01 или 1% битов в последовательности, большинство битов будут равны нулю. В среднем будет установлен только один бит из ста. С другой стороны, даже при такой низкой вероятности есть шанс встретить длинную подпоследовательность последовательных последовательностей, которая выходит за пределы 8, 16, 32, 64 бит.
Вопрос заключается в том, как эффективно создать такую последовательность с использованием обычного PRNG.
Edit
Игрушечная реализация рациональной выборки переменных Бернулли в javascript, предложенная Питером О. ::11012*
// Based on
// https://arxiv.org/abs/1304.1916
// https://arxiv.org/pdf/1304.1916.pdf (page 21, figure 6)
class Xor128 {
constructor(x, y, z, w) {
this.x = x;
this.y = y;
this.z = z;
this.w = w;
}
prev() {
var t = this.w ^ this.z ^ (this.z >>> 19);
t ^= t >>> 8;
t ^= t >>> 16;
this.w = this.z;
this.z = this.y;
this.y = this.x;
t ^= t << 11;
t ^= t << 22;
this.x = t;
return this.w;
}
curr() {
return this.w;
}
next() {
var t = this.x ^ (this.x << 11);
this.x = this.y;
this.y = this.z;
this.z = this.w;
return this.w = this.w ^ (this.w >>> 19) ^ (t ^ (t >>> 8));
}
}
function* flip(xor128) {
while (true) {
var value = xor128.next();
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
yield value & 1; value >>>= 1;
}
}
function* bernoulli(flip, k, n) {
var b;
var v = k
for (const bit of flip) {
v <<= 1;
if (v >= n) {
v -= n;
b = 1;
} else {
b = 0;
}
if (bit === 1) {
yield b;
v = k;
}
}
}
var xor128 = new Xor128(1, 2, 3, 4);
var z = 0, o = 0;
var then = Date.now();
for (const value of bernoulli(flip(xor128), 5, 1000)) {
if (value === 0) {
z++;
} else {
o++;
}
if (Date.now() - then > 1000) {
console.log(`${z} ${o}`);
}
}
// Pieces of code to test out xor128:
//
// for (let index = 0; index < 100; index++) {
// console.log(xor128.curr())
// xor128.next();
// }
// console.log('-----------------------------------')
// for (let index = 0; index < 100; index++) {
// xor128.prev();
// console.log(xor128.curr())
// }
Другое редактирование
Приведенный ниже код реализован в C#, генерирующем 91,2 миллиона бит в секунду, упакованных в тип данных UInt32 (MacBookPro 2019 Core I9 2,4 ГГц). Я думаю, что в C можно было бы получить более 100 миллионов битов, и такое чувство, что возможно дальнейшее использование двоичного арифмети c для генерации всех 32 битов случайного числа параллельно, некоторое развертывание l oop или, может быть, SIMD не уверен, в любом случае вот код:
public class Bernoulli
{
public UInt32 X { get; set; }
public UInt32 Y { get; set; }
public UInt32 Z { get; set; }
public UInt32 W { get; set; }
public Bernoulli()
: this(Guid.NewGuid())
{
}
public Bernoulli(Guid guid)
{
var index = 0;
var bytes = guid.ToByteArray();
X = (UInt32)((bytes[index++] << 24) | (bytes[index++] << 16) | (bytes[index++] << 8) | bytes[index++]);
Y = (UInt32)((bytes[index++] << 24) | (bytes[index++] << 16) | (bytes[index++] << 8) | bytes[index++]);
Z = (UInt32)((bytes[index++] << 24) | (bytes[index++] << 16) | (bytes[index++] << 8) | bytes[index++]);
W = (UInt32)((bytes[index++] << 24) | (bytes[index++] << 16) | (bytes[index++] << 8) | bytes[index++]);
}
public Bernoulli(UInt32 x, UInt32 y, UInt32 z, UInt32 w)
{
X = x;
Y = y;
Z = z;
W = w;
}
UInt64 bits = 0;
UInt32 bitsCount = 0;
public UInt32 Next(UInt32 k, UInt32 n)
{
UInt32 b;
var c = 0;
var v = k;
var r = 0u;
// ------------------------
do
{
while (bitsCount <= 32)
{
b = X ^ (X << 11);
X = Y;
Y = Z;
Z = W;
bits <<= 32;
bits |= ((UInt64)(W = W ^ (W >> 19) ^ (b ^ (b >> 8))));
bitsCount += 32;
}
while (c < 32 && 0 < bitsCount)
{
v <<= 1;
// Two lines of code below is a two step optimization:
// First we optimize the following statement:
//
// if (v >= n)
// {
// v -= n;
// b = 1;
// }
// else
// {
// b = 0;
// }
//
// into the following:
//
// var b = v < n ? 0u : 1u;
// v -= b * n
//
// thus reducing branching, but we would like also to omit
// multiplication, which we can do through:
b = v < n ? 0u : 0xFFFFFFFFu;
v -= b & n;
if ((bits & 1) == 1)
{
r |= b & 1;
r <<= 1;
v = k;
c++;
}
bits >>= 1;
bitsCount--;
}
} while (c < 32);
return r;
}
}