max(ctz(x), ctz(y)) 是否有更快的算法?

33
对于min(ctz(x), ctz(y)),我们可以使用ctz(x | y)来获得更好的性能。但是对于max(ctz(x), ctz(y))呢? ctz表示“计算末尾零位数”。
C++版本(Compiler Explorer
#include <algorithm>
#include <bit>
#include <cstdint>

int32_t test2(uint64_t x, uint64_t y) {
    return std::max(std::countr_zero(x), std::countr_zero(y));
}

Rust 版本(编译器浏览器

pub fn test2(x: u64, y: u64) -> u32 {
    x.trailing_zeros().max(y.trailing_zeros())
}

3
单元测试:https://godbolt.org/z/1hY4ch9sh - Marek R
2
请注意,指定处理器架构会使代码变得更加美观。在这种情况下,clang可以将其优化为无分支代码:https://godbolt.org/z/dWse6hxbY - Marek R
3
在ARM上,ctz(x)被实现为clz(rbit(x))。由于我们有max(clz(x), clz(y)) = clz(min(x,y)),这使得我们可以执行clz(min(rbit(x), rbit(y))),从而节省了一个clz。(并且在这种架构上,min很容易无分支地实现。)因此,了解您的架构如何实际执行ctz可能会有所帮助。 - Nate Eldredge
2
@PeterCordes 在我的实际工作中,我主要关注x86_64和aarch64,并使用默认目标标志和本地目标标志。但我很高兴看到人们讨论不同的情况。我不希望这个问题过于具体,以至于对其他查看此页面的人没有帮助。 - QuarticCat
1
@PeterCordes 嗯,我点赞了那条评论。感谢你教我英语。 :) - QuarticCat
显示剩余2条评论
4个回答

24
我认为没有什么比最大化的天真方法更好了。一个尝试是使用这个等式
x + y = min(x, y) + max(x, y)

因此
max(ctz(x), ctz(y)) = ctz(x) + ctz(y) - min(ctz(x), ctz(y))

通过这种方式,我们可以将最大函数简化为已经优化过的最小函数,尽管需要进行一些额外的操作。
以下是不同方法的Rust实现示例:
pub fn naive(x: u64, y: u64) -> u32 {
    x.trailing_zeros().max(y.trailing_zeros())
}

pub fn sum_minus_min(x: u64, y: u64) -> u32 {
    x.trailing_zeros() + y.trailing_zeros() - (x | y).trailing_zeros()
}

pub fn nielsen(x: u64, y: u64) -> u32 {
    let x_lsb = x & x.wrapping_neg();
    let y_lsb = y & y.wrapping_neg();
    let xy_lsb = x_lsb | y_lsb;
    let lsb = xy_lsb & xy_lsb.wrapping_neg();
    let xy_max_lsb = if xy_lsb == lsb { lsb } else { xy_lsb ^ lsb };
    xy_max_lsb.trailing_zeros()
}

pub fn timmermans(x: u64, y: u64) -> u32 {
    let loxs = !x & x.wrapping_sub(1);
    let loys = !y & y.wrapping_sub(1);
    return (loxs | loys).count_ones();
}

pub fn kealey(x: u64, y: u64) -> u32 {
    ((x | x.wrapping_neg()) & (y | y.wrapping_neg())).trailing_zeros()
}

我的机器上的结果:
ctz_max/naive           time:   [279.09 ns 279.55 ns 280.10 ns]
ctz_max/sum_minus_min   time:   [738.91 ns 742.87 ns 748.61 ns]
ctz_max/nielsen         time:   [935.35 ns 937.63 ns 940.40 ns]
ctz_max/timmermans      time:   [803.39 ns 806.98 ns 810.76 ns]
ctz_max/kealey          time:   [295.03 ns 295.93 ns 297.03 ns]

天真的实现击败了所有其他实现。唯一能与天真实现竞争的是马丁·基利提出的方法。请注意,由于测试工具的一些开销,实际因素可能比时间指示的要高。
很明显,你只有几个CPU指令可以用来优化天真实现,所以我认为你无法做任何事情。作为参考,这里是Rust编译器在现代x86_64处理器上将这些实现编译为独立函数时生成的汇编代码。
example::naive:
        tzcnt   rcx, rdi
        tzcnt   rax, rsi
        cmp     ecx, eax
        cmova   eax, ecx
        ret

example::sum_minus_min:
        tzcnt   rcx, rdi
        tzcnt   rax, rsi
        add     eax, ecx
        or      rsi, rdi
        tzcnt   rcx, rsi
        sub     eax, ecx
        ret

example::nielsen:
        blsi    rax, rdi
        blsi    rcx, rsi
        or      rcx, rax
        blsi    rax, rcx
        xor     edx, edx
        cmp     rcx, rax
        cmovne  rdx, rcx
        xor     rdx, rax
        tzcnt   rax, rdx
        ret

example::timmermans:
        lea     rax, [rdi - 1]
        andn    rax, rdi, rax
        lea     rcx, [rsi - 1]
        andn    rcx, rsi, rcx
        or      rcx, rax
        xor     eax, eax
        popcnt  rax, rcx
        ret

example::kealey:
        mov     rax, rdi
        neg     rax
        or      rax, rdi
        mov     rcx, rsi
        neg     rcx
        or      rcx, rsi
        and     rcx, rax
        tzcnt   rax, rcx
        ret

在我运行的基准测试中,函数被内联,循环部分展开,并且一些子表达式被提取出内部循环,所以汇编代码看起来比上面的要复杂得多。
为了进行测试,我使用了Criterion。以下是额外的代码:
use criterion::{black_box, criterion_group, criterion_main, Criterion};

const NUMBERS: [u64; 32] = [
    ...
];

fn bench<F>(func: F)
where
    F: Fn(u64, u64) -> u32,
{
    for x in NUMBERS {
        for y in NUMBERS {
            black_box(func(x, y));
        }
    }
}

fn compare(c: &mut Criterion) {
    let mut group = c.benchmark_group("ctz_max");
    group.bench_function("naive", |b| b.iter(|| bench(naive)));
    group.bench_function("sum_minus_min", |b| b.iter(|| bench(sum_minus_min)));
    group.bench_function("nielsen", |b| b.iter(|| bench(nielsen)));
    group.bench_function("timmermans", |b| b.iter(|| bench(timmermans)));
    group.bench_function("kealey", |b| b.iter(|| bench(kealey)));
}

criterion_group!(benches, compare);
criterion_main!(benches);

NUMBERS是使用这段Python代码生成的,目的是尽可能地增加min()函数的分支预测难度:

[
    random.randrange(2 ** 32) * 2 ** random.randrange(32)
    for dummy in range(32)
]

我正在运行基准测试使用的内容。
RUSTFLAGS='-C target-cpu=native -C opt-level=3' cargo bench

在第八代i7处理器(威士忌湖)上。

你可能想累加所有结果的总和,并在不正确时抛出异常,以确保没有重要内容被优化掉。此外,请使用 -O3 以及任何你需要做的事情来启用 Rust 中的内联。 - Matt Timmermans
cargo bench 会自动进行优化构建。默认情况下,使用 -O 选项给 rustc 编译器,这等同于 clang 的 -O2 选项。我也尝试了 -O opt-level=3,这会使原始实现减少 5%,并使所有其他版本提高了 5%。我使用 black_box() 避免函数返回值被优化掉。如果我删除 black_box(),整个代码都会被优化掉,所有时间记录都是精确的 0。在优化构建中,内联会自动发生,并且我已经验证了汇编以确保这些函数实际上得到了内联。 - Sven Marnach
很不幸,Rustc/LLVM选择了cmova,它需要2个微操作(因为它需要4个输入,包括CF和ZF的SPAZO组),而没有选择在Broadwell及其后续版本(包括Skylake系列)上只需要1个微操作的cmovbcmovae。(它们只需要CF。)在AMD CPU或Skylake或更高版本中,特别是在tzcnt没有错误依赖关系的情况下,要比2x tzcnt / cmp/cmov更难。在Intel上,它的每时钟吞吐量几乎肯定是可以接受的。 - Peter Cordes
鉴于时间的变化和LLVM对错误依赖关系的普遍不负责任(除非完全看到包含错误依赖关系的循环,否则不愿意花费uops进行xor零操作),在某些测试中可能会出现tzcnt延迟而不是吞吐量的瓶颈?但是不,你的Whiskey Lake CPU没有tzcnt错误依赖关系,所以这不可能是原因。 - Peter Cordes
1
@PeterCordes 实际的基准测试时间相当嘈杂,将函数完全内联到基准测试循环中的汇编非常复杂且难以理解。仅从孤立函数的机器代码来看,无法解释我观察到的时间计算,并且计时会因多种因素而变化,例如函数是否在同一 crate 中定义,即使它们被内联。然而,有一个结果是一致的:无论我做了什么,原始实现在我的机器上都是最快的。 - Sven Marnach
显示剩余5条评论

18

以下三者等效:

  • max(ctz(a),ctz(b))
  • ctz((a|-a)&(b|-b))
  • ctz(a)+ctz(b)-ctz(a|b)

数学公式 ctz(a)+ctz(b)-ctz(a|b) 需要6个CPU指令,可以在3路超标量CPU上并行化为3个步骤:

  • 3× ctz
  • 1× 位或
  • 1× 加法
  • 1× 减法

位运算 ctz((a|-a)&(b|-b)) 需要6个CPU指令,可以在2路超标量CPU上并行化为4个步骤:

  • 2× 取反
  • 2× 位或
  • 1× 位与
  • 1× ctz

朴素的 max(ctz(a),ctz(b)) 需要5个CPU指令,可以在2路超标量CPU上并行化为4个步骤:

  • 2× ctz
  • 1× comparison
  • 1× conditional branch
  • 1× load/move (so that the "output" is always in the same register)

...但请注意,分支指令可能非常昂贵。

如果您的CPU具有条件加载/移动指令,则可以将其减少为4个CPU指令,需要3个超标量步骤。

如果您的CPU具有max指令(例如SSE4),则可以将其减少为3个CPU指令,需要2个超标量步骤。

所有这些都说了,超标量操作的机会取决于您要放置在一起的指令。通常,通过并行放置不同的指令,可以获得最多的效果,因为它们同时使用CPU的不同部分。通常会有更多的“加”和“按位或”单元,而不是“ctz”单元,因此进行多个ctz指令实际上可能是限制因素,特别是对于“数学恒等式”版本。

如果“比较和分支”太昂贵,则可以在4个CPU指令中进行非分支“max”。假设A和B是正整数:

  1. 将B从A中减去,得到C
  2. 将前一次计算时的进位从D中减去,并加上D本身(不管D原先是多少,现在它要么是0要么是-1)
  3. C &= D (C现在为min(0,A-B))
  4. A -= C (A'现在为max(A,B))

还要注意的是,每个带有SSE4的CPU都具有本地64位整数的最大指令。 - Sven Marnach
1
第二个选项在Haswell和Skylake上与默认编译标志(即没有tzcnt)的naive选项相当,根据llvm-mca https://godbolt.org/z/a81ceGWPc。尽管llvm-mca显示naive选项的指令成本略低,但这是因为它无法预测分支成本。我相信这是我们能够达到的最远的地方,所以我会接受这个答案。使用`tzcnt`,也许没有代码可以击败naive选项。 - QuarticCat
而且,我很高兴学到了一些位运算技巧,这是我之前不知道的(a | -a)。 - QuarticCat
我同意这个回答。但是,当操作数比CPU上的单条指令处理的宽度更宽时,选项之间如何进行比较还让我有些好奇。例如,在8位CPU上使用32位操作数。 - nielsen
2
请注意,非分支最大值通常使用条件移动来实现,例如 x86_64 上的 cmov - Sven Marnach
显示剩余6条评论

11
你可以这样做:
#include <algorithm>
#include <bit>
#include <cstdint>

int32_t maxr_zero(uint64_t x, uint64_t y) {
    uint64_t loxs = ~x & (x-1); // low zeros of x
    uint64_t loys = ~y & (y-1); // low zeros of y
    return std::countr_zero((loxs|loys)+1);
}

6
即使是这样简单的东西,也会使用太多的CPU指令来与朴素实现竞争。在现代CPU上,CTZ是一条单一且快速的机器指令,因此朴素实现真的很难被超越。 - Sven Marnach
2
我对这个问题进行了 Rust 版本的基准测试,结果比朴素实现要慢得多。 - Sven Marnach
1
GCC和Clang都使用cmov来实现max(但是GCC还会疯狂地重新引入一个多余的分支来测试y是否为零,以及一个多余的test\cmov对来测试x是否为零)。 - harold
1
哦,对了。我不习惯思考x86汇编语言。使用cmov进行max的天真版本可能会更快。 - Matt Timmermans
1
我认为您可以通过使用 std::popcount(loxs | loys) 稍微改进一下。这样只省去了一个加法操作,但是这也算是些许改进。 - harold
显示剩余12条评论

1

我不确定这个函数是否更快,但是它会接受 xy 并计算输入到 ctz 获取最大值:

uint64_t getMaxTzInput(uint64_t x, uint64_t y)
{
   uint64_t x_lsb = x & (~x + 1);  // Least significant 1 of x
   uint64_t y_lsb = y & (~y + 1);  // Least significant 1 of y
   uint64_t xy_lsb = x_lsb | y_lsb;  // Least significant 1s of x and y (could be the same)
   uint64_t lsb = (xy_lsb) & (~(xy_lsb)+1);  // Least significant 1 among x and y

   // If the least significant 1s are different for x and y, remove the least significant 1
   // to get the second least significant 1.
   uint64_t xy_max_lsb = (xy_lsb == lsb) ? lsb : xy_lsb ^ lsb;
   return xy_max_lsb;
}

因此,ctz(getMaxTzInput(x,y))应该只需调用一次ctz即可给出正确的值。

1
...并且它通过了Marek的单元测试 - Ted Lyngmo
...而且它也通过了我增强版的Marek的单元测试,其中包括案例{0, 0, 64},并检查UB(我的解决方案失败了)。 - Ted Lyngmo
1
但它仍然比朴素实现要慢得多且更加复杂。(我使用 Rust 版本的此代码进行了测量。) - Sven Marnach
4
请注意,(~x + 1) 只是写作 -x 的花式方式。 - Sven Marnach
你的代码假设这两个值都不为零。如果按照朴素的方式,max_ctz(2,0)应该是64,但是你的函数返回2,所以ctz(2)==1。但对于非零输入的情况,我们可以简化最后一步。lsb = xy_lsb & (xy_lsb - 1);(清除最低位) return lsb ? lsb : xy_lsb。如果清除OR结果的最低位得到了零,则这些位在同一位置,因此返回执行此操作之前的值。即使用来自and或blsr的标志进行cmov或csel。(5条指令与x86 BMI1相比,你的8条指令或AArch64相比,8条指令与10条指令:https://godbolt.org/z/73j7xzedf) - Peter Cordes
63-clz(xy_lsb)也可以工作,如果我们知道结果将在0..63范围内,对于非零输入,可以使用clz(xy_lsb) ^ 63。因此,在x、y、orlzcntxor reg,63上使用BMI1 blsi。(https://godbolt.org/z/5ePPzdfh6) 或者只需使用bsr(返回最高位设置的位索引,如果非零)而不是lzcnt(返回前导零计数),因为我们只关心非零输入,如果我们能说服编译器在可用时发出bsr。但是,bsr在AMD上很慢。我知道XOR技巧来自GCC使用bsr / xor来实现countr_zero / __builtin_clz - Peter Cordes

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接