使用 192/256 位整数对无符号 64 位整数向量的点积求和的最快方法?

Posted

技术标签:

【中文标题】使用 192/256 位整数对无符号 64 位整数向量的点积求和的最快方法?【英文标题】:Fastest way to sum dot product of vector of unsigned 64 bit integers using 192/256 bit integer? 【发布时间】:2020-01-03 08:38:23 【问题描述】:

我需要计算两个向量的点积:uint64_t a[N], b[N]; (N<=60) 包含 64 位无符号整数。正是这个循环:

unsigned __int128 ans = 0;
for(int i=0;i<N;++i)
    ans += (unsigned __int128) a[i] * (unsigned __int128) b[i];

ans 会溢出,因此结果必须保存在一个宽整数中,比如 256 位。但是由于 N

缓慢的解决方案:

    手动处理溢出:
unsigned __int128 ans = 0;
uint64_t overflow = 0;
for(int i=0;i<N;++i)
    auto temp = ans;
    ans += (unsigned __int128) a[i] * (unsigned __int128) b[i];
    if(ans<temp) overflow++;

这很慢,因为添加 if 会使程序减慢 ~ 2.2 倍。

    使用 boost::multiprecision::uint256_t 或 GMP 之类的库。

可能很快的解决方案: 如果我们在 64 位机器上使用汇编编程,则可以使用 3 个 64 位寄存器执行加法,方法是使用 add 后跟 adcadc 从低位到高位。

但我不想求助于 ASM,因为它很难维护,而且不可移植。

我的目标是让它快速且可维护。

编辑:Peter 在他的评论中指出,clang 支持我使用adc 的想法,而 gcc 仍然使用分支(手动处理溢出)。

【问题讨论】:

今天在另一个线程上提出了一个稍微相关的讨论,涉及点积计算。 this question 上的 OP 有大约六个实现,他进行了基准测试。在他的问题中也有一个指向他的源代码的链接。 您可以通过简单地屏蔽结果和右移来检测溢出。你不需要if 声明。 @selbie:另一个问题是关于uint8_t 点积,它可以通过转换为float 使用SIMD 自动矢量化。甚至 AVX512 也没有 64x64 => 128 位整数乘法或压缩 128 位整数加法。并且没有内置类型可以容纳两个 __int128 的总和而不会溢出,否则让编译器生成一个 3x adc 链会很容易。 clang 识别通常的 sum &lt; x 习惯用法,用于执行无符号 x+y,即使是 128 位整数。 godbolt.org/z/nZyQcD。 GCC 对if()overflow += sum&lt;prod; 有重大遗漏优化,但至少后者避免了分支。 @user8469759:嗯? ans 的上半部分的低位不会告诉你任何事情。 OP 将有效的 128 位产品累加到一个 160 或 192 位累加器中,由 128 位 ans 和一个单独的高块组成。 ans&gt;&gt;64 会告诉你总和是否溢出了低 64 位。 【参考方案1】:

您绝对不需要 256 位累加器。 N 整数的总和只能产生 大多数 N 进位,因此 64 位 overflow 进位计数器足以满足 2^64 个元素的向量的点积长。即 128 位乘积之和的总宽度只需 192 = 64*3 或 160 = 128 + 32 位。即一个额外的寄存器。


是的,最佳的 x86-64 asm 是 mov-从阵列加载,mulmulx 使用来自另一个阵列的内存源,然后将 add + adc 加载到 @987654337 @, 和adc reg, 0 累积结转。

(可能有一些循环展开,可能有 2 个或更多累加器(3 个寄存器组)。如果为 Intel CPU 展开,可能避免mul 内存源的索引寻址模式,因此它可以微融合不幸的是,GCC / clang 不会创建相对于另一个数组索引一个数组的循环,以最小化循环开销(1 个指针增量),同时还避免其中一个数组的索引寻址模式,因此您不会获得最佳效果来自编译器的 asm。但是 clang 非常好。)

