快速随机整数

2022年5月12日 · 作者:Raimo Niskanen

当您需要“随机”整数,并且快速且廉价地生成它们至关重要时;那么 rand 模块中功能齐全的伪随机数生成器可能就有点过头了。这篇博文将深入探讨该模块的新增功能、即时编译器如何优化它们、已知技巧,并尝试比较这些苹果和土豆。

目录 #

速度重于质量?#

rand 模块中实现的伪随机数生成器提供了许多有用的功能,例如可重复的序列、无偏差的范围生成、任意大小的范围、不重叠的序列、生成浮点数、正态分布浮点数等等。其中许多功能是通过插件框架实现的,但会带来性能成本。

rand 模块提供的不同算法经过选择,以具有出色的统计质量并在严格的 PRNG 测试中表现良好(请参阅 PRNG 测试部分)。

大多数这些算法是为具有 64 位算术(无符号)的机器设计的,但在 Erlang 中,此类整数会变成大数,其处理速度几乎比即时整数慢一个数量级。

64 位 VM 中的 Erlang 项是标记的 64 位字。即时整数的标签为 4 位,为有符号整数值留下 60 位。因此,最大的正即时整数值为 259-1。

许多算法处理无符号整数,因此我们有 59 位可用于此目的。理论上可以使用负值和正值的分裂代码路径来假装 60 位无符号数,但这非常不切实际。

我们决定在这种情况下选择 58 位无符号整数,因为这样我们就可以例如将两个整数相加,并检查溢出或简单地掩码回 58 位,而中间结果不会变成大数。使用 59 位整数需要在甚至进行加法之前检查溢出,因此避免大数的代码会消耗掉避免大数而获得的大部分速度。所以就用 58 位整数!

在 Erlang 中表现良好的算法是那些经过重新设计以处理 58 位整数的算法。但即便如此,在 Erlang 中执行时,它们也远不如其 C 语言起源的速度快。在 Erlang 中获得良好的 PRNG 质量比在 C 中花费更多。在 测量结果部分,我们看到在 C 中以亚纳秒速度运行的 exsp 算法在 Erlang 中需要 17 纳秒。

在这方面,32 位 Erlang 是一个悲伤的故事。这种 Erlang 系统上的大数限制非常低,计算必须使用 26 位整数,因此设计一个不使用大数的 PRNG 必须在周期和大小上都非常小,以至于变得太差而无法使用。已知的技巧 erlang:phash2(erlang:unique_integer(), Range) 仍然相当快,但所有 rand 生成器的工作方式与 64 位系统完全相同,因此它们处理大数,因此速度慢得多。

如果您的应用程序需要一个“随机”整数用于非关键目的,例如选择一个工作进程、选择一条路线等,并且性能比可重复性和统计质量重要得多,那么有哪些选择?

建议的解决方案 #

推理和测量结果在以下部分中,但简而言之

  • 我们认为编写 NIF 并不能获得值得付出的努力的性能。
  • 编写 BIF 也是如此。但是,……BIF(也许还有 NIF)可以实现性能和质量的组合,这是以任何其他方式都无法实现的。如果对这种组合的需求很高,我们可以重新考虑这个决定。
  • 使用系统时间是一个坏主意。
  • erlang:phash2(erlang:unique_integer(), Range) 有其用例。
  • 我们已经实现了一个简单的 PRNG 来填补非关键但速度非常快的数字生成的空白:mwc59

无论如何都使用 rand #

rand 真的慢吗?嗯,也许考虑到它所做的事情,就不是了。

在本文章末尾的 测量结果中,它显示使用 rand 模块的默认算法生成高质量的随机数需要 45 纳秒。

尽可能快地生成一个数字(rand:mwc59/1)可以在不到 4 纳秒的时间内完成,但该算法在统计质量方面存在问题。请参阅 PRNG 测试实现 PRNG 部分。

如果可以使用循环变量存储生成器的状态,则使用高质量的算法(rand:exsp_next/1)需要 16 纳秒。

如果您不能在循环变量中存储生成器状态,则会有更多开销,请参阅 存储状态部分。

现在,如果您还需要一个范围不方便的数字,例如不比生成器大小小很多,您可能必须实现一个拒绝和重新采样的循环,甚至连接数字。

必须实现 rand 模块已经提供的这么多功能的代码开销很容易接近其 26 纳秒的开销,因此通常没有必要重新实现这个轮子……

编写一个 BIF #

