快速生成每个比特具有0或1的概率的伪随机比特的方法

通常，一个随机数发生器返回一个位stream，在每个位置观察0或1的概率是相等的（即50％）。 我们称之为一个公正的PRNG。

我需要生成一串具有以下属性的伪随机比特：在每个位置上看到1的概率是p（即看到0的概率是1-p）。 参数p是0到1之间的实数; 在我的问题中它恰好有0.5％的分辨率，即它可以取值0％，0.5％，1％，1.5％，…，99.5％，100％。

请注意，p是一个概率，而不是一个确切的分数。 在n位数据stream中设置为1的实际位数必须遵循二项式分布B（n，p）。

有一种天真的方法，可以使用无偏PRNG来生成每个位的值（伪代码）：

generate_biased_stream(n, p): result = [] for i in 1 to n: if random_uniform(0, 1) < p: result.append(1) else: result.append(0) return result 

这样的实现比产生无偏stream的实现要慢得多，因为它每位调用一次随机数生成器函数; 而无偏stream生成器每个字大小调用一次（例如，一次调用就可以产生32或64个随机位）。

我想要一个更快的实现，即使它稍微牺牲随机性。 想到一个想法是预先计算一个查找表：对于p的200个可能的值中的每一个，使用较慢的algorithm计算C 8位值并将它们保存在表中。 然后快速algorithm将随机select其中的一个来产生8个偏斜位。

在包络计算的后面，看看需要多less内存：C应该至less为256（可能的8位值的数量），可能更多是为了避免采样效应; 假设1024.也许这个数字应该取决于p，但是让我们保持简单，并说平均值是1024.因为有200个值=>总的内存使用量是200 KB。 这并不坏，可能适合二级caching（256 KB）。 我仍然需要评估它是否存在引入偏差的抽样效应，在这种情况下，C必须增加。

这个解决scheme的不足之处在于，它只能一次产生8位数据，即使有很多工作也是如此，而一个无偏PRNG只需要一些算术指令就可以产生64个数据。

我想知道是否有一个更快的方法，基于位操作，而不是查找表。 例如直接修改随机数生成代码，为每一位引入一个偏差。 这将达到与无偏见的PRNG相同的性能。

---

3月5日编辑

谢谢大家的build议，我有很多有趣的想法和build议。 这里是最重要的：

• 更改问题要求，以便p的分辨率为1/256而不是1/200。 这允许更有效地使用比特，并且还提供了更多优化的机会。 我想我可以做这个改变。
• 使用算术编码来高效地消耗来自无偏置发生器的比特。 随着上述分辨率的变化，这变得更容易。
• 有less数人认为PRNG速度非常快，因此使用算术编码可能会导致代码变慢，这是由于引入的开销。 相反，我应该总是消耗最坏情况的位数并优化代码。 看下面的基准。
• @ricibuild议使用SIMD。 这是一个好主意，只有在我们总是消耗固定的比特数的情况下才有效。

基准（无算术解码）

注意：正如你们许多人所build议的那样，我将分辨率从1/200更改为1/256。

我写了几个朴素方法的实现 ，它只需要8个随机无偏位，并产生1个有偏位：

• 没有SIMD
• SIMD使用Agner Fog的vector类库， 如@rici所示
• SIMD使用内部函数

我使用两个无偏伪随机数发生器：

• xorshift128plus
• 来自Agner Fog的图书馆的Ranvec1（Mersenne Twister-like） 。

我还测量了无偏PRNG的速度进行比较。 结果如下：

 RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry) Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 16.081 16.125 16.093 [Gb/s] Number of ones: 536,875,204 536,875,204 536,875,204 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 0.778 0.783 0.812 [Gb/s] Number of ones: 104,867,269 104,867,269 104,867,269 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 2.176 2.184 2.145 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 2.129 2.151 2.183 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600 

与标量方法相比，SIMD将性能提高了3倍。 正如所料，它比无偏差发电机慢了8倍。

