c++ - 用192/256位整数求和无符号64位整数 vector 的点积的最快方法？

Question

我需要计算两个向量的点积：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 <= 60，我们甚至可以保留160（64 * 2 + 32）位整数的结果。

慢速解决方案：


手动处理溢出：

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，然后是adc和adc从低位到高位。

但是我不想求助于ASM，因为它很难维护并且不能移植。

我的目标是使其快速且可维护。

编辑：彼得在评论中指出，当gcc仍使用分支（手动处理溢出）时，clang支持我使用adc的想法。

Answer 1

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



是的，为此的最佳x86-64 asm是mov-从阵列上装入mul或mulx，并从另一个阵列装入内存，然后add + adc装入ans，然后adc reg, 0累积进位。

（也许有一些循环展开，可能有2个或更多的累加器（一组3个寄存器）。如果为Intel CPU展开，则可能避免使用mul内存源的索引寻址模式，因此它可以对负载进行微熔丝。GCC不幸的是，/ clang不创建相对于另一个数组索引一个数组的循环，以最大程度地减少循环开销（1指针增量），同时还避免了其中一个数组的索引寻址模式，因此您不会从一个数组中获取最佳的asm。编译器。但是clang很好。）

clang9.0使用-O3 -march=broadwell -fno-unroll-loops为您的代码生成代码。 Clang识别出通常的a<b习惯用法，即使没有128位整数也可以执行无符号a+=b，导致add reg,reg / adc reg,reg / adc reg,0（不幸的是，clang的循环展开取消优化为mov ； setc； movzx； add而不是adc reg,0，这不利于循环展开的好处！这是应该报告的错过优化的错误。）

GCC实际上在您的if()上分支，您可以通过将其无分支地编写为
overflow += sum<prod;。但这对GCC的其他主要优化遗漏无济于事：实际上对cmp/sbb与sum < prod进行了128位比较，而不是仅使用最后一个adc的CF输出。 GCC知道如何优化uint64_t（Efficient 128-bit addition using carry flag）的方法，但显然不适合__uint128_t的进位。

可能无法通过更改源代码来保持gcc来避免错过优化，这是GCC开发人员必须在GCC中修复的错误。 （因此，您应该在其Bugzilla https://gcc.gnu.org/bugzilla/上将其报告为未优化的错误；并包括指向此问答的链接）。自己尝试完全手动使用uint64_t块会更糟：中间的块具有随身携带和随身携带的功能。这很难用C正确编写，并且GCC不会将其优化为adc。

即使使用overflow += __builtin_add_overflow(sum, prod, &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即可提高效率。乱序的执行程序可以交错多个短FLAGS依赖链。

您可以将另外3个寄存器用于另一个192位累加器（总和和进位），但这可能是不必要的。

对于前端，这看起来像9块（假设mul不能保持微寻址（分层）的索引寻址模式，但是cmp/jcc会进行宏熔断），因此每2.25最多运行1乘时钟周期。 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输入是短的非循环执行的依赖项链的一部分，乱序执行将在早餐时吃掉。

假设它仍与mul的内存操作数分层，那么具有5个宽度管道的Ice Lake可以更快地运行，例如每个迭代1.8个周期。

Zen的流水线宽度为5条指令，如果包含任何多uop指令，则流水线宽度为6 uop。因此，它可能会以每次迭代2个周期运行，这是其2c吞吐量的瓶颈，无法实现完全乘法。 （与Intel相同，正常的imul r64, r64非扩展乘法是1 /时钟，Zen上的延迟为3c。）

Zen 2将mul m64和mulx加速到1 / clock（https://www.uops.info/html-tp/ZEN2/MUL_M64-Measurements.html），并且可能是逐个时钟执行此循环的最快CPU。



通过对一个数组进行相对于另一个数组的展开和索引编制，手动优化版本的成本可能是每乘6微克=移动负载（1微克）+ mul mem（2微克）+加+ 2x adc（3 ），因此在Broadwell上约为1.5个周期/乘法。

它仍然会在前端成为瓶颈，最小的2 uop循环开销（指针增量和cmp / jcc）意味着展开为4可以使我们每6.5个周期乘以4，而不是每6个周期乘以92，即完全展开的理论最大值（假设没有用于内存或次优调度的停顿，或者至少是无序的exec足以吸收内存并且不会使前端停顿。后端应该能够跟上前端的步伐在Haswell及更高版本上，因此ROB应该有空间吸收一些小气泡，但可能不会出现L3或TLB遗漏。）



我认为在这里使用SIMD没有任何好处。添加元素的最大宽度是64位，甚至AVX512也只有64x64 => 128位乘法。除非您可以转换为double，否则AVX512可以以更快的速度运行。