Erlang 论坛上有一个讨论主题:寻找更快的 RNG。受此启发,Andrew Bennett(又名 potatosalad)编写了一个 实验性的 BIF

建议的 BIF erlang:random_integer(Range) 不提供可重复性、每个调度程序的生成器状态、调度程序之间保证的序列分离以及高生成器质量。所有这些都归功于使用了 rand 模块中良好的生成器之一,但现在是用其原始编程语言 C 编写的,在 BIF 中。

性能比 mwc59 生成器状态更新慢一点,但质量一流。请参阅 测量结果部分。

出现了关于维护负担、还要实现什么等等的问题。例如,我们可能还需要 erlang:random_integer/0erlang:random_float/0 和一些系统信息来获取生成器位数……

如果 BIF 在那里返回一个 27 位整数,则 BIF 可以在 32 位系统上实现良好的性能,这又成为另一个悬而未决的问题。BIF 生成器在生成的数字方面还是在性能方面应该是平台独立的?

编写一个 NIF #

potatosalad 还编写了一个 NIF,因为我们(Erlang/OTP 团队)认为它可能具有足够好的性能。

然而,测量结果表明,开销明显大于 BIF。尽管 NIF 使用了与 BIF 相同的技巧将状态存储在线程特定数据中,但最终其性能与 erlang:phash2(erlang:unique_integer(), Range) 相同,后者比 BIF 慢约 2 到 3 倍。

我们尝试的一种速度改进是让 NIF 生成一个数字列表,并将该列表用作 Erlang 中的缓存。使用此类缓存的性能与 BIF 一样快,但引入了一些问题,例如您必须决定缓存大小,应用程序必须将缓存保留在堆上,并且当在数字范围内生成时,应用程序必须知道在整个缓存中生成相同范围内的数字。

像 BIF 一样,NIF 也可以在 32 位系统上实现良好的性能,但同样存在开放性问题——平台独立的数字还是性能?

使用系统时间 #

一种建议的技巧是使用 os:system_time(microseconds) 来获取数字。这个技巧有一些特殊性

  • 重复调用时,您可能会多次获得相同的数字。
  • 分辨率取决于系统,因此在某些系统上,您甚至会多次获得相同的数字。
  • 在某些情况下,时间可能会向后跳跃并重复。
  • 从历史上看,它一直是一个瓶颈,尤其是在虚拟化平台上。获取操作系统时间比预期的要困难。

请参阅 测量结果 部分,了解此“解决方案”的性能。

哈希一个“唯一”值 #

最佳组合无疑是 erlang:phash2(erlang:unique_integer(), Range)erlang:phash2(erlang:unique_integer()),后者稍微快一点。

erlang:unique_integer/0 旨在返回一个开销非常小的唯一整数。很难找到一个更好的哈希整数的候选者。

erlang:phash2/1,2 是 Erlang 项的当前通用哈希函数。它具有非常适合 32 位 Erlang 系统的默认返回大小,并且具有 Range 参数。范围上限是用 C 中的简单 rem (%) 完成的,这比在 Erlang 中快得多。这仅适用于远小于 32 位的范围,因为如果范围大于 16 位,则范围上限中的偏差开始变得明显。

唉,此解决方案在 PRNG 测试中表现不佳。

请参阅 测量结果 部分,了解此解决方案的性能。

编写一个简单的 PRNG #

为了快速,PRNG 算法的实现不能执行太多操作。操作必须在即时值(而不是大数)上进行,并且函数的返回值必须是即时值(复合项会给垃圾收集器带来负担)。这严重限制了可以使用的强大算法。

我们编写了一个并将其命名为 mwc59,因为它具有 59 位状态,并且最彻底的加扰函数返回一个 59 位的值。还有一个更快的中间加扰函数,它返回一个 32 位值,这是 MWC 生成器的“数字”大小。也可以直接使用状态的低 16 位而不进行加扰。请参阅 实现 PRNG 部分,了解如何设计此生成器以及原因。

作为真正快速且质量低的算法与功能齐全的算法之间的另一个差距填充,rand 中的一个内部函数已导出:rand:exsp_next/1。此函数实现了 Xoroshiro116+,它作为算法 exsp 存在于 rand 插件框架中。它已被导出,因此对于不需要任何框架功能的应用程序,可以在没有插件框架开销的情况下获得良好的质量。

请参阅 测量结果 部分,了解速度比较。

质量 #

PRNG 的质量有很多不同的方面。以下是一些。

