可向量化的 CLZ 与 CTZ

# 向量化 CLZ 与 CTZ 来源：https://purplesyringa.moe/blog/vectorized-clz-and-ctz/ 2026年10月9日 `clz` 和 `ctz` 是用于计算定长整数中前导（尾随）零位数量的指令。现代 CPU 原生支持这些指令，但它们并不总是很快，例如在 Arrow Lake 上，`tzcnt` 的延迟为 `3`。 我在一个正在开发的 FPU 模拟器中使用了 `ctz`，但后来发现可以借助浮点技巧来避免它，并且我刚刚意识到这可以以一种迂回的方式推广为一个可向量化的 `ctz` polyfill。为了完整性，我还实现了 `clz`。 ### clz 我们先从较简单的 `clz` 开始： ``` fn clz(x: u32) -> u32 { let a = 2.0f64.powi(-970); 32 - ((f64::from_bits(a.to_bits() | x as u64) - a).to_bits() >> 52) as u32 } ``` 基本思路如下： 浮点数的指数是其数值的带偏移对数。将一个 32 位数字 x 替换到某个值 2^k 的尾数中，我们得到一个表示 2^k(1 + 2^-52·x) 的双精度浮点数。然后以双精度浮点数减去 2^k，得到 2^(k-52)·x。提取指数可得 k-52+31-clz(x)，从中通过一次按位减法即可算出 `clz`。k 的取值需确保 `clz` 在 x=0 时也能正确工作。 我们需要 64 位双精度浮点数来处理 32 位输入；不幸的是，这意味着该技巧不能用于任意 64 位输入，最多只能处理 52 位。 假设输入和输出存储在 `u64x4` 中，编译结果如下： ``` vpbroadcastq ymm1, [rip + bias] vorpd ymm0, ymm0, [rip + a] vsubpd ymm0, ymm0, [rip + a] vpsrlq ymm0, ymm0, 52 vpsubq ymm0, ymm1, ymm0 a: .quad 0x350000000000000, 0x350000000000000, 0x350000000000000, 0x350000000000000 bias: .quad 32 ``` 在我的 Haswell 上，这以每次迭代 0.45ns 运行，而标量版本为 1ns。当受延迟限制时，数字变为 2ns 对比 1ns（但如果你在向量化 `ctz` 上受延迟限制，那你大概哪里做错了）。 Ian Qvist 在 Alder Lake 上测试了这一点（感谢！），结果为每次迭代 0.29ns，而标量版本为 0.85ns；受延迟限制时为 1.3ns 对比 0.85ns。在现代 Intel CPU 上，结果应该相同或更好。 AMD CPU 使得 `lzcnt` 非常廉价，标量版本可能会胜出。但请记住 Zen 系列 CPU 支持 AVX-512，其中包含 `vplzcntd`，这也是一个选项。 ### ctz 接下来是 `ctz`： ``` fn ctz(x: u32) -> u32 { let a = f64::from_bits((0x340000100000001 ^ x as u64) ^ (x as u64 + u32::MAX as u64)) - f64::from_bits(0x340000000000000); (a.to_bits() >> 52) as u32 } ``` 我们从 x ⊕ (x-1) 开始，以隔离最低的置位比特。`ctz` 等于该值的对数，我们通过按位加法加上 2^k 再以双精度浮点数减去 2^k 来确定它，然后检查指数——在精心选择的 k 下，指数中包含未偏移的 `ctz`。我们预先将 2^32 混入尾数，并用 x+2^32-1 代替 x-1 以正确处理 x=0 的情况，还预先将 1 混入尾数以确保奇数 x 产生 a=0 而非一个缓慢的非规格化数。（你能*想象*我花了多少时间来安排这一切吗？） 该函数编译为： ``` vpxor ymm1, ymm0, [rip + c1] vpaddq ymm0, ymm0, [rip + u32_max] vpxor ymm0, ymm1, ymm0 vsubpd ymm0, ymm0, [rip + c2] vpsrlq ymm0, ymm0, 52 c1: .qword 0x340000100000001, 0x340000100000001, 0x340000100000001, 0x340000100000001 u32_max: .qword 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff c2: .qword 0x340000000000000, 0x340000000000000, 0x340000000000000, 0x340000000000000 ``` 在 Haswell 上，这以每次迭代 0.49ns 运行，受延迟限制时为 2.3ns。在 Alder Lake 上，为每次迭代 0.35ns，受延迟限制时为 1.3ns。与 `clz` 相比的减速是由于多用了一条指令。如果存在 AVX-512，可以通过使用 `vpternlogq` 来避免，但到了那个地步，你还不如直接对 `(x-1) & !x` 运行 `vpopcntd`。标量版本的表现与 `clz` 无异。 > **后补：** Nikolay Malkovsky (https://t.me/a_zachem_eto_nuzhno) 指出 de Bruijn 序列 (https://en.wikipedia.org/wiki/De_Bruijn_sequence) 提供了另一种可向量化的方法。经过一些测试，我得出了以下代码： ``` const char table[32] = { 0, 4, 5, 6, 11, 9, 7, 12, 15, 3, 10, 8, 14, 2, 13, 1, 0, 4, 5, 6, 11, 9, 7, 12, 15, 3, 10, 8, 14, 2, 13, 1, }; __m256i bit = _mm256_andnot_si256(x, _mm256_sub_epi32(x, _mm256_set1_epi32(1))); __m256i high = _mm256_madd_epi16( _mm256_cmpeq_epi16(bit, _mm256_set1_epi16(-1)), _mm256_set1_epi16(-16) ); __m256i index = _mm256_srli_epi32(_mm256_mullo_epi32(bit, _mm256_set1_epi32(0xf0a6f0a7)), 28); __m256i low = _mm256_shuffle_epi8(_mm256_loadu_si256((__m256i*)table), index); return _mm256_add_epi32(low, high); ``` 我们无法使用真正的 32 字节查找表，因为 `vpshufb` 不能跨 16 字节 lane。我所采用的方法很难解释，但本质上我们使用一个重复两次的 16 位 de Bruijn 序列来计算 `ctz` 的第 0-3 位，然后根据哪一半全为零来加上 16 或 32。`0xf0a6f0a7` 是仅有的四个能使此方法有效的魔法常数之一。 在 Haswell 上这需要 1ns（在 Alder Lake 上 0.7ns），但吞吐量是其两倍，因此如果它有助于避免数据洗牌，可能比基于浮点的方法稍微快一些。