我想高效地实现将一个无符号a整数按照任意2的幂进行向上舍入的除法。因此,从数学上讲,我需要的是ceiling(p/q)
0。在C语言中,以下函数是一个不考虑q
受限域的草稿实现1:
/** q must be a power of 2, although this version works for any q */
uint64_t divide(uint64_t p, uint64_t q) {
uint64_t res = p / q;
return p % q == 0 ? res : res + 1;
}
当然,我不想在机器级别使用除法或取模运算,因为即使在现代硬件上,这需要很多循环。我正在寻找一种使用移位和/或其他廉价操作的强度缩减 - 利用了q
是2的幂次方的事实。
您可以假设我们拥有高效的lg(unsigned int x)
函数,如果x
是2的幂,则返回x
的以2为底的对数。
如果q
为零,则未定义的行为是可以接受的。
请注意,简单的解决方案:(p+q-1) >> lg(q)
在一般情况下不起作用 - 例如,请尝试p == 2 ^ 64-100 and q == 256
2。
平台细节
我对C语言中表现良好且可能在各种平台上表现良好的解决方案感兴趣,但是为了具体化、奖励以及任何关于性能的明确讨论都需要包含目标架构,我将具体说明如何测试它们:
- Skylake CPU
gcc 5.4.0
,编译标志为-O3 -march=haswell
使用gcc内置函数(例如bitscan/leading zero内置函数)是可以的,特别是我已经按照以下方式实现了可用的lg()
函数:
inline uint64_t lg(uint64_t x) {
return 63U - (uint64_t)__builtin_clzl(x);
}
inline uint32_t lg32(uint32_t x) {
return 31U - (uint32_t)__builtin_clz(x);
}
我验证了这些代码至少在使用-march=haswell
时会被编译成单个 bsr
指令,尽管看起来可能涉及减法。当然,你可以忽略这些并在解决方案中使用任何其他内置函数。
基准测试
我为现有答案编写了一个基准测试,并将根据更改进行共享和更新结果。
为小型、可能被内联的操作编写一个良好的基准测试非常困难。当代码被内联到调用站点时,函数的大部分工作可能会消失,特别是当它在循环中时3。
您可以通过确保代码不被内联来简单地避免整个内联问题:在另一个编译单元中声明它。我试图使用bench
二进制文件来实现这一点,但实际上结果相当无意义。几乎所有实现都以4或5个时钟周期的速度平局结束,但即使只执行一个什么都不做的虚拟方法,例如return 0
,它也需要相同的时间。因此,你大多数时候只是测量call + ret
的开销。此外,你几乎永远不会以这种方式使用函数——除非你出错了,否则它们将可用于内联,这会改变一切。
因此,我将重点关注两个基准测试,它们反复调用要测试的方法,允许内联、跨函数优化、循环提升甚至矢量化。
总体上有两种基准测试类型:延迟和吞吐量。关键区别在于,在延迟基准测试中,每个对divide
的调用都依赖于前一个调用,因此通常不能轻松重叠调用4:
uint32_t bench_divide_latency(uint32_t p, uint32_t q) {
uint32_t total = p;
for (unsigned i=0; i < ITERS; i++) {
total += divide_algo(total, q);
q = rotl1(q);
}
return total;
}
注意,运行的total
取决于每个除法调用的输出,并且它也是divide
调用的输入之一。另一方面,吞吐量变体不会将一个除法的输出馈送到后续的除法中。这允许一个调用的工作与后续调用重叠(由编译器和CPU特别是实现),甚至允许矢量化:
uint32_t bench_divide_throughput(uint32_t p, uint32_t q) {
uint32_t total = p;
for (unsigned i=0; i < ITERS; i++) {
total += fname(i, q);
q = rotl1(q);
}
return total;
}
请注意,在这里我们将循环计数器作为被除数输入 - 这是一个变量,但它不依赖于前一个除法调用。
此外,每个基准测试都有三种除数q
的行为方式:
- 编译时常数除数。例如,对
divide(p, 8)
的调用。在实践中很常见,当除数在编译时已知时,代码可以简单得多。 - 不变的除数。这里的除数在编译时不知道,但对于整个基准测试而言是恒定的。这允许部分与编译时常数相同的优化。
- 可变的除数。除数在每次循环迭代中发生变化。上面的基准函数显示了这种变体,使用“向左旋转1”指令来变换除数。
将所有内容组合起来,您共有6个不同的基准测试。
结果
总的来说
为了选出最佳算法,我查看了每个提议算法的12个子集:(延迟,吞吐量) x (常数a,不变q,可变q) x (32位,64位)
,并根据以下方式每个子测试分配2、1或0分数:
- 最佳算法(在5%的容差范围内)获得2分。
- “足够接近”的算法(比最佳算法慢不超过50%)获得1分。
- 其余算法得分为零。
因此,最大总分是24分,但没有算法达到该得分。这是总体结果:
╔═══════════════════════╦═══════╗
║ Algorithm ║ Score ║
╠═══════════════════════╬═══════╣
║ divide_user23_variant ║ 20 ║
║ divide_chux ║ 20 ║
║ divide_user23 ║ 15 ║
║ divide_peter ║ 14 ║
║ divide_chrisdodd ║ 12 ║
║ stoke32 ║ 11 ║
║ divide_chris ║ 0 ║
║ divide_weather ║ 0 ║
╚═══════════════════════╩═══════╝
因此,针对这个特定的测试代码、这个特定的编译器和这个平台,用户2357112的“variant” ... + (p & mask) != 0
表现最佳,与chux的建议(实际上是相同的代码)并列。
以下是所有子分数,它们合计为以上得分:
╔══════════════════════════╦═══════╦════╦════╦════╦════╦════╦════╗
║ ║ Total ║ LC ║ LI ║ LV ║ TC ║ TI ║ TV ║
╠══════════════════════════╬═══════╬════╬════╬════╬════╬════╬════╣
║ divide_peter ║ 6 ║ 1 ║ 1 ║ 1 ║ 1 ║ 1 ║ 1 ║
║ stoke32 ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║
║ divide_chux ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_user23 ║ 8 ║ 1 ║ 1 ║ 2 ║ 2 ║ 1 ║ 1 ║
║ divide_user23_variant ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_chrisdodd ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║
║ divide_chris ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
║ divide_weather ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
║ ║ ║ ║ ║ ║ ║ ║ ║
║ 64-bit Algorithm ║ ║ ║ ║ ║ ║ ║ ║
║ divide_peter_64 ║ 8 ║ 1 ║ 1 ║ 1 ║ 2 ║ 2 ║ 1 ║
║ div_stoke_64 ║ 5 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 1 ║
║ divide_chux_64 ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_user23_64 ║ 7 ║ 1 ║ 1 ║ 2 ║ 1 ║ 1 ║ 1 ║
║ divide_user23_variant_64 ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║
║ divide_chrisdodd_64 ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║
║ divide_chris_64 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
║ divide_weather_64 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║
╚══════════════════════════╩═══════╩════╩════╩════╩════╩════╩════╝
每个测试都像XY这样命名,其中X属于{延迟, 吞吐量},Y属于{Constant Q, Invariant Q, Variable Q}。例如,LC是“具有恒定q的延迟测试”。
分析
在最高级别上,解决方案可以大致分为两类:快速(前6名)和慢速(后两名)。差异很大:所有的快速算法都至少在两个子测试中最快,通常当它们没有第一时,它们会落入“足够接近”的范畴(只有在stoke
和chrisdodd
的向量化失败的情况下才有例外)。然而,慢速算法在每个测试中都得0分(甚至都不接近)。因此,您可以基本上淘汰慢速算法。
自动向量化
在快速算法中,一个重要的区分因素是能否进行自动向量化。
在延迟测试中,没有算法能够自动向量化,这是有道理的,因为延迟测试旨在直接将其结果馈送到下一次迭代中。因此,您只能按顺序计算结果。
然而,在吞吐量测试中,许多算法能够自动向量化常量Q和不变Q的情况。在这两个测试中,除数q
是循环不变的(并且在前一种情况下它是编译时常量)。被除数是循环计数器,因此它是可变的,但是可预测的(尤其是通过将8添加到前一个输入向量可以轻松地计算出被除数的向量:[0, 1, 2,...,7] + [8, 8,...,8] == [8, 9, 10,...,15]
)。
在这种情况下,gcc
能够向量化peter
,stoke
,chux
,user23
和user23_variant
。它无法向量化chrisdodd
,可能是因为它包含了一个分支(但条件语句并不严格防止向量化,因为许多其他解决方案具有条件元素但仍然向量化)。影响非常大:向量化的算法在吞吐量上显示出约8倍的改进,而不是其他快速算法。
然而,向量化并不是免费的!以下是每个函数“constant”变体的函数大小,Vec?
列显示函数是否向量化:
Size Vec? Name
045 N bench_c_div_stoke_64
049 N bench_c_divide_chrisdodd_64
059 N bench_c_stoke32_64
212 Y bench_c_divide_chux_64
227 Y bench_c_divide_peter_64
220 Y bench_c_divide_user23_64
212 Y bench_c_divide_user23_variant_64
The trend is clear - 向量化函数的大小约为非向量化函数的4倍。这是因为核心循环本身更大(向量指令倾向于更大,而且有更多的指令),以及循环设置和特别是后循环代码要大得多:例如,向量化版本需要缩减以汇总所有向量中的部分和。循环计数是固定的并且是8的倍数,因此不会生成尾代码-但是如果可变,则生成的代码将更大。此外,尽管运行时间有大幅改善,但gcc的向量化实际上很差。以下是Peter例程向量化版本的摘录:
on entry: ymm4 == all zeros
on entry: ymm5 == 0x00000001 0x00000001 0x00000001 ...
4007a4: c5 ed 76 c4 vpcmpeqd ymm0,ymm2,ymm4
4007ad: c5 fd df c5 vpandn ymm0,ymm0,ymm5
4007b1: c5 dd fa c0 vpsubd ymm0,ymm4,ymm0
4007b5: c5 f5 db c0 vpand ymm0,ymm1,ymm0
这段代码独立于8个DWORD
元素,在ymm2
中生成。如果我们将x
作为ymm2
中的单个DWORD
元素,将y
作为ymm1
的传入值,则这四个指令对应于:
x == 0 x != 0
x = x ? 0 : -1; // -1 0
x = x & 1; // 1 0
x = 0 - x; // -1 0
x = y1 & x; // y1 0
因此,前三个指令可以简单地替换为第一个指令,因为两种情况下的状态都是相同的。这样在该依赖链中添加了两个周期(它没有循环执行),并增加了两个额外的uops。显然 gcc
的优化阶段与向量化代码存在某种不良交互,在标量代码中很少会错过如此微不足道的优化。检查其他向量化版本 同样显示性能有很大损失。
分支 vs 无分支
几乎所有的解决方案都编译为无分支代码,即使C代码中有条件语句或显式分支。条件部分足够小,以至于编译器通常决定使用条件移动或某个变体。一个例外是 chrisdodd
,在所有吞吐量测试中都编译了一个分支(检查是否 p == 0
),但在延迟测试中没有编译出来。以下是来自常量q
吞吐量测试的一个典型示例:
0000000000400e60 <bench_c_divide_chrisdodd_32>:
400e60: 89 f8 mov eax,edi
400e62: ba 01 00 00 00 mov edx,0x1
400e67: eb 0a jmp 400e73 <bench_c_divide_chrisdodd_32+0x13>
400e69: 0f 1f 80 00 00 00 00 nop DWORD PTR [rax+0x0]
400e70: 83 c2 01 add edx,0x1
400e73: 83 fa 01 cmp edx,0x1
400e76: 74 f8 je 400e70 <bench_c_divide_chrisdodd_32+0x10>
400e78: 8d 4a fe lea ecx,[rdx-0x2]
400e7b: c1 e9 03 shr ecx,0x3
400e7e: 8d 44 08 01 lea eax,[rax+rcx*1+0x1]
400e82: 81 fa 00 ca 9a 3b cmp edx,0x3b9aca00
400e88: 75 e6 jne 400e70 <bench_c_divide_chrisdodd_32+0x10>
400e8a: c3 ret
400e8b: 0f 1f 44 00 00 nop DWORD PTR [rax+rax*1+0x0]
400e76
分支跳过了p == 0
的情况。实际上,编译器可以将第一次迭代“削减”出来(显式计算其结果),然后完全避免跳转,因为在此之后它可以证明p != 0
。在这些测试中,该分支是完全可预测的,这可能会给实际使用分支编译的代码带来优势(因为比较和分支代码本质上是离线的,且非常接近免费),并且这也是为什么chrisdodd
能够在吞吐量、变量q的情况下获胜的一个重要原因。
详细测试结果
在这里您可以找到一些详细的测试结果以及一些测试本身的细节。
延迟
下面的结果针对每个算法进行了10亿次迭代的测试。周期数是通过将时间/调用乘以时钟频率来计算的。您通常可以假设类似于4.01
与4.00
相同,但像5.11
这样更大的偏差似乎是真实且可重复的。
divide_plusq_32
的结果使用(p + q - 1) >> lg(q)
,但仅供参考,因为此函数对于大的p + q
会失败。 dummy
的结果是一个非常简单的函数:return p + q
,并允许您估算基准测试开销5(加法本身最多需要一个周期)。
==============================
Bench: Compile-time constant Q
==============================
Function ns/call cycles
divide_peter_32 2.19 5.67
divide_peter_64 2.18 5.64
stoke32_32 1.93 5.00
stoke32_64 1.97 5.09
stoke_mul_32 2.75 7.13
stoke_mul_64 2.34 6.06
div_stoke_32 1.94 5.03
div_stoke_64 1.94 5.03
divide_chux_32 1.55 4.01
divide_chux_64 1.55 4.01
divide_user23_32 1.97 5.11
divide_user23_64 1.93 5.00
divide_user23_variant_32 1.55 4.01
divide_user23_variant_64 1.55 4.01
divide_chrisdodd_32 1.95 5.04
divide_chrisdodd_64 1.93 5.00
divide_chris_32 4.63 11.99
divide_chris_64 4.52 11.72
divide_weather_32 2.72 7.04
divide_weather_64 2.78 7.20
divide_plusq_32 1.16 3.00
divide_plusq_64 1.16 3.00
divide_dummy_32 1.16 3.00
divide_dummy_64 1.16 3.00
==============================
Bench: Invariant Q
==============================
Function ns/call cycles
divide_peter_32 2.19 5.67
divide_peter_64 2.18 5.65
stoke32_32 1.93 5.00
stoke32_64 1.93 5.00
stoke_mul_32 2.73 7.08
stoke_mul_64 2.34 6.06
div_stoke_32 1.93 5.00
div_stoke_64 1.93 5.00
divide_chux_32 1.55 4.02
divide_chux_64 1.55 4.02
divide_user23_32 1.95 5.05
divide_user23_64 2.00 5.17
divide_user23_variant_32 1.55 4.02
divide_user23_variant_64 1.55 4.02
divide_chrisdodd_32 1.95 5.04
divide_chrisdodd_64 1.93 4.99
divide_chris_32 4.60 11.91
divide_chris_64 4.58 11.85
divide_weather_32 12.54 32.49
divide_weather_64 17.51 45.35
divide_plusq_32 1.16 3.00
divide_plusq_64 1.16 3.00
divide_dummy_32 0.39 1.00
divide_dummy_64 0.39 1.00
==============================
Bench: Variable Q
==============================
Function ns/call cycles
divide_peter_32 2.31 5.98
divide_peter_64 2.26 5.86
stoke32_32 2.06 5.33
stoke32_64 1.99 5.16
stoke_mul_32 2.73 7.06
stoke_mul_64 2.32 6.00
div_stoke_32 2.00 5.19
div_stoke_64 2.00 5.19
divide_chux_32 2.04 5.28
divide_chux_64 2.05 5.30
divide_user23_32 2.05 5.30
divide_user23_64 2.06 5.33
divide_user23_variant_32 2.04 5.29
divide_user23_variant_64 2.05 5.30
divide_chrisdodd_32 2.04 5.30
divide_chrisdodd_64 2.05 5.31
divide_chris_32 4.65 12.04
divide_chris_64 4.64 12.01
divide_weather_32 12.46 32.28
divide_weather_64 19.46 50.40
divide_plusq_32 1.93 5.00
divide_plusq_64 1.99 5.16
divide_dummy_32 0.40 1.05
divide_dummy_64 0.40 1.04
吞吐量
以下是吞吐量测试的结果。请注意,这里的许多算法都是自动向量化的,因此对于那些算法来说,性能相对非常好:在许多情况下只需要不到一个周期的时间。其中一个结果是,与大多数延迟结果不同,64位函数的速度要慢得多,因为向量化在较小的元素大小上更有效(尽管差距比我预期的要大)。
==============================
Bench: Compile-time constant Q
==============================
Function ns/call cycles
stoke32_32 0.39 1.00
divide_chux_32 0.15 0.39
divide_chux_64 0.53 1.37
divide_user23_32 0.14 0.36
divide_user23_64 0.53 1.37
divide_user23_variant_32 0.15 0.39
divide_user23_variant_64 0.53 1.37
divide_chrisdodd_32 1.16 3.00
divide_chrisdodd_64 1.16 3.00
divide_chris_32 4.34 11.23
divide_chris_64 4.34 11.24
divide_weather_32 1.35 3.50
divide_weather_64 1.35 3.50
divide_plusq_32 0.10 0.26
divide_plusq_64 0.39 1.00
divide_dummy_32 0.08 0.20
divide_dummy_64 0.39 1.00
==============================
Bench: Invariant Q
==============================
Function ns/call cycles
stoke32_32 0.48 1.25
divide_chux_32 0.15 0.39
divide_chux_64 0.48 1.25
divide_user23_32 0.17 0.43
divide_user23_64 0.58 1.50
divide_user23_variant_32 0.15 0.38
divide_user23_variant_64 0.48 1.25
divide_chrisdodd_32 1.16 3.00
divide_chrisdodd_64 1.16 3.00
divide_chris_32 4.35 11.26
divide_chris_64 4.36 11.28
divide_weather_32 5.79 14.99
divide_weather_64 17.00 44.02
divide_plusq_32 0.12 0.31
divide_plusq_64 0.48 1.25
divide_dummy_32 0.09 0.23
divide_dummy_64 0.09 0.23
==============================
Bench: Variable Q
==============================
Function ns/call cycles
stoke32_32 1.16 3.00
divide_chux_32 1.36 3.51
divide_chux_64 1.35 3.50
divide_user23_32 1.54 4.00
divide_user23_64 1.54 4.00
divide_user23_variant_32 1.36 3.51
divide_user23_variant_64 1.55 4.01
divide_chrisdodd_32 1.16 3.00
divide_chrisdodd_64 1.16 3.00
divide_chris_32 4.02 10.41
divide_chris_64 3.84 9.95
divide_weather_32 5.40 13.98
divide_weather_64 19.04 49.30
divide_plusq_32 1.03 2.66
divide_plusq_64 1.03 2.68
divide_dummy_32 0.63 1.63
divide_dummy_64 0.66 1.71
a指定unsigned,我们避免了与C和C++中有符号整数的右移行为相关的所有问题。
0当然,在C中,这种记法实际上不起作用,因为/
截断结果,所以ceiling
什么都没做。因此,请将其视为伪记法,而不是直接使用C。
1我还对所有类型均为uint32_t
而不是uint64_t
的解决方案感兴趣。
2一般来说,任何使p + q >= 2^64
成立的p
和q
都会导致问题,因为会溢出。
3尽管如此,这个函数应该放在一个循环中,因为只有当它在一个相当紧密的循环中调用时,才真正需要半打周期的微小函数的性能。
4这有点简化了 - 只有被除数p
依赖于前一次迭代的输出,因此一些与q
相关的处理工作仍然可以重叠。
5然而,在使用这样的估算时要谨慎 - 开销不仅仅是加法。如果开销显示为4个周期,而某个函数f
需要5个周期,则可能不准确地说f
中实际工作的代价为5-4 == 1
,因为执行方式是重叠的。
double
类型,当然了。你说得对,它不是2的幂次方。我的例子和符号是误导性的。我现在已经更加清晰地表达了。 - BeeOnRopedividend==0
能够正确工作。但我也对特殊情况的解决方案感兴趣,例如,如果它们比一般版本有显著改进,那么这些特殊情况也很重要。但为了赏金的目的,它应该对所有输入都能正确工作。对于赏金,我使用了-03 -march=haswell
,并在Skylake上运行它。对于赏金,请使用C和gcc内置函数,但当然我也对能够击败C等效版本的汇编版本感兴趣。被接受的答案应该是全面优秀的,并且可能与赏金获奖者不同! - BeeOnRope