周期 #

erlang:phash2(erlang:unique_integer(), Range) 在概念上具有无限周期,因为假设它重复的时间将比 Erlang 节点存活的时间长。

对于新的快速 mwc59 生成器,其周期约为 259。对于 rand 中的常规生成器,它至少为 2116 - 1,这是一个巨大的差异。在 Erlang 节点的生命周期内可能可以消耗 259 个数字,但不能消耗 2116 个。

rand 中也有一些生成器的周期为 2928 - 1,这似乎荒谬地长,但这有助于生成许多保证不重叠的并行子序列。

例如,在物理模拟中,通常只使用生成器周期的一部分,包括生成数字的数量和生成数字的范围,或者某些特定数字不再重复可能会影响模拟。如果你从一副牌中抽出了 3 张 A,你就知道只剩下一张了。

有些应用程序可能对生成器的周期敏感,而另一些则不然,这需要考虑。

大小 #

新的快速 mwc59 生成器的值大小为 59、32 或 16 位,具体取决于所使用的加扰函数。rand 模块中的大多数常规生成器的值大小为 58 位。

如果你需要 2 的幂次范围内的数字,则可以简单地屏蔽掉低位。

V = X band ((1 bsl RangeBits) - 1).

或者向下移动所需的位数。

V = X bsr (GeneratorBits - RangeBits).

这取决于是否已知生成器具有弱高位或弱低位。

如果需要的范围不是 2 的幂次,但仍然远小于生成器的大小,则可以使用 rem

V = X rem Range.

经验法则是 Range 应该小于生成器大小的平方根。这比按位运算慢得多,并且该操作会传播低位,如果已知生成器具有弱低位,则可能会出现问题。

另一种方法是使用截断乘法。

V = (X * Range) bsr GeneratorBits

此处的经验法则是 Range 应小于 2GeneratorBits 的平方根,即 2GeneratorBits/2。此外,X * Range 不应创建大数,因此不应超过 59 位。此方法会传播高位,如果已知生成器具有弱高位,则可能会出现问题。

其他技巧也是可能的,例如,如果你需要 0 到 999 范围内的数字,则可以使用按位运算来获取 0 到 1023 之间的数字,如果太高,则重新尝试,这实际上可能比使用 rem 平均更快。此方法还完全消除了生成数字中的偏差。前面的方法有经验法则来获得如此小的偏差,以至于变得难以察觉。

频谱得分 #

生成器的频谱得分衡量来自生成器的数字序列之间的不相关程度。 N 个数字的序列被解释为 N 维向量,并且维度 N 的频谱得分是对这些向量在 N 维(超)立方体中分布的均匀程度的度量。

os:system_time(microseconds) 只是简单地递增,因此它应该具有很差的频谱得分。

erlang:phash2(erlang:unique_integer(), Range) 的频谱得分未知,因为这不属于哈希函数背后的数学原理。但是,哈希函数旨在为任何输入良好地分配哈希值,因此人们可以希望数字的统计分布无论如何都是体面的且“随机”的。不幸的是,这似乎在PRNG 测试中并不成立。

rand 模块中的所有常规 PRNG 都具有良好的频谱得分。新的 mwc59 生成器在大多数情况下都是这样,但由于其不平衡的设计和 2 的幂次乘数,在 2 维和 3 维中不是。加扰器用于弥补这些缺陷。

PRNG 测试 #

有一些测试框架可以测试 PRNG 的统计属性,例如 TestU01 框架或 PractRand

rand 模块中的常规生成器在这些测试中表现良好,并通过了全面的测试套件。

尽管 mcg59 生成器通过了 PractRand 2 TB 测试和 TestU01 及其低 16 位,没有任何加扰,但当稍微调整测试参数时,其统计问题就会显现出来。为了在更多情况下和更多位中表现良好,需要加扰函数。尽管如此,较小的状态空间和基本生成器的缺陷使其难以完全通过所有测试。但是,通过彻底的双 Xorshift 加扰器,它变得非常好。

erlang:phash2(N, Range) 在递增序列上的 TestU01 中表现不佳,这表明哈希函数具有与 PRNG 不同的设计标准。

但是,这些类型的测试可能与你的应用程序完全无关。

可预测性 #

对于某些应用程序,生成的数字可能必须是加密不可预测的,而对于另一些应用程序,则没有严格的要求。

