o90*_*000 57 c c++ random algorithm bit-manipulation
通常,随机数发生器返回比特流,对于该比特流,在每个位置观察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比特流中设置为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
这种实现比生成无偏流的实现要慢得多,因为它每个位调用一次随机数生成器函数; 而无偏流生成器每字大小调用一次(例如,它可以通过一次调用生成32或64个随机位).
我想要更快的实现,即使它稍微牺牲了随机性.想到的一个想法是预先计算查找表:对于p的200个可能值中的每一个,使用较慢的算法计算C 8位值并将它们保存在表中.然后快速算法将随机选择其中一个以生成8个偏斜位.
包络计算的背面,以查看需要多少内存:C应至少为256(可能的8位值的数量),可能更多以避免采样效果; 比方说,1024也许数量应具体取决于p各不相同,但让我们保持它的简单,说平均为1024由于是p的200个值=>总内存为200 KB.这不错,可能适合L2缓存(256 KB).我仍然需要对它进行评估,以确定是否存在引入偏差的采样效应,在这种情况下,必须增加C.
该解决方案的不足之处是,它可以一次生成只有8位,甚至有大量的工作,而一个不带偏见PRNG可以生成64只有几个算术指令一次.
我想知道是否有一种更快的方法,基于位操作而不是查找表.例如,直接修改随机数生成代码以为每个位引入偏差.这将实现与无偏PRNG相同的性能.
谢谢大家的建议,我收到了很多有趣的想法和建议.以下是最重要的:
注意:正如你们许多人建议的那样,我将分辨率从1/200更改为1/256.
我写了几个朴素方法的实现,只需要8个随机无偏位并生成1个偏置位:
我使用两个无偏的伪随机数生成器:
我还测量了无偏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.
简单的测试表明算术解码代码是瓶颈,而不是PRNG.所以我只是对最昂贵的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倍.可能存在进一步优化和/或并行化算术解码方法的方法,但由于其复杂性,我决定停在这里并选择朴素的SIMD实现.
谢谢大家的帮助.
min*_*iot 28
您可以做的一件事是从底层无偏生成器中多次采样,得到几个32位或64位字,然后执行按位布尔算法.例如,对于4个单词b1,b2,b3,b4,您可以获得以下分布:
expression | p(bit is 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%
可以针对更精细的分辨率进行类似的构造.它有点乏味,仍然需要更多的发电机调用,但每位至少不需要一个.这类似于a3f的答案,但可能更容易实现,而且我怀疑,比扫描0xFnybbles的单词更快.
请注意,对于您所需的0.5%分辨率,您需要8个无偏字,用于一个偏置字,这样可以得到(0.5 ^ 8)= 0.390625%的分辨率.
ric*_*ici 25
如果您准备p基于256个可能的值进行近似,并且您有一个PRNG可以生成单个位彼此独立的统一值,那么您可以使用向量化比较从单个随机数生成多个偏置位.
如果(1)您担心随机数质量并且(2)您可能需要具有相同偏差的大量位,那么这是值得做的.第二个要求似乎暗示了原始问题,它批评了一个提出的解决方案,如下:"这个解决方案的一个缺点是它一次只能生成8位,即使有很多工作,而无偏见的PRNG只需几条算术指令即可一次生成64个." 这里的暗示似乎是在单个调用中生成大块偏置位是有用的.
随机数质量是一个难题.即使不是不可能测量也很难,因此不同的人会提出强调和/或贬低"随机性"不同方面的不同指标.通常可以权衡随机数生成的速度以降低"质量"; 这是否值得做,取决于您的确切应用.
随机数质量的最简单的可能测试涉及单个值的分布和发电机的周期长度.C库rand和Posix random函数的标准实现通常会通过分布测试,但循环长度不适合长时间运行的应用程序.
然而,这些发生器通常非常快:glibc实现random仅需要几个周期,而经典线性同余发生器(LCG)需要乘法和加法.(或者,在glibc实现的情况下,上面的三个产生31位.)如果这对于您的质量要求是足够的,那么尝试优化没有什么意义,特别是如果偏差概率经常变化.
请记住,循环长度应比预期的样本数长很多; 理想情况下,它应该大于该数字的平方,因此如果您希望生成千兆字节的随机数据,则周期长度为2 31的线性同余生成器(LCG)是不合适的.即使是Gnu三项非线性附加反馈发生器,其周期长度据称约为2 35,也不应用于需要数百万个样品的应用中.
另一个质量问题,更难以测试,与连续样本的独立性有关.短周期长度完全失败,因为一旦重复开始,生成的随机数就会与历史值精确相关.Gnu三项式算法虽然周期较长,但由于产生的第i 个随机数r i始终是两个值r i -3 + r i -31或r之一,因此具有明确的相关性.我 -3 + r i -31 +1.这可能会产生令人惊讶或至少令人费解的后果,尤其是伯努利实验.
这是一个使用Agner Fog有用的矢量类库的实现,它抽象了SSE内在函数中的许多恼人细节,并且还有一个快速矢量化随机数生成器(在存档special.zip内部找到vectorclass.zip),它允许我们从中生成256位八次调用256位PRNG.你可以阅读福格博士的解释,为什么他甚至发现梅森捻线机有质量问题,以及他提出的解决方案; 我没有资格发表评论,但是至少它似乎在我尝试过的伯努利实验中给出了预期的结果.
#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;
};
在我的测试中,它在260毫秒或每纳秒一位产生并计数268435456位.测试机是i5,所以它没有AVX2; 因人而异.
在实际使用情况中,对于201个可能的值p,8位阈值的计算将令人烦恼地不精确.如果不期望这种不精确,您可以调整上述内容以使用16位阈值,但代价是产生两倍的随机数.
或者,您可以手动滚动基于10位阈值的矢量化,这将使您获得0.5%增量的非常好的近似值,使用通过检查每10位借位进行矢量化阈值比较的标准位操作黑客减去值向量和重复阈值.例如,结合std::mt19937_64每个64位随机数平均为6位.
Mar*_*son 17
从一个角度信息理论点,(具有比特的偏置流p != 0.5)具有小于在它比无偏流的信息,因此在理论上它应当采取(平均)小于无偏输入的大于1个比特以产生单一比特有偏见的输出流.例如,伯努利随机变量的熵p = 0.1是-0.1 * log2(0.1) - 0.9 * log2(0.9)位,它位于0.469比特周围.这表明,对于这种情况,p = 0.1我们应该能够为每个无偏输入位产生略高于两位的输出流.
下面,我给出了两种产生偏置位的方法.在需要尽可能少的输入无偏位的意义上,两者都实现接近最佳效率.
一种实用的方法是使用算术(de)编码解码您的无偏输入流,如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
正如所承诺的那样,我们已经从源无偏流中使用少于五十万的输出偏置流产生了100万比特.
出于优化目的,在将其转换为C/C++时,使用基于整数的定点算法而不是浮点算法对其进行编码可能是有意义的.
这是一种更简单的方法,而不是试图将算术解码方法转换为直接使用整数.它不再是算术解码,但它并非完全不相关,并且它实现了与上面的浮点版本相同的输出偏置位/输入无偏比特率.它的组织使得所有数量都适合无符号的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循环的开始时,value均匀地分布在所有整数中[0, high),并且独立于所有先前输出的位.如果你关心速度而不是完美的正确性,你可以摆脱discard和value >= discard分支:这就是确保我们输出0并1具有恰当的正确概率.把这个复杂性留下来,你将获得几乎正确的概率.此外,如果您将分辨率设置为p等于1/256而不是1/200,则可能会使用位操作替换可能耗时的除法和模运算.
使用与以前相同的测试代码,但使用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.
现在,只需像以前一样生成一个随机数,然后用1111a 替换每个半字节边界,1并将其他所有内容转换为a 0.
根据其他概率进行相应调整.
您将获得理论上的最佳行为,即,如果您使用算术编码接近此值,则可以对随机数生成器进行真正最小的使用,并能够精确地对任何概率p进行建模.
算术编码是一种数据压缩形式,它将消息表示为数字范围的子间隔.它提供理论上的最佳编码,并且可以为每个输入符号使用小数位.
我们的想法是这样的:假设您有随机位,这是序列1的概率为p.为方便起见,我将使用q该位为零的概率.(q = 1-p).算术编码分配给数字范围的每个比特部分.对于第一位,如果输入为0则分配间隔[0,q),如果输入为1则分配间隔[q,1].后续位分配当前范围的比例子间隔.例如,假设q = 1/3输入1 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
第一个数字1选择范围的前三分之二(1-q); 第二位数字,0,选择的底部三分之一即,依此类推.在第一步和第二步之后,间隔跨越中点; 但是在第三步之后它完全低于中点,所以可以输出第一个压缩数字:0.该过程继续,并EOF添加一个特殊符号作为终止符.
这与你的问题有什么关系？压缩输出将具有随机零和具有相等概率的输出.因此,为了获得具有概率p的比特,只需假装您的RNG的输出是如上所述的算术编码的结果,并将解码器处理应用于它.也就是说,读取位就好像它们将行间隔细分为更小和更小的片段一样.例如,在我们01从RNG 读取之后,我们将在[0.25,0.5]范围内.继续读取位,直到"解码"足够的输出.由于你正在模仿解压缩,你会得到比你输入更多的随机比特.因为算术编码在理论上是最优的,所以没有可能的方法将RNG输出变成更偏向的比特而不牺牲随机性:你得到的是真的最大值.
问题是你不能在几行代码中做到这一点,我不知道我可以指向你的库(虽然必须有一些你可以使用).不过,这很简单.在上面的文章提供了一种通用的编码器和解码器的代码,在C.这是相当简单的,它支持具有任意概率多个输入符号; 在你的情况下,一个更简单的实现是可能的(正如Mark Dickinson的答案现在所示),因为概率模型是微不足道的.为了扩展使用,需要做更多的工作来产生一个不为每个位做很多浮点计算的健壮实现.
维基百科也有一个有趣的讨论算术编码被视为基数的变化,这是查看你的任务的另一种方式.
小智 8
呃,伪随机数发生器一般都很快.我不确定这是什么语言(也许是Python),但是"result.append"(几乎肯定包含内存分配)可能比"random_uniform"慢(这只是做一些数学运算).
如果要优化此代码的性能:
如果您使用的是编译语言(甚至JIT编译),则每次控制转移(如果,同时,函数调用等)都会受到性能影响.消除你所能做的.内存分配也(通常)非常昂贵.
如果您使用的是解释性语言,则所有投注都将被取消.最简单的代码很可能是最好的.解释器的开销会使你正在做的事情变得相形见绌,所以尽可能减少它的工作量.
我只能猜到你的性能问题在哪里:
顺便说一下,如果你真的想尽可能使用最少的随机性,可以使用"算术编码"来解码随机流.它不会很快.
如果p接近0,则可以计算第n位为1的第一位的概率; 然后你计算一个0到1之间的随机数,并相应地选择n.例如,如果p = 0.005(0.5%),并且随机数是0.638128,您可以计算(我猜这里)n = 321,因此您填充321 0位和一位设置.
如果p接近1,则使用1-p而不是p,并设置1位加1位0.
如果p不接近1或0,则制作一个包含所有256个8位序列的表,计算它们的累积概率,然后得到一个随机数,在累积概率数组中进行二分搜索,你可以设置8位.
假设您可以访问随机位的生成器,您可以生成一个值以p逐位进行比较,并在您证明生成的值小于或大于或等于时中止p.
按以下步骤在给定概率的流中创建一个项目p:
0.二进制开始1已被绘制,你会得到0.1p,则输出a1p,则输出a0我们假设p二进制表示法是0.1001101...; 如果这个过程产生的任何的0.0,0.1000,0.10010,...,值不能变得大于或等于p了; 如果有的话0.11,0.101,0.100111,...生成,该值不能成为小于p.
对我来说,看起来这个方法在期望中使用了大约两个随机位.算术编码(如Mark Dickinson的答案所示)每个偏置位(平均)最多消耗一个随机位用于固定p; 修改的成本p尚不清楚.
该实现通过"/ dev/urandom"特殊字符文件的接口对随机设备内核模块进行单次调用,以获得表示给定分辨率中所有值所需的随机数据的数量.最大可能分辨率为1/256 ^ 2,因此0.005可以表示为:
256分之328^ 2,
即:
分辨率:256*256
x:328
错误0.000004883.
该实现计算位数bits_per_byte,该位数是处理给定分辨率所需的均匀分布的位数,即表示所有@resolution值.它然后对随机化设备进行单次调用(如果URANDOM_DEVICE定义了"/ dev/urandom" ,否则它将使用来自设备驱动程序的额外噪声通过调用"/ dev/random",如果位中没有足够的熵,则可能会阻塞)获得所需数量的均匀分布字节并填充rnd_bytes随机字节数组.最后,它从rnd_bytes数组的每个bytes_per_byte字节读取每个伯努利样本所需的比特数,并将这些比特的整数值与单个伯努利结果的成功概率进行比较x/resolution.如果值命中,即它落在x/resolution我们任意选择为[0,x /分辨率] 段的长度段中,那么我们注意到成功并将1插入到结果数组中.
从随机设备中读取:
/* 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;
}
填写伯努利样本的矢量:
/*
* @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;
}