最快的偏置发生器达到2.1 Gb / s。

 RNG: xorshift128plus Method: Unbiased with 1/1 efficiency (incorrect, baseline) Gbps/s: 18.300 21.486 21.483 [Gb/s] Number of ones: 536,867,655 536,867,655 536,867,655 Theoretical : 104,857,600 Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 22.660 22.661 24.662 [Gb/s] Number of ones: 536,867,655 536,867,655 536,867,655 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 1.065 1.102 1.078 [Gb/s] Number of ones: 104,868,930 104,868,930 104,868,930 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 4.972 4.971 4.970 [Gb/s] Number of ones: 104,869,407 104,869,407 104,869,407 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 4.955 4.971 4.971 [Gb/s] Number of ones: 104,869,407 104,869,407 104,869,407 Theoretical : 104,857,600 

对于xorshift，与标量方法相比，SIMD将性能提高了5倍。 它比无偏差发生器慢4倍。 请注意，这是xorshift的标量实现。

最快的偏置发生器达到了4.9 Gb / s。

 RNG: xorshift128plus_avx2 Method: Unbiased with 1/1 efficiency (incorrect, baseline) Gbps/s: 18.754 21.494 21.878 [Gb/s] Number of ones: 536,867,655 536,867,655 536,867,655 Theoretical : 104,857,600 Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 54.126 54.071 54.145 [Gb/s] Number of ones: 536,874,540 536,880,718 536,891,316 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 1.093 1.103 1.063 [Gb/s] Number of ones: 104,868,930 104,868,930 104,868,930 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 19.567 19.578 19.555 [Gb/s] Number of ones: 104,836,115 104,846,215 104,835,129 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 19.551 19.589 19.557 [Gb/s] Number of ones: 104,831,396 104,837,429 104,851,100 Theoretical : 104,857,600 

这个实现使用AVX2并行运行4个无偏差的xorshift生成器。

最快的偏置发生器达到19.5 Gb / s。

算术解码的基准

简单的testing表明，算术解码代码是瓶颈，而不是PRNG。 所以我只是testing最贵的PRNG。

 RNG: Ranvec1(Mersenne Twister for Graphics Processors + Multiply with Carry) Method: Arithmetic decoding (floating point) Gbps/s: 0.068 0.068 0.069 [Gb/s] Number of ones: 10,235,580 10,235,580 10,235,580 Theoretical : 10,240,000 Method: Arithmetic decoding (fixed point) Gbps/s: 0.263 0.263 0.263 [Gb/s] Number of ones: 10,239,367 10,239,367 10,239,367 Theoretical : 10,240,000 Method: Unbiased with 1/1 efficiency (incorrect, baseline) Gbps/s: 12.687 12.686 12.684 [Gb/s] Number of ones: 536,875,204 536,875,204 536,875,204 Theoretical : 104,857,600 Method: Unbiased with 1/1 efficiency, SIMD=vectorclass (incorrect, baseline) Gbps/s: 14.536 14.536 14.536 [Gb/s] Number of ones: 536,875,204 536,875,204 536,875,204 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency Gbps/s: 0.754 0.754 0.754 [Gb/s] Number of ones: 104,867,269 104,867,269 104,867,269 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=vectorclass Gbps/s: 2.094 2.095 2.094 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600 Method: Biased with 1/8 efficiency, SIMD=intrinsics Gbps/s: 2.094 2.094 2.095 [Gb/s] Number of ones: 104,859,067 104,859,067 104,859,067 Theoretical : 104,857,600 

简单的不动点方法达到了0.25 Gb / s，而朴素标量法快了3倍，而简单的SIMD方法快了8倍。 可能有办法进一步优化和/或并行化算术解码方法，但是由于其复杂性，我决定在这里停下来select简单的SIMD实现。

谢谢大家的帮助。

Answer 1

如果您准备根据256个可能的值来近似p ，并且您有一个PRNG可以生成统一的值，其中各个位彼此独立，那么可以使用向量化比较来从单个随机生成多个有偏位数。