对于“非关键”应用程序存在一个灰色地带,例如,一个流氓方可能会影响输入数据,并且如果它知道 PRNG 序列,则可以将所有数据引导到哈希表槽、使一个特定的工作进程过载或类似的操作,从而攻击应用程序。而且,一开始是“非关键”的应用程序可能有一天会悄悄地变得对业务至关重要……

这是一个需要考虑的方面。

存储状态 #

如果 PRNG 的状态可以保存在循环变量中,则成本几乎为零。但是,一旦必须将其存储在堆变量中,由于堆数据分配、术语构建和垃圾回收,它将花费性能。

测量结果 部分中,我们看到最快的 PRNG 可以在不到 4 纳秒的时间内生成一个新状态,该状态也是生成的整数。不幸的是,仅在 2 元组中返回该值和新状态大约增加了 10 纳秒。

必须存储 PRNG 状态的应用程序状态通常更加复杂,因此更新它的成本可能会更大。

播种 #

播种与可预测性相关。如果你可以猜测种子,你就知道生成器的输出。

种子与生成器相关,如何创建好的种子通常比生成一个数字需要更长的时间。有时,种子及其可预测性并不那么重要,可以使用常量。如果生成器实例每次播种只生成少量数字,则播种可能会成为更困难的问题。

erlang:phash2(erlang:unique_integer(), Range) 是预先播种的,或者更确切地说,不能播种的,因此它没有播种成本,但另一方面,如果可以估计自节点启动以来已生成多少个唯一整数,则它可能具有相当的可预测性。

rand 模块中的默认播种使用节点名称的哈希值、系统时间和 erlang:unique_integer() 的组合来创建种子,希望该种子足够不可预测。

建议的 NIF 和 BIF 解决方案还需要一种创建足够好的种子的方法,其中“足够好”很难用数字来衡量。

JIT 优化 #

新实现的 mwc59 生成器的速度部分归功于编译器和即时编译 BEAM 代码加载器中最近的基于类型的优化

没有基于类型的优化 #

这是 mwc59 生成器的 Erlang 代码

mwc59(CX) ->
    C = CX band ((1 bsl 32)-1),
    X = CX bsr 32,
    16#7fa6502 * X + C.

使用 no_type_opt 标志禁用基于类型的优化,代码编译为以下 Erlang BEAM 汇编程序 (erlc -S rand.erl)

    {gc_bif,'bsr',{f,0},1,[{x,0},{integer,32}],{x,1}}.
    {gc_bif,'band',{f,0},2,[{x,0},{integer,4294967295}],{x,0}}.
    {gc_bif,'*',{f,0},2,[{x,0},{integer,133850370}],{x,0}}.
    {gc_bif,'+',{f,0},2,[{x,0},{x,1}],{x,0}}.

当由 JIT (x86) 加载 (erl +JDdump true) 时,机器代码变为

# i_bsr_ssjd
    mov rsi, qword ptr [rbx]
# is the operand small?
    mov edi, esi
    and edi, 15
    cmp edi, 15
    short jne L2271

上面是一个测试,判断 {x,0} 是否是一个小整数,如果不是,则调用 L2271 的回退来处理任何术语。

接下来是用于右移的机器代码,Erlang bsr 32,x86 sar rax, 32,以及跳过回退代码

    mov rax, rsi
    sar rax, 32
    or rax, 15
    short jmp L2272
L2271:
    mov eax, 527
    call 140439031217336
L2272:
    mov qword ptr [rbx+8], rax
# line_I

这里是具有类似测试和回退代码的 band

# i_band_ssjd
    mov rsi, qword ptr [rbx]
    mov rax, 68719476735
# is the operand small?
    mov edi, esi
    and edi, 15
    cmp edi, 15
    short jne L2273
    and rax, rsi
    short jmp L2274
L2273:
    call 140439031216768
L2274:
    mov qword ptr [rbx], rax

下面是带有测试、回退代码和溢出检查的 *

# line_I
# i_times_jssd
    mov rsi, qword ptr [rbx]
    mov edx, 2141605935
# is the operand small?
    mov edi, esi
    and edi, 15
    cmp edi, 15
    short jne L2276
# mul with overflow check, imm RHS
    mov rax, rsi
    mov rcx, 133850370
    and rax, -16
    imul rax, rcx
    short jo L2276
    or rax, 15
    short jmp L2275
L2276:
    call 140439031220000
L2275:
    mov qword ptr [rbx], rax

以下是带有测试、回退代码和溢出检查的 +

