【发布时间】：2020-03-21 19:34:10
【问题描述】：
假设您有一个常规数字生成器，它能够生成均匀分布的随机 32 位数字。假设您正在寻找一种方法来生成伪随机位序列，其中位（即设置位）以预定义的概率出现在序列中。
产生这种序列的一种简单方法是在每个位级别上运行数字生成器，但是对于序列中大多数位为零的小概率（例如 0.01 或 1% 的位）来说，这是非常低效的。平均而言，只有百分之一会被设置。另一方面，即使概率如此之低，也有机会遇到超过 8、16、32、64 位的连续子序列的长子序列。
问题是如何使用常规 PRNG 有效地生成这样的序列。
编辑
Peter O. 建议的 javascript 中理性伯努利变量采样的玩具实现：
// 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# 实现，每秒生成 9120 万位，打包成 UInt32 数据类型（MacBookPro 2019 Core I9 2.4 Ghz）。我认为在 C 中可能会获得超过 1 亿位，而且感觉可以进一步利用二进制算术并行生成所有 32 位随机数，一些循环展开或者 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;
}
}
【问题讨论】：
标签： random pseudocode