如果（1）您担心随机数字的质量，以及（2）您可能需要相同偏差的大量比特，那么这样做是值得的。 第二个要求似乎被原来的问题所暗示，它批评了一个提议的解决scheme，如下所示：“这个解决scheme的不足之处在于，它只能一次产生8位数据，即使有很多工作，而没有偏见的PRNG可以用一些算术指令一次生成64个。“ 在这里，暗示似乎是在一次调用中产生一大块有偏位的是有用的 。

随机数质量是一个困难的课题。 如果不是不可能的话，很难衡量，因此不同的人会提出不同的指标，强调和/或贬低“随机性”的不同方面。 一般来说，为了降低“质量”，可以权衡随机数生成的速度; 这是否值得做，取决于你的确切应用。

随机数质量的最简单可能的testing涉及个体值的分布和发生器的周期长度。 C库rand和Posix random函数的标准实现通常会通过分配testing，但是对于长时间运行的应用程序来说，周期长度是不够的。

然而，这些发生器通常非常快速： random的glibc实现只需要几个周期，而经典的线性同余发生器（LCG）需要一个乘法和一个加法。 （或者，在glibc实现的情况下，上述三个生成31位）。如果这足以满足您的质量要求，那么尝试优化是没有意义的，特别是在偏差概率频繁变化的情况下。

请记住，周期长度应该比预期的样本数量长很多; 理想情况下，它应该大于那个数的平方，所以如果你期望产生千兆字节的随机数据，那么周期长度为2 ³¹的线性同余发生器（LCG）是不合适的。 即使是Gnu三项式非线性加性反馈发生器，其周期长度约为2 ³⁵ ，不应用于需要数百万个样本的应用。

另一个质量问题难以testing，涉及连续样品的独立性。 短周期长度在此度量标准上完全失败，因为一旦重复开始，生成的随机数与历史值精确相关。 Gnu三项algorithm虽然周期较长，但由于第i ^个随机数r _i总是r _{i -3} + r _{i -31}或r _{i -3} + _{我 – 31} ＆1; 1。 这可能会令人惊讶或至less令人费解的后果，特别是伯努利实验。

下面是一个使用Agner Fog的有用的向量类库的实现 ，它抽象出SSE内在函数中很多烦人的细节，还有一个快速向量化的随机数生成器 （可在vectorclass.zip归档文件的vectorclass.zipfind）让我们从8个调用产生256位到256位PRNG。 您可以阅读Fog博士解释为什么他甚至find了Mersenne扭曲者质量问题，以及他提出的解决scheme; 我没有资格发表评论，但它确实至less在我尝试过的伯努利实验中给出了预期的结果。

 #include "vectorclass/vectorclass.h" #include "vectorclass/ranvec1.h" class BiasedBits { public: // Default constructor, seeded with fixed values BiasedBits() : BiasedBits(1) {} // Seed with a single seed; other possibilities exist. BiasedBits(int seed) : rng(3) { rng.init(seed); } // Generate 256 random bits, each with probability `p/256` of being 1. Vec8ui random256(unsigned p) { if (p >= 256) return Vec8ui{ 0xFFFFFFFF }; Vec32c output{ 0 }; Vec32c threshold{ 127 - p }; for (int i = 0; i < 8; ++i) { output += output; output -= Vec32c(Vec32c(rng.uniform256()) > threshold); } return Vec8ui(output); } private: Ranvec1 rng; }; 

在我的testing中，在260毫秒内产生并计数268435456位，或者每纳秒一位。 testing机是i5，所以没有AVX2; 因人而异。

在实际使用情况下，对于p有201个可能的值，8位阈值的计算将令人生厌地不准确。 如果这种不精确性是不希望的，那么可以使用上面的16位阈值，代价是产生两倍的随机数。

或者，您可以手动滚动基于10位阈值的向量化，这会使您以0.5％的增量给出一个非常好的近似值，使用标准的位操作黑客做vector化的阈值比较，通过检查每10位减去向量的值和重复的阈值。 结合，比如std::mt19937_64 ，会给你一个平均每个64位随机数的6位。