# i_plus_ssjd
    mov rsi, qword ptr [rbx]
    mov rdx, qword ptr [rbx+8]
# are both operands small?
    mov eax, esi
    and eax, edx
    and al, 15
    cmp al, 15
    short jne L2278
# add with overflow check
    mov rax, rsi
    mov rcx, rdx
    and rcx, -16
    add rax, rcx
    short jno L2277
L2278:
    call 140439031219296
L2277:
    mov qword ptr [rbx], rax

基于类型的优化 #

当编译器可以推断出有关参数的类型信息时,它可以发出更有效的代码。人们希望添加一个限制参数为 59 位整数的保护,但不幸的是,编译器还不能使用这样的保护测试。

但是向 Erlang 代码添加冗余的输入位掩码会将编译器置于正确的轨道上。这是一个蹩脚的解决方案,并且只会在编译器得到改进以从保护中推断出相同信息之前使用。

现在,Erlang 代码具有第一个冗余的 59 位掩码

mwc59(CX0) ->
    CX = CX0 band ((1 bsl 59)-1),
    C = CX band ((1 bsl 32)-1),
    X = CX bsr 32,
    16#7fa6502 * X + C.

然后,BEAM 汇编程序变为,在编译器中使用 OTP-25.0 版本中的默认基于类型的优化

    {gc_bif,'band',{f,0},1,[{x,0},{integer,576460752303423487}],{x,0}}.
    {gc_bif,'bsr',{f,0},1,[{tr,{x,0},{t_integer,{0,576460752303423487}}},
             {integer,32}],{x,1}}.
    {gc_bif,'band',{f,0},2,[{tr,{x,0},{t_integer,{0,576460752303423487}}},
             {integer,4294967295}],{x,0}}.
    {gc_bif,'*',{f,0},2,[{tr,{x,0},{t_integer,{0,4294967295}}},
             {integer,133850370}],{x,0}}.
    {gc_bif,'+',{f,0},2,[{tr,{x,0},{t_integer,{0,572367635452168875}}},
             {tr,{x,1},{t_integer,{0,134217727}}}],{x,0}}.

请注意,在初始输入 band 操作之后,类型信息 {tr,{x_},{t_integer,Range}} 已一直向下传播。

现在,JIT 编译的代码明显缩短了。

输入掩码操作对该值一无所知,因此它具有操作数测试和回退到任何术语代码

# i_band_ssjd
    mov rsi, qword ptr [rbx]
    mov rax, 9223372036854775807
# is the operand small?
    mov edi, esi
    and edi, 15
    cmp edi, 15
    short jne L1816
    and rax, rsi
    short jmp L1817
L1816:
    call 139812177115776
L1817:
    mov qword ptr [rbx], rax

对于所有后续操作,操作数测试和回退代码已优化为成为直接的机器代码序列

# line_I
# i_bsr_ssjd
    mov rsi, qword ptr [rbx]
# skipped test for small left operand because it is always small
    mov rax, rsi
    sar rax, 32
    or rax, 15
L1818:
L1819:
    mov qword ptr [rbx+8], rax
# line_I
# i_band_ssjd
    mov rsi, qword ptr [rbx]
    mov rax, 68719476735
# skipped test for small operands since they are always small
    and rax, rsi
    mov qword ptr [rbx], rax
# line_I
# i_times_jssd
# multiplication without overflow check
    mov rax, qword ptr [rbx]
    mov esi, 2141605935
    and rax, -16
    sar rsi, 4
    imul rax, rsi
    or rax, 15
    mov qword ptr [rbx], rax
# i_plus_ssjd
# add without overflow check
    mov rax, qword ptr [rbx]
    mov rsi, qword ptr [rbx+8]
    and rax, -16
    add rax, rsi
    mov qword ptr [rbx], rax

执行时间从 3.7 纳秒降至 3.3 纳秒,仅通过避免冗余检查和测试就快了 10%,尽管增加了不需要的初始输入掩码操作。

并且还有改进的空间。这些值在操作之间来回移动到 BEAM {x,_} 寄存器 (qword ptr [rbx])。JIT 可以避免从 {x,_} 寄存器移回,因为它有可能知道该值在进程寄存器中。如果编译器发出信息,说明该值在操作之后不会从 {x,_} 寄存器中使用,则可以优化移出到 {x,_} 寄存器。

实现 PRNG #

