快速生成每个比特具有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实现。

谢谢大家的帮助。

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

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

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

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

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

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

另一个质量问题难以testing,涉及连续样品的独立性。 短周期长度在此度量标准上完全失败,因为一旦重复开始,生成的随机数与历史值精确相关。 Gnu三项algorithm虽然周期较长,但由于第i 随机数r i总是r i -3 + r i -31r 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位。

你可以做的一件事是从底层的无偏生成器多次采样,得到几个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%的分辨率。

从信息论的观点来看,有偏位的比特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 ++。 该代码专门用于p1/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循环的开始whilevalue[0, high)中的所有整数中均匀分布,并且独立于所有先前输出的比特。 如果你关心速度而不是完美的正确性,你可以摆脱discardvalue >= discard分支:这只是为了确保我们输出01正确的概率。 把这个并发症抛出去,你就会得到几乎正确的概率。 此外,如果将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 

假设1出现的概率是6.25%(1/16)。 4位数有16种可能的位模式: 0000,0001, ..., 1110,1111

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

相应地调整其他概率。

你会得到理论上的最佳行为,即使真正最小限度的使用随机数发生器,并能够正确地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符号。 在你的情况下,一个更简单的实现是可能的(正如马克·迪金森的答案现在所示),因为概率模型是微不足道的。 对于扩展的使用,需要做更多的工作来产生一个健壮的实现,而不需要为每一位做大量的浮点计算。

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

呃,伪随机数发生器通常是相当快的。 我不确定这是什么语言(可能是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。 它不会很快。

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

如果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位。

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

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

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

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

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

它能做什么

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 .