Answer 2

你可以做的一件事是从底层的无偏生成器多次采样，得到几个32位或64位的字，然后执行按位布尔运算。 例如，对于4个单词b1,b2,b3,b4 ，您可以得到以下分布：

    expression式|  p（位是1）
     ----------------------- + -------------
     b1＆b2＆b3＆b4 |  6.25％
     b1＆b2＆b3 |  12.50％
     b1＆b2＆（b3 | b4）|  18.75％
     b1＆b2 |  25.00％
     b1 |  （b2＆（b3 | b4））|  31.25％
     b1＆（b2 | b3）|  37.50％
     b1＆（b2 | b3 | b4））|  43.75％
     b1 |  50.00％

类似的结构可以做出更好的分辨率。 它有点繁琐，仍然需要更多的发生器调用，但至less不是每位一个。 这与a3f的答案类似，但可能更容易实现，而且我怀疑它比扫描0xF nybbles的单词更快。

请注意，对于你想要的0.5％的分辨率，你需要8个无偏的单词来处理一个有偏见的单词，这会给你一个（0.5 ^ 8）= 0.390625％的分辨率。

Answer 3

从信息论的观点来看，有偏位的比特stream（ p != 0.5 ）比无偏stream有更less的信息，所以理论上它应该（平均） less于 1比特的无偏input产生一个偏置输出stream的位。 例如， p = 0.1的伯努利随机variables的熵是-0.1 * log2(0.1) - 0.9 * log2(0.9)位，大约为0.469位。 这表明，对于p = 0.1的情况，我们应该能够在每个无偏input位上产生略多于两位的输出stream。

下面我给出两种产生偏置位的方法。 从需要尽可能less的input无偏位的angular度来看，两者都达到接近最佳效率。

方法1：算术（解）编码

一个实用的方法是使用算术（de）编码来解码您的无偏inputstream，正如在alexis的答案中已经描述的那样。 对于这种简单的情况，编码起来并不难。 这里有一些未经优化的伪代码（ 咳嗽，Python ）：

 import random def random_bits(): """ Infinite generator generating a stream of random bits, with 0 and 1 having equal probability. """ global bit_count # keep track of how many bits were produced while True: bit_count += 1 yield random.choice([0, 1]) def bernoulli(p): """ Infinite generator generating 1-bits with probability p and 0-bits with probability 1 - p. """ bits = random_bits() low, high = 0.0, 1.0 while True: if high <= p: # Generate 1, rescale to map [0, p) to [0, 1) yield 1 low, high = low / p, high / p elif low >= p: # Generate 0, rescale to map [p, 1) to [0, 1) yield 0 low, high = (low - p) / (1 - p), (high - p) / (1 - p) else: # Use the next random bit to halve the current interval. mid = 0.5 * (low + high) if next(bits): low = mid else: high = mid 

这是一个示例用法：

 import itertools bit_count = 0 # Generate a million deviates. results = list(itertools.islice(bernoulli(0.1), 10**6)) print("First 50:", ''.join(map(str, results[:50]))) print("Biased bits generated:", len(results)) print("Unbiased bits used:", bit_count) print("mean:", sum(results) / len(results)) 

以上给出了以下示例输出：

 First 50: 00000000000001000000000110010000001000000100010000 Biased bits generated: 1000000 Unbiased bits used: 469036 mean: 0.100012 

如所承诺的，我们已经使用来自源无偏stream的不到五十万个产生了100万比特的输出偏置stream。

为了进行优化，在将其转换为C / C ++时，使用基于整数的定点algorithm而不是浮点代码可能是有意义的。

方法2：基于整数的algorithm