为了在 Erlang 中创建一个真正快速的 PRNG,语言实现存在一些限制

  • 如果生成器状态是一个复杂项,即堆项,而不是立即值,则状态更新会慢得多。因此,状态应该是一个最大 59 位的整数。
  • 如果中间结果创建了一个大数(bignum),即溢出 59 位,则算术运算会慢得多,因此中间结果必须生成适合 59 位的值。
  • 如果生成器在一个复合项中同时返回生成的值和新状态,那么,再次,更新堆数据会使其慢得多。因此,生成器应该只返回一个立即整数状态。
  • 如果返回的状态整数不能用作生成的数字,则可以使用对状态进行操作的单独值函数。但是,两次调用会使调用开销加倍。

LCG 和 MCG #

第一次尝试是使用经典的 2 的幂线性同余生成器。

X1 = (A * X0 + C) band (P-1)

以及乘法同余生成器。

X1 = (A * X0) rem P

为了避免大数运算,乘积 A * X0 必须适合 59 位。Pierre L’Ecuyer 的经典论文“不同大小和良好晶格结构的线性同余生成器表”列出了两个 35 位生成器,即一个 P = 235 的 LCG 和一个 P 为略低于 235 的质数的 MCG。这些是找到的最大的生成器,其乘法运算不会溢出 59 位。

LCG 的速度非常好。MCG 的速度稍慢,因为它必须执行整数除法 rem,但是由于 P 接近 235,因此可以进行优化,使得速度仅比 LCG 慢约 50%。

不幸的是,2 的幂 LCG 的短周期和已知的怪癖在 PRNG 测试中显示出来。

它们彻底失败了。

MWC #

米兰大学的 Sebastiano Vigna 也帮助设计了我们当前的 58 位 Xorshift 系列生成器,他建议使用带进位的乘法生成器 (Multiply With Carry generator) 代替。

T  = A * X0 + C0,
X1 = T band ((1 bsl Bits)-1),
C0 = T bsr Bits.

此生成器对大小为 Bits 的“数字”进行操作,如果一个数字是机器字的一半,则乘法不会溢出。无需将状态作为数字 X 和进位 C,可以将它们合并为以 T 作为状态。我们得到

X  = T0 band ((1 bsl Bits)-1),
C  = T0 bsr Bits,
T1 = A * X + C

MWC 生成器实际上是具有 2 的幂乘数的 MCG 生成器的另一种形式,因此这是一个等效的生成器。

T0 = (T1 bsl Bits) rem ((A bsl Bits) - 1)

在此形式中,生成器以相反的顺序更新状态,因此交换了 T0T1。模数 (A bsl Bits) - 1 必须是一个安全素数,否则生成器没有最大周期。

基本生成器 #

由于乘数(或其乘法逆元)是 2 的幂,因此 MWC 生成器在 3 维中获得较差的频谱得分,因此必须对状态使用加扰函数才能提高质量。

开始搜索合适的数字大小和乘数,主要通过使用程序来尝试安全素数的乘数并估计频谱得分,例如 CPRNG

当生成器平衡时,即乘数 A 接近 Bits 位时,频谱得分是最好的,除了 3 维中已知的问题。但是,由于无论如何都需要加扰函数,因此有机会尝试使用 27 位乘数生成一个舒适的 32 位数字。使用这些大小,乘积 A * X0 不会创建大数,并且使用 32 位数字,可以使用标准的PRNG 测试在开发期间测试生成器。

不幸的是,由于使用了这种稍微不平衡的参数,2 维的频谱得分也很差,但是加扰器也可以解决这个问题...

最终的生成器是

mwc59(T) ->
    C = T bsr 32,
    X = T band ((1 bsl 32)-1),
    16#7fa6502 * X + C.

此基本生成器的 32 位数字在 PRNG 测试中表现不佳,但实际上低 16 位在 PractRand 中通过了 2 TB,在位反转的情况下通过了 1 TB,这非常好。2 维和 3 维的频谱得分差的问题在于 MWC 数字的较高位。

加扰 #

加扰器必须快速,在使用时仅使用少量快速操作。对于像这样的算术生成器,Xorshift 是合适的加扰器。我们研究了单 Xorshift、双 Xorshift 和双 XorRot。双 XorRot 比双 Xorshift 慢,但效果不佳,可能是因为生成器具有良好的低位,因此需要将它们向上移动以改善高位。将高位向下旋转到低位没有任何改善。

这是一个单 Xorshift 加扰器

V = T bxor (T bsl Shift)