clang9.0 使用 -O3 -march=broadwell -fno-unroll-loops 为您的代码生成它。 Clang 识别通常的 a&lt;b 习惯用法,用于执行无符号 a+=b,即使使用 128 位整数,导致 @ 987654343@ / adc reg,reg / adc reg,0 (不幸的是,clang 的循环展开优化为 mov ; setc ; movzx ; add 而不是 adc reg,0,破坏了循环展开的好处!应该报告的错过优化错误。)

GCC 实际上在您的if() 上分支,您可以通过将其无分支地编写为overflow += sum&lt;prod; 来修复它。但这对 GCC 的其他主要遗漏优化没有帮助:实际上对 sum &lt; prodcmp/sbb 进行 128 位比较,而不是仅使用最后一个 adc 的 CF 输出。 GCC 知道如何为 uint64_t (Efficient 128-bit addition using carry flag) 优化它,但显然不适合从 __uint128_t 执行。

可能无法通过更改源来手动控制 gcc 以避免错过优化,这是 GCC 开发人员必须在 GCC 中修复的错误。(因此您应该将其报告为他们的 bugzilla https://gcc.gnu.org/bugzilla/ 上的一个错过的优化错误;包括指向此 Q&A 的链接)。尝试自己使用 uint64_t 块进行全手动操作会更糟:中间块有进位和进位。这很难在 C 中正确编写,而且 GCC 不会将其优化为 adc

即使使用 overflow += __builtin_add_overflow(sum, prod, &amp;sum); 也没有完全帮助,给我们同样的 cmp/sbb/adc GCC 代码生成! https://gcc.gnu.org/onlinedocs/gcc/Integer-Overflow-Builtins.html 说“编译器将尝试使用硬件指令在可能的情况下实现这些内置函数,例如加法后溢出条件跳转,进位条件跳转等。“但显然它只是不知道如何处理 128 位整数。

#include <stdint.h>

static const int N = 2048;
using u128 = unsigned __int128;
u128 dp(uint64_t a[], uint64_t b[], uint64_t *overflow_out)

    u128 sum = 0;
    uint64_t overflow = 0;
    for(int i=0;i<N;++i)
        u128 prod = a[i] * (unsigned __int128) b[i];
        //sum += prod;
        //overflow += sum<prod;       // gcc uses adc after a 3x mov / cmp / sbb; clang is still good
        overflow += __builtin_add_overflow(sum, prod, &sum);  // gcc less bad, clang still good
    

    *overflow_out = overflow;
    return sum;

Clang 很好地编译了这个 (Godbolt):

# clang9.0 -O3 -march=broadwell -fno-unroll-loops  -Wall -Wextra
     ... zero some regs, copy RDX to R8 (because MUL uses that implicitly)

.LBB0_1:                                # =>This Inner Loop Header: Depth=1
        mov     rax, qword ptr [rsi + 8*rcx]
        mul     qword ptr [rdi + 8*rcx]
        add     r10, rax
        adc     r9, rdx                 # sum += prod
        adc     r11, 0                  # overflow += carry_out
        inc     rcx
        cmp     rcx, 2048
        jne     .LBB0_1               # while(i < N);

        mov     qword ptr [r8], r11   # store the carry count
        mov     rax, r10
        mov     rdx, r9               # return sum in RDX:RAX
        ret

请注意,您不需要需要 ADOX / ADCX 或 MULX 来提高效率。乱序 exec 可以交错多个短 FLAGS 依赖链。

您可以将另外 3 个寄存器用于另一个 192 位累加器(求和和进位),但这可能没有必要。

这看起来像前端的 9 微指令(假设 mul 不能保持索引寻址模式微熔(非层压),但 cmp/jcc 将宏熔)所以它最多只能运行 1每 2.25 个时钟周期相乘。 Haswell 和更早的版本将 adc reg,reg 运行为 2 uop,但 Broadwell 将 3 输入指令(FLAGS + 2 regs)改进为 1 uop。 adc reg,0 在 SnB 系列上特例为 1 uop。

循环携带的依赖链每个只有 1 个循环:add 到 r10,adc 到 r9,adc 到 r11。这些 ADC 指令的 FLAGS 输入是一个短的非循环携带依赖链的一部分,乱序执行将在早餐时吃掉它。

具有 5 宽管道的 Ice Lake 运行速度会更快一些,例如每次迭代可能需要 1.8 个周期,假设它仍然未分层 mul 的内存操作数。

Zen 有一条 5 条指令宽的流水线,如果它包含任何多 uop 指令,则为 6 uop 宽。因此,它可能会在每次迭代 2 个周期时运行它,在其 2c 吞吐量方面遇到瓶颈以进行完全乘法。 (正常imul r64, r64 非扩展乘法是 1/clock,Zen 上的延迟为 3c,与 Intel 相同。)

Zen 2 将 mul m64mulx 加速到 1 个/时钟 (https://www.uops.info/html-tp/ZEN2/MUL_M64-Measurements.html),并且在时钟换时钟的基础上可能是该循环中最快的 CPU。


通过相对于另一个数组展开和索引一个数组,手动优化版本的每次乘法成本可能接近 6 uop = mov-load (1 uop) + mul mem ( 2 uops)+ add + 2x adc(3 uops),所以在 Broadwell 上大约 1.5 个周期/乘法。

它仍然会成为前端的瓶颈,并且 2 uop(指针增量和 cmp/jcc)的最小循环开销意味着展开 4 可以为我们每 6.5 个周期而不是每 6 个周期提供 4 次乘法,或者完全展开理论最大值的 92%。 (假设内存或次优调度没有停顿,或者至少无序 exec 足以吸收这一点并且不会停顿前端。后端应该能够跟上前端在这里,在 Haswell 和更高版本上,所以 ROB 应该有空间来吸收一些小气泡,但可能不是 L3 或 TLB 未命中。)


我认为在这里使用 SIMD 没有任何好处;加法的最宽元素大小是 64 位,甚至 AVX512 也只有 64x64 => 128 位乘法。除非您可以转换为 double,在这种情况下,AVX512 可以非常更快地运行。

【讨论】:

我在 GCC 中提交了一个关于错过优化的错误:gcc.gnu.org/bugzilla/show_bug.cgi?id=93141 你能解释一下这个说法吗:gcc doesn't actually *branch* unless you use an if()。 AFAIK 分支是一个比较,可以跟随一个跳转。参考:参见条件跳转:汇编中的分支cs.uaf.edu/2012/fall/cs301/lecture/09_12_jmp.html部分 我还在上大学,还有很多关于汇编编程的知识:/ 感谢您的评论。另外,我认为您也应该在 clang 中提交错误。我可能会错过细节。但请在此处分享 BUG ID。 @madhur4127:不,分支是条件跳转指令(或任何跳转,取决于上下文)。特别是JCC instruction。 compare-and-branch 是一种常见的模式,但例如你可以 add ; jc 从添加执行分支。无论如何,cmp 不是分支;它对 RIP(程序计数器)没有任何奇怪的作用。它只写标志。在我们的例子中,使用这些标志的指令是 SBB,而不是 JCC。 @BrettHale:我最初对 MULX 是否会有所帮助也有同样的想法,但后来我意识到通过 FLAGS 的 dep 链很短且独立。 MULX + ADOX/CX 适用于一个大乘法,其中所有加法都是一个长链进位的一部分(或 2 个带有 ADX)。但是在这里,每个进位链都在 add/adc/adc 0 之后结束,非常短,足以让乱序 exec 在几次乘法迭代中为我们交错。如果您正在为有序 CPU 进行软件流水线处理,手动交错内容,则只需要 ADCX/ADOX。

以上是关于使用 192/256 位整数对无符号 64 位整数向量的点积求和的最快方法?的主要内容,如果未能解决你的问题,请参考以下文章

使用 32 位无符号整数乘以 64 位数字的算法

如何在 64 位 Perl 中解压缩(64 位)无符号长整数?

如何在 IA32 上将带符号的整数相加成更大的和。 32位有符号整数的64位总和?

将 8 个字节转换为有符号长整数(64 位)

NSDecimalNumber 和大的无符号长长(64 位)整数

在 32 位 iPhone 上处理 64 位整数