而不是试图将算术解码方法直接使用整数，这是一个更简单的方法。 它不再是算术解码，但它不完全不相关，并且达到了与上述浮点版本接近相同的输出偏置位/input无偏位比。 它的组织结构使得所有的数量都适合一个无符号的32位整数，所以应该很容易转换成C / C ++。 该代码专门用于p是1/200的精确倍数的情况，但是这种方法适用于任何可以表示为具有合理的小分母的有理数的p 。

 def bernoulli_int(p): """ Infinite generator generating 1-bits with probability p and 0-bits with probability 1 - p. p should be an integer multiple of 1/200. """ bits = random_bits() # Assuming that p has a resolution of 0.05, find p / 0.05. p_int = int(round(200*p)) value, high = 0, 1 while True: if high < 2**31: high = 2 * high value = 2 * value + next(bits) else: # Throw out everything beyond the last multiple of 200, to # avoid introducing a bias. discard = high - high % 200 split = high // 200 * p_int if value >= discard: # rarer than 1 time in 10 million value -= discard high -= discard elif value >= split: yield 0 value -= split high = discard - split else: yield 1 high = split 

关键的观察是，每当我们到达while循环的开始while ， value在[0, high)中的所有整数中均匀分布，并且独立于所有先前输出的比特。 如果你关心速度而不是完美的正确性，你可以摆脱discard和value >= discard分支：这只是为了确保我们输出0和1正确的概率。 把这个并发症抛出去，你就会得到几乎正确的概率。 此外，如果将p的分辨率设置为1/256而不是1/200 ，那么可能会耗费时间的分割和模运算可以用位运算来replace。

使用与之前相同的testing代码，但是使用bernoulli_int代替bernoulli ，得到以下结果p=0.1 ：

 First 50: 00000010000000000100000000000000000000000110000100 Biased bits generated: 1000000 Unbiased bits used: 467997 mean: 0.099675 

Answer 4

假设1出现的概率是6.25％（1/16）。 4位数有16种可能的位模式： 0000,0001, ..., 1110,1111 。

现在，只需要像以前那样生成一个随机数，并用1代替半个边界上的每个1111 ，然后将其他所有内容都replace为0 。

相应地调整其他概率。

Answer 5

你会得到理论上的最佳行为，即使真正最小限度的使用随机数发生器，并能够正确地build模任何概率p ，如果你使用算术编码来解决这个问题。

算术编码是数据压缩的一种forms，它将消息表示为数字范围的子间隔。 它提供了理论上的最佳编码，并且可以使用每个input符号的小数位数。