在尝试 Shift 常数时,结果表明,使用较大的移位常数,生成器在 PractRand 中表现更好,而使用较小的移位常数,生成器在生日间隔测试(例如 TestU01 BigCrush 中)和冲突测试中表现更好。唉,不可能找到一个对两者都有效的常数。

选择的单 Xorshift 常数为 8,它在 PractRand 中通过了 4 TB,在 TestU01 中通过了 BigCrush,但在更彻底的生日间隔测试中失败。失败很少,例如 8 维和 9 维中的最低位,以及 2 维和 3 维中的一些中间位。这不太可能影响大多数应用程序,并且如果使用生成的 32 位的高位,这些缺陷应该会隐藏起来。

最终的加扰器必须避免大数运算,并将值屏蔽为 32 位,因此看起来像这样

mwc59_value32(T) ->
    V0 = T  band ((1 bsl 32)-1),
    V1 = V0 band ((1 bsl (32-8))-1),
    V0 bxor (V1 bsl 8).

更好的加扰器是双 Xorshift,它可以同时具有小的移位和大的移位。使用小移位 4 可以使组合生成器在生日间隔和冲突测试中表现出色,然后使用大移位 27 将整个改进的 32 位 MWC 数字一直移动到生成器 59 位状态的最高位。这就是想法,结果证明效果很好。

双 Xorshift 加扰器生成一个 59 位数字,其中低位、高位、反向低位、反向高位等……在 PractRandTestU01 BigCrush 以及详尽的生日间隔和冲突测试中都表现良好。它也没有比单 Xorshift 加扰器慢很多。

这是一个双 Xorshift 加扰器,先移位 4 位,再移位 27 位

V1 = T bxor (T bsl 4),
V  = V1 bxor (V1 bsl 27).

避免了大数运算并产生 59 位值,成为最终的加扰器

mwc59_value(T) ->
    V0 = T  band ((1 bsl (59-4))),
    V1 = T  bxor (V0 bsl 4),
    V2 = V1 band ((1 bsl (59-27))),
    V1 bxor (V2 bsl 27).

非常感谢 Sebastiano Vigna,他完成了大部分(实际上是全部)参数搜索以及生成器和加扰器的广泛测试,并以对可行性的了解为后盾。以这种特殊方式使用 MWC 生成器在数学方面是相当未知的领域,因此广泛的测试是信任生成器质量的方式。

rand_SUITE:measure/1 #

rand 模块的测试套件 — rand_SUITE,在 Erlang/OTP 源代码树中,包含一个测试用例 measure/1。此测试用例是对 rand 模块中所有算法的微基准测试,还包括其他一些算法。它测量每个生成的数字的执行时间(以纳秒为单位),并以绝对值和相对于默认算法 exsss(被认为是 100%)的相对值呈现时间。请参见测量结果

measure/1 即使没有测试框架也可以运行。只要 rand_SUITE.beam 在代码路径中,rand_SUITE:measure(N) 将以 N 作为努力因子运行基准测试。N = 1 是默认值,例如 N = 5 会进行更慢更彻底的测量。

测试用例分为几部分,其中每个部分首先使用默认生成器运行预热,然后运行一个空的基准测试生成器以测量基准测试开销,之后运行特定部分的所有生成器。在开销运行之后,从呈现的结果中减去基准测试开销。

预热和开销测量与补偿是 measure/1 测试用例的最新改进。通过在每个测试用例循环迭代中内联 10 个 PRNG 迭代,也降低了开销,这使得开销降低到没有此类内联的三分之一,并且开销现在与最快的生成器本身一样大,接近 Erlang 中的函数调用开销。

不同的 measure/1 部分是不同的用例,例如“均匀整数半范围 + 1”等。其中许多测试了插件框架功能的性能。本文感兴趣的测试部分是“均匀整数范围 10000”、“均匀整数 32 位”和“均匀整数全范围”。

测量结果 #

以下是从作者的笔记本电脑运行 rand_SUITE:measure(20) 中选择的一些结果

{mwc59,Tag} 生成器是 rand:mwc59/1,其中 Tag 指示是否使用了 raw 生成器、rand:mwc59_value32/1rand:mwc59_value/1 加扰器。

{exsp,_} 生成器是 rand:exsp_next/1,这是一个新导出的内部函数,不使用插件框架。从插件框架调用时,它在下面被调用为 exsp

unique_phash2erlang:phash2(erlang:unique_integer(), Range)

system_timeos:system_time(microsecond)

