快速、公正、整数伪随机生成器，具有任意边界。

Question

对于蒙特卡罗积分过程，我需要从具有N个桶的直方图中提取大量随机样本，其中N是任意的（即不是2的幂），但在计算过程中不会改变。

通过很多，我的意思是大约10^10，100亿，因此面对如此庞大的样本数量，任何一种漫长的预处理都可能是值得的。

我可以使用非常快速的均匀伪随机数生成器，通常会产生64位无符号整数（以下讨论中的所有整数都是无符号的）。

拉取样本的朴素方法：histogram [prng()％histogram.size()]

这种朴素方法非常慢：模运算使用整数除法（IDIV），非常昂贵，编译器不知道histogram.size()的值在编译时，不能像平常那样进行优化（即http://www.azillionmonkeys.com/qed/adiv.html）。

事实上，我的大部分计算时间都花在了提取该可恶的模数上。

稍微不那么天真的方法：我使用libdivide（http://libdivide.com/），它能够快速执行“除以编译时不知道的常量”。这给了我一个非常好的优势（约25%），但我有一种困扰的感觉，认为我可以做得更好，原因如下：

- 第一直觉：libdivide计算除法。我需要的是模数，为了到达目的地，我必须进行额外的乘法和减法：mod = dividend - divisor*(uint64_t)(dividend/divisor)。我怀疑可能有一个小胜利，在使用libdivide类型的技术直接生成模数方面。 - 第二直觉：我实际上并不关心模数本身。我真正想要的是有效地生成一个保证严格小于N的均匀分布的整数值。

模数是实现这一目标的一种相当标准的方法，因为它具有两个属性：

• A) 如果prng()满足要求，那么mod(prng()，N)保证是均匀分布的