这个想法是这样的：假设你有一个随机位序列，它的概率是1 。 为了方便起见，我将使用q代表位为零的概率。 （ q = 1-p ）。 算术编码分配给数字范围的每个位部分。 对于第一位，如果input为0，则分配区间[0，q）;如果input为1，则分配区间[q，1）。随后的位分配当前范围的比例子区间。 例如，假设q = 1/3input1 0 0将被编码如下：

 Initially [0, 1), range = 1 After 1 [0.333, 1), range = 0.6666 After 0 [0.333, 0.5555), range = 0.2222 After 0 [0.333, 0.407407), range = 0.074074 

第一个数字1select范围的前三分之二（1-q） 第二个数字0select下面的三分之一，依此类推。 第一步和第二步之后，间隔横跨中点; 但是在第三步之后它完全低于中点，所以可以输出第一个压缩数字： 0 。 该过程继续，并添加一个特殊的EOF符号作为终结符。

这与你的问题有什么关系？ 压缩的输出将具有随机的零和具有相等概率的零。 因此，要获得概率为p的比特，只需假设RNG的输出是如上所述的算术编码的结果，并对其应用解码器处理。 也就是说，读取位就像将行间隔细分为更小和更小的片段一样。 例如，从RNG读取01之后，我们将在[0.25,0.5]的范围内。 保持读取位，直到足够的输出被“解码”。 由于你正在模拟解压缩，所以你会得到比你放入的更多的随机数。因为算术编码在理论上是最优的，所以没有任何可能的方法把RNG输出变成更偏向的位而不牺牲随机性：最大。

问题是，你不能用几行代码来做这件事，而且我不知道我可以指向你的一个库（尽pipe必须有一些你可以使用）。 不过，这很简单。 上面的文章以C语言为通用编码器和解码器提供了代码。它非常简单，它支持任意概率的多个input符号。 在你的情况下，一个更简单的实现是可能的（正如马克·迪金森的答案现在所示），因为概率模型是微不足道的。 对于扩展的使用，需要做更多的工作来产生一个健壮的实现，而不需要为每一位做大量的浮点计算。

维基百科也有一个有趣的讨论算术编码被认为是基数的变化，这是另一种方式来查看您的任务。

Answer 6

呃，伪随机数发生器通常是相当快的。 我不确定这是什么语言（可能是Python），但是“result.append”（几乎可以肯定包含内存分配）可能比“random_uniform”（这只是一个小math）慢。

如果你想优化这个代码的性能：

1. 确认这是一个问题。 优化是一个工作，使代码更难以维护。 除非必要，否则不要这样做。
2. configuration文件。 运行一些testing来确定哪些代码部分实际上是最慢的。 这些是你需要加快的部分。
3. 进行更改，并确认其实际速度更快。 编译器非常聪明; 通常清晰的代码会被编译成更好的代码，而这些代码可能比看上去更复杂。

如果你使用的是编译语言（甚至是JIT编译），那么每一次控制转移（如果，函数调用等）都会带来性能上的冲击。 消除你可以。 内存分配也是（通常）非常昂贵。

如果您使用的是解释性语言，则所有投注均为closures。 最简单的代码很可能是最好的。 翻译的开销会使你所做的任何事情都变得渺小，所以尽可能地减less它的工作量。

我只能猜测你的性能问题在哪里：

1. 内存分配。 以全尺寸预先分配数组，并稍后填写条目。 这可以确保在添加条目时不需要重新分配内存。
2. 分行。 你可以通过投射结果或类似的东西来避免“如果”。 这将取决于编译器。 检查程序集（或configuration文件），以validation它是否做你想要的。
3. 数字types。 找出你的随机数发生器原生使用的types，并在那种types中进行算术。 例如，如果生成器自然返回32位无符号整数，则首先将“p”缩放到该范围，然后将其用于比较。

顺便说一句，如果你真的想使用尽可能less的随机性，使用“算术编码”来解码你的随机stream。 它不会很快。

Answer 7

一种精确结果的方法是首先随机产生一个k位块，跟随二项分布的1位数，然后使用这里的一种方法产生一个正好多位的k位字。 例如mic006的方法只需要log k k位的随机数，而我的只需要一个。

Answer 8

如果p接近0，则可以计算出第n位是1的第一位的概率; 那么你计算一个0和1之间的随机数，并相应地selectn。 例如，如果p = 0.005（0.5％），而随机数是0.638128，那么你可以计算（我在这里猜测）n = 321，所以你填充321 0位和一个位设置。

如果p接近于1，则使用1-p而不是p，并设置1位加上一个0位。

如果p不接近于1或0，则build立一个包含8位全部256个序列的表格，计算它们的累积概率，然后得到一个随机数，在累积概率数组中进行二进制search，可以设置8位。

Answer 9

假设您可以访问一个随机位的生成器，您可以生成一个值来逐位比较，一旦certificate生成的值小于或大于或等于p ，就立即中止。

按如下步骤在给定概率p的stream中创build一个项目：

1. 从二进制开始
2. 追加一个随机位; 假设1已经被绘制，你会得到0.1
3. 如果结果（二进制表示法）比p小得多，则输出1
4. 如果结果可能大于或等于p ，则输出0
5. 否则（如果两者都不能排除），继续步骤2。

假设二进制符号p是0.1001101... ; 如果该过程产生0.10010 ，…中的任何一个，则该值不能大于或等于p ; 如果产生0.100111 ，…中的任何一个，则该值不能小于p 。

对我来说，看起来这种方法在预期中使用了大约两个随机比特。 算术编码（如Mark Dickinson的答案所示）对于固定的p ，每个偏置位（平均） 最多消耗一个随机位; 修改p的代价不清楚。

Answer 10

它能做什么

This implementation makes single call to random device kernel module via interface of "/dev/urandom" special character file to get number of random data needed to represent all values in given resolution. Maximum possible resolution is 1/256^2 so that 0.005 can be represented by:

328/256^2,

即：

resolution: 256*256

x: 328

with error 0.000004883.

How it does that

The implementation calculates the number of bits bits_per_byte which is number of uniformly distributed bits needed to handle given resolution, ie represent all @resolution values. It makes then a single call to randomization device ("/dev/urandom" if URANDOM_DEVICE is defined, otherwise it will use additional noise from device drivers via call to "/dev/random" which may block if there is not enough entropy in bits) to get required number of uniformly distributed bytes and fills in array rnd_bytes of random bytes. Finally it reads number of needed bits per each Bernoulli sample from each bytes_per_byte bytes of rnd_bytes array and compares the integer value of these bits to probability of success in single Bernoulli outcome given by x/resolution . If value hits, ie it falls in segment of x/resolution length which we arbitrarily choose to be [0, x/resolution) segment then we note success and insert 1 into resulting array.

---

Read from random device:

 /* if defined use /dev/urandom (will not block), * if not defined use /dev/random (may block)*/ #define URANDOM_DEVICE 1 /* * @brief Read @outlen bytes from random device * to array @out. */ int get_random_samples(char *out, size_t outlen) { ssize_t res; #ifdef URANDOM_DEVICE int fd = open("/dev/urandom", O_RDONLY); if (fd == -1) return -1; res = read(fd, out, outlen); if (res < 0) { close(fd); return -2; } #else size_t read_n; int fd = open("/dev/random", O_RDONLY); if (fd == -1) return -1; read_n = 0; while (read_n < outlen) { res = read(fd, out + read_n, outlen - read_n); if (res < 0) { close(fd); return -3; } read_n += res; } #endif /* URANDOM_DEVICE */ close(fd); return 0; } 

Fill in vector of Bernoulli samples:

 /* * @brief Draw vector of Bernoulli samples. * @details @x and @resolution determines probability * of success in Bernoulli distribution * and accuracy of results: p = x/resolution. * @param resolution: number of segments per sample of output array * as power of 2: max resolution supported is 2^24=16777216 * @param x: determines used probability, x = [0, resolution - 1] * @param n: number of samples in result vector */ int get_bernoulli_samples(char *out, uint32_t n, uint32_t resolution, uint32_t x) { int res; size_t i, j; uint32_t bytes_per_byte, word; unsigned char *rnd_bytes; uint32_t uniform_byte; uint8_t bits_per_byte; if (out == NULL || n == 0 || resolution == 0 || x > (resolution - 1)) return -1; bits_per_byte = log_int(resolution); bytes_per_byte = bits_per_byte / BITS_PER_BYTE + (bits_per_byte % BITS_PER_BYTE ? 1 : 0); rnd_bytes = malloc(n * bytes_per_byte); if (rnd_bytes == NULL) return -2; res = get_random_samples(rnd_bytes, n * bytes_per_byte); if (res < 0) { free(rnd_bytes); return -3; } i = 0; while (i < n) { /* get Bernoulli sample */ /* read byte */ j = 0; word = 0; while (j < bytes_per_byte) { word |= (rnd_bytes[i * bytes_per_byte + j] << (BITS_PER_BYTE * j)); ++j; } uniform_byte = word & ((1u << bits_per_byte) - 1); /* decision */ if (uniform_byte < x) out[i] = 1; else out[i] = 0; ++i; } free(rnd_bytes); return 0; } 

用法：

 int main(void) { int res; char c[256]; res = get_bernoulli_samples(c, sizeof(c), 256*256, 328); /* 328/(256^2) = 0.0050 */ if (res < 0) return -1; return 0; } 

Complete code , results .