RNG uniform integer range 10000 performance
                   exsss:     57.5 ns (warm-up)
                overhead:      3.9 ns      6.8%
                   exsss:     53.7 ns    100.0%
                    exsp:     49.2 ns     91.7%
         {mwc59,raw_mod}:      9.8 ns     18.2%
       {mwc59,value_mod}:     18.8 ns     35.0%
              {exsp,mod}:     22.5 ns     41.9%
          {mwc59,raw_tm}:      3.5 ns      6.5%
      {mwc59,value32_tm}:      8.0 ns     15.0%
        {mwc59,value_tm}:     11.7 ns     21.8%
               {exsp,tm}:     18.1 ns     33.7%
           unique_phash2:     23.6 ns     44.0%
             system_time:     30.7 ns     57.2%

前两个是预热和开销测量。测量的开销会从 “overhead:” 行之后的所有测量中减去。此处测得的开销为 3.9 ns,这与 exsss 在预热运行期间比在 overhead 之后多测量 3.8 ns 的情况非常吻合。但是,预热运行有点不可预测。

{_,*mod}system_time 都使用 (X rem 10000) + 1 来实现所需的范围。rem 操作非常昂贵,我们将在下一节中看到。

{_,*tm} 使用截断的乘法运算来实现范围,即 ((X * 10000) bsr GeneratorBits) + 1,这比使用 rem 快得多。

erlang:phash2/2 有一个范围参数,该参数在 BIF 中执行 rem 10000 操作,这相当便宜,我们将在下一节中看到。

RNG uniform integer 32 bit performance
                   exsss:     55.3 ns    100.0%
                    exsp:     51.4 ns     93.0%
        {mwc59,raw_mask}:      2.7 ns      4.9%
         {mwc59,value32}:      6.6 ns     12.0%
     {mwc59,value_shift}:      8.6 ns     15.5%
            {exsp,shift}:     16.6 ns     30.0%
           unique_phash2:     22.1 ns     40.0%
             system_time:     23.5 ns     42.6%

在本节中,为了生成 32 位范围内的数字,{mwc59,raw_mask}system_time 使用位掩码 X band 16#ffffffff{_,*shift} 使用 bsr 来移出低位,而 {mwc59_value32} 本身具有正确的范围。在这里,我们看到位运算比上一节中的 rem 操作快 10 ns。{mwc59,raw_*} 快 3 倍以上。

与上一节中的截断乘法变体相比,这里的位运算快 3 ns。

unique_phash2 仍然使用 BIF 编码的整数除法来实现范围,这使其速度与上一节大致相同,但似乎 2 的幂的整数除法速度更快。

RNG uniform integer full range performance
                   exsss:     45.1 ns    100.0%
                    exsp:     39.8 ns     88.3%
                   dummy:     25.5 ns     56.6%
             {mwc59,raw}:      3.7 ns      8.3%
         {mwc59,value32}:      6.9 ns     15.2%
           {mwc59,value}:      8.5 ns     18.8%
             {exsp,next}:     16.8 ns     37.2%
       {splitmix64,next}:    331.1 ns    734.3%
           unique_phash2:     21.1 ns     46.8%
                procdict:     75.2 ns    166.7%
        {mwc59,procdict}:     16.6 ns     36.8%

在本节中,不进行范围上限。使用原始生成器输出。

在这里,我们有 dummy 生成器,它是 rand 插件框架中一个未记录的生成器,仅执行最小的状态更新并返回一个常数。它在这里用于测量插件框架开销。

插件框架的开销经测量为 25.5 纳秒,与 exsp - {exsp,next} = 23.0 纳秒相当吻合,这在插件框架内外使用了相同的算法,提供了对框架开销的另一个衡量标准。

procdict 是默认算法 exsss,但它使插件框架将生成器状态存储在进程字典中,这里花费了 30 纳秒。

{mwc59,procdict} 将生成器状态存储在进程字典中,这里花费了 12.9 纳秒。存储的状态项比插件框架的小得多。与上一段中的 procdict 进行比较。

总结 #

新的快速生成器在 rand 模块中的函数填补了在质量之上追求速度的空白,其中基于类型的 JIT 优化提升了性能。

高速和高质量的结合只能通过 BIF 实现来实现,但我们希望这是我们不需要处理的组合……

实现一个 PRNG 是一项棘手的工作。

rand_SUITE:measure/1 最近的改进突出了宝贵的 CPU 周期都用在了哪里。