• B) mod(prgn()，N)保证属于[0，N [

但是取模运算除了满足上述两个限制条件外，还有更多作用，事实上它可能做了太多工作。

我们需要的只是一个函数，任何函数都符合A）和B）的约束条件，并且具有快速性。

所以，长话短说，这里有两个问题：

• 是否存在与libdivide等价的计算整数取模的方法?

• 是否存在整数X和N的函数F（X，N），它遵循以下两个约束条件：

  • 如果X是均匀分布的随机变量，则F(X，N)也是均匀分布的
  • F(X，N)保证在[0，N [内

（PS：我知道如果N很小，我不需要消耗所有从PRNG中输出的64位。实际上，我已经这样做了。但是像我说的那样，即使是这种优化与计算模数相比也只是一个小小的胜利。）

编辑：prng() % N确实不是完全均匀分布的。但对于足够大的N，我认为这不是什么大问题（或者是吗？）

编辑2：prng() % N可能非常差地分布。我从未意识到它会有多糟糕。哎呀。我在这里找到了一篇好文章：http://ericlippert.com/2013/12/16/how-much-bias-is-introduced-by-the-remainder-technique

- blondiepassesby

---

(1) - 在大多数平台上，余数是通过硬件级别的除法计算“免费”得出的结果。 (2) - 模运算和除法都不能给出一个无偏的结果。 - Oliver Charlesworth

只有当N能够整除M时，才能保证 prng() 返回的数字在 [0, M[ 范围内进行取模操作后得到的结果 mod(prng(), N) 均匀分布。请注意，prng() 返回的是均匀分布的随机数。 - huon

你尝试过使用 std::uniform_int_distribution 吗？ - Jarod42

正如@Jarod42所说，什么是典型的N？假设您在“现代”x86上，N > (256kB/8B) = 32k将导致L2缓存溢出，这必定会成为主要的性能影响。 - Oliver Charlesworth

你可以尝试使用 histogram[ (int)(prng() * (HISTOGRAM_SIZE / (PRNG_MAX + 1.0))) ]。每个直方图预先计算常量一次。这将编译为一个浮点乘法和一个整数转换。使用MMX或GPU实现可以同时处理多个。但是，如果您使用随机直方图访问来清除缓存，我同意@OliCharlesworth的看法，这是一个巨大的代价。 - Gene

Answer 1

在这种情况下，最简单的方法可能是最好的。如果您的伪随机数生成器足够快，那么一种极其简单的方法是预先计算比您的N大的下一个2的幂次方减1，然后使用它作为掩码。也就是说，给定一个看起来像0001xxxxxxxx的数字（其中x表示我们不关心它是1还是0），我们需要一个像000111111111这样的掩码。

从那里开始，我们按照以下方式生成数字：

1. 生成一个数字
2. AND它与您的掩码
3. 如果结果> n，则返回第1步

这种方法的确切有效性取决于N与2的幂次方之间的接近程度。每个连续的2的幂次方都是（显然）比其前任加倍。因此，在最佳情况下，N恰好比2的幂次方少1，而我们在步骤3中的测试始终通过。我们只增加了一个掩码和一个比较到PRNG本身所需的时间。

在最坏的情况下，N恰好等于2的幂次方。在这种情况下，我们预计会丢弃大约一半生成的数字。

平均而言，N大约位于2的幂次方之间的中心位置。这意味着平均而言，我们将抛弃四个输入中的一个。我们几乎可以忽略掩码和比较本身，因此与“原始”生成器相比，我们的速度损失基本上等于丢弃的输出数量，或者平均为25%。

- Jerry Coffin

---

我已经实现了你描述的拒绝方法。它的效果相当不错，无论是在质量（均匀性）还是速度方面都优于取模运算。但对于大多数N值来说，即使使用分支预测提示，它仍然比“通过双精度归一化”方法慢。 - blondiepassesby

Answer 2

如果您可以快速访问所需的指令，您可以将prng()与N相乘并返回128位结果的高64位。这有点像将[0,1)中的均匀实数乘以N并截断，偏差大约为模运算版本（即，实际上可以忽略不计；这个答案的32位版本可能有轻微但是可感知的偏差）。

另一个值得探索的可能性是，在单个位上运行无分支模算法的字并行运算，以批量获取随机数。

- David Eisenstat

---

我已经实现了这个。它运行良好，比浮点数方法更快（FP需要20.34秒，而您的方法只需要18.66秒，因此快了约9%）。不错。 - blondiepassesby

Answer 3

Libdivide，或者任何其他复杂的优化取模的方法都是过度的。在您的情况下，唯一明智的方法是：

1. 确保你的表格大小是2的幂次方（如果必须，添加填充！）

2. 用位掩码操作替换模运算。像这样：

   size_t tableSize = 1 << 16;
   size_t tableMask = tableSize - 1;
   
   ...
   
   histogram[prng() & tableMask]
   

位运算是任何一款性能出色的CPU上的单个周期，其速度无与伦比。

--

注意：
我不知道你的随机数生成器的质量如何，但使用随机数的最后几位可能不是一个好主意。一些随机数生成器在最后几位上产生较差的随机性，在高位上产生更好的随机性。如果你的随机数生成器是这种情况，请使用位移来获取最显著的位。

size_t bitCount = 16;

...

histogram[prng() >> (64 - bitCount)]

这种方法和位掩码一样快，但使用的是不同的位。

- cmaster - reinstate monica

---

如果你不需要2的幂次方呢？ - Oliver Charlesworth

@cmaster：25% 的增益并不是我所谓的过度杀伤力。就我而言，这已经相当不错了。至于四舍五入……嗯，我的直方图很大，所以下一个二次幂“还很远”。此外，直方图由来自现实世界的观察值组成，因此你的填充想法 - 至少可以说 - 不是微不足道的：你到底要用什么来“填充”呢？ - blondiepassesby

@OliCharlesworth 我明确提到了：如果您的用例限制您使用非2次幂的表大小，请填充它以使其成为2的幂。如果必须的话，丢弃所有将进入填充的随机数，但是让物理表成为2的幂以避免除法。 - cmaster - reinstate monica

那个填充可以是任何东西，甚至可以是将结果地址与正确的最大值进行比较，并在不在范围内时将其丢弃。这将节省空间开销。然而，如果条件也很关键，则尽量避免使用它。在大多数情况下，只需在填充部分进行一些虚假计算即可摆脱它，并且没有回报值。当然，由于填充而产生额外开销，但该开销少于所有情况下有用工作量，即在大多数情况下是明智的权衡。 - cmaster - reinstate monica

Answer 4

你可以通过循环填充末尾空格的方式，将你的直方图扩展到“大”二次幂。填充时使用一些虚拟值（保证在真实数据中永远不会出现）。例如，给定一个直方图。

[10, 5, 6]

将其扩展到长度为16，如下所示（假设-1是适当的哨兵）：

[10, 5, 6, 10, 5, 6, 10, 5, 6, 10, 5, 6, 10, 5, 6, -1]

那么，可以通过二进制掩码 histogram[prng() & mask] 进行抽样，其中 mask = (1 << new_length) - 1，并检查哨兵值以进行重试，即：

int value;
do {
    value = histogram[prng() & mask];
} while (value == SENTINEL);

// use `value` here

该扩展长度超出必要范围，以确保大多数元素有效（例如，在上面的示例中，只有1/16的查找会“失败”，通过将其扩展到例如64，可以进一步减少此比率）。您甚至可以对检查使用“分支预测”提示（例如 GCC中的__builtin_expect），以便编译器为value != SENTINEL的情况优化代码顺序，这通常是普遍情况。这非常注重内存与速度的平衡。

- huon

---

现代x86实现不支持分支预测提示。哦，算了，你说的是编译器提示。 - Oliver Charlesworth

然而，对我来说并不清楚（仅仅因为我手头没有证据）分支预测失误率是否超过使用模数的成本。 - Oliver Charlesworth

@dbaupp ：这是和 cmaster 提到的一样的想法（填充）。但不幸的是，我不知道如何在 N 不是 2 的幂的情况下“填充”一个 4 兆直方图，而又不失去原来直方图的统计特性。拒绝采样是个好主意，但由于数据体积变大，可能会慢。虽然我会尝试，因为它很有趣。编辑 - 我真的很喜欢这个想法适用于小直方图大小。 - blondiepassesby

@blondiepassesby：好的，明白了。你使用的处理器带有4MB以上的L2缓存，请问是哪种呢？ - Oliver Charlesworth

@blondiepassesby，是的，我想缓存效应可能会很糟糕；在实施这样的方案之前，一定值得进行测量。（另外，你真的需要64位整数吗？使用32位整数可以将数据大小减半，或者如果你非常幸运，16位整数也可以工作。） - huon

Answer 5

以下是一些补充其他优秀答案的想法：

1. 百分之多少的时间花在取模运算上？你怎么知道这个百分比是多少？我之所以问，是因为有时候人们会说某个操作非常缓慢，但实际上，它只占用不到10%的时间，并且他们认为它很大，只是因为他们使用的是一个愚蠢的自己计时的分析器。(我很难想象一个取模运算会比一个随机数生成器占用更多时间。)

2. 何时确定桶的数量？如果它不经常更改，您可以编写一个程序生成器。当桶的数量发生变化时，自动打印出一个新的程序，编译，链接，然后用于您的大规模执行。这样，编译器将知道桶的数量。

3. 您是否考虑过使用准随机数生成器，而不是伪随机数生成器？它可以在更少的样本中提供更高的积分精度。

4. 桶的数量是否可以减少而不会对积分的准确性造成太大影响？

- Mike Dunlavey

---

Monte-Carlo积分的总体时间中有25%用于取模。使用perf和gprof进行验证。另外25%的时间用于访问直方图，当它不适合L2时。我使用的随机数生成器是AVX优化版本的Mersenne Twister，可以生成大批量的随机数。它只占总计算时间的7%。取模是致命的。 - blondiepassesby

这是一个有趣的想法，我将进一步尝试实验（假设它比AVX-MT更快，并且我不必对QRNG的输出取模以查看我的直方图）。 - blondiepassesby

4. 简短回答：我不知道。长篇回答：直方图数据来自现实世界，我没有足够的统计背景知识，不知道在我的计算受到严重影响之前可以承受多少数据丢失。

- blondiepassesby

@blondiepassesby：1. 我有点像个疯子，告诉人们 gprof 的所有错误，比如它不能提供行或指令级别的分辨率。（我会惊讶如果 perf 更好。）在我信任百分比之前，我会使用该帖子中推荐的诊断方法。无论如何，我都不会假设只有一个要修复的问题。如果你能用模数获得 30% 的加速，那很好，但可能还有更多。这就是为什么。即其他 68% 的时间是什么？ - Mike Dunlavey

perf比gprof好得多。实际上，它甚至比shark更好。它确实可以给出指令级分辨率，并且在给定正确标志的情况下可以正确处理流水线（即您不会看到movq支付上面idiv 3 5指令的实际成本：归因更或多或少地完成）。如果您对基准测试充满热情，我有点惊讶您以前没有看过perf。 - blondiepassesby

Answer 6

非均匀性是指dbaupp所警告的，可以通过拒绝和重新绘制值不少于M*(2^64/M)（在取模之前）来避免。
如果M可以用不超过32位表示，则可以通过重复乘法（参见David Eisenstat的答案）或divmod获得多个小于M的值；或者，您可以使用位运算来挑选足够长的位模式以适应M，同样拒绝不少于M的值。
（如果随机数生成没有被模数所淹没，我会感到惊讶的，因为它需要更多的时间/周期/能量消耗。）

- greybeard

---

特别是针对您的最后一句话。 - Mike Dunlavey

@老程序员：模数运算比随机数发生器开销更大。当批处理和向量化时，Mersenne Twister的成本接近于零。它通常可以在不到3个周期内产生32位随机数。 - blondiepassesby

@blondiepassesby：由于无法在上面进行评论，那么反复从直方图中取样有什么用处呢？我认为从中取样填充直方图是有价值的。 - greybeard

@老程序员：直方图是从真实世界（物理）过程的观察值构建的。在某种意义上，直方图是真实世界过程的“模型”。我正在尝试将这个物理过程纳入到一个更大系统的模拟中，并尝试通过蒙特卡罗积分计算最终状态。为了做到这一点，我需要从直方图中抽取许多随机样本。 - blondiepassesby

Answer 7

为了填充桶，您可以使用std::二项分布直接将每个桶填满，而不是一个样本一个样本地填充桶：

以下可能会有所帮助：

int nrolls = 60; // number of experiments
const std::size_t N = 6;
unsigned int bucket[N] = {};

std::mt19937 generator(time(nullptr));

for (int i = 0; i != N; ++i) {
    double proba = 1. / static_cast<double>(N - i);
    std::binomial_distribution<int> distribution (nrolls, proba);
    bucket[i] = distribution(generator);
    nrolls -= bucket[i];
}

实时示例

- Jarod42

---

我不确定我理解你的想法是如何运作的。你能解释一下“喂桶”是什么意思，二项分布是如何发挥作用的吗？bucket[i]代表什么？我必须承认：你的代码让我有点困惑。 - blondiepassesby

我误读了，使用了feed而不是pull，很抱歉。据我所知，您循环选择每次要选择哪个直方图（桶），我建议计算您采取给定直方图的次数。与其掷骰子60次并查看分布，我通过二项式分布计算每个数字被掷出的次数。 - Jarod42

Answer 8

感谢大家的建议。

首先，我现在完全相信取模运算确实是非常恶劣的。
它既非常慢，而且在大多数情况下产生的结果都是不正确的。

在实施和测试了很多建议之后，
@Gene提出的解决方案似乎是速度/质量折衷中最好的：

1. 预先计算normalizer：

   auto normalizer = histogram.size() / (1.0+urng.max());

2. 使用以下代码绘制样本：

   return histogram[ (uint32_t)floor(urng() * normalizer);

这是我迄今尝试过的所有方法中最快的，而且据我所知，
它产生的分布要好得多，即使可能不像拒绝法那样完美。

编辑：我实现了David Eisenstat的方法，这与Jarkkol的建议基本相同：index = (rng() * N) >> 32。它的效果与浮点规范化一样好，而且速度稍快（事实上快9%）。所以现在这是我的首选方式。

- blondiepassesby

---

将浮点数和定点数来回转换真的比直接使用定点数更快吗？让我感到惊讶。 - David Eisenstat

@DavidEisenstat：我还没有尝试过固定点方法（拟随机积分路线似乎更有前途）。你是指你建议的方法，还是Jarkkol的方法？ - blondiepassesby

区别就在于一个（便宜）的移位指令，所以哪个都可以？ - David Eisenstat

Answer 9

你可以使用定点数学，即整数乘法和位移，而不是整数除法。例如，如果你的prng()返回0-65535范围内的值，并且你想将其量化为0-99范围内的值，则可以执行(prng()*100)>>16操作。只需确保乘法不会溢出你的整数类型，因此你可能需要将prng()的结果向右移动。请注意，与取模相比，此映射更好，因为它保留了均匀分布。

- JarkkoL