为什么 _umul128 比 mul128x64x2 函数的标量代码运行得慢？
Why _umul128 works slower than scalar code for mul128x64x2 function?
我第二次尝试实现快速 mul128x64x2 函数。 First time I ask the question 不与 _umul128 MSVC 版本进行比较。现在我做了这样的比较，我得到的结果表明 _umul128 函数比本地标量和手工制作的 simd AVX 1.0 代码慢。
下面是我的测试代码：
#include <iostream>
#include <chrono>
#include <intrin.h>
#include <emmintrin.h>
#include <immintrin.h>
#pragma intrinsic(_umul128)
constexpr uint32_t LOW[4] = { 4294967295u, 0u, 4294967295u, 0u };
__forceinline void multiply128x128( const uint32_t ABCD[4], const uint32_t EFGH[4], uint32_t OUT[2][4] ) noexcept
{
__m128i L = _mm_lddqu_si128( reinterpret_cast< __m128i const* >( LOW ) );
__m128i IN = _mm_lddqu_si128( reinterpret_cast< __m128i const* >( EFGH ) );
__m128i A = _mm_set1_epi32( ABCD[0] );
__m128i B = _mm_set1_epi32( ABCD[1] );
__m128i C = _mm_set1_epi32( ABCD[2] );
__m128i D = _mm_set1_epi32( ABCD[3] );
__m128i ED = _mm_mul_epu32( IN, D );
__m128i EC = _mm_mul_epu32( IN, C );
__m128i EB = _mm_mul_epu32( IN, B );
__m128i EA = _mm_mul_epu32( IN, A );
IN = _mm_srli_epi64( IN, 32 );
__m128i FD = _mm_mul_epu32( IN, D );
__m128i FC = _mm_mul_epu32( IN, C );
__m128i FB = _mm_mul_epu32( IN, B );
__m128i FA = _mm_mul_epu32( IN, A );
__m128i FD_H = _mm_srli_epi64( FD, 32 );
__m128i FD_L = _mm_and_si128 ( L, FD );
__m128i FC_H = _mm_srli_epi64( FC, 32 );
__m128i FC_L = _mm_and_si128 ( L, FC );
__m128i FB_H = _mm_srli_epi64( FB, 32 );
__m128i FB_L = _mm_and_si128 ( L, FB );
__m128i FA_H = _mm_srli_epi64( FA, 32 );
__m128i FA_L = _mm_and_si128 ( L, FA );
__m128i ED_H = _mm_srli_epi64( ED, 32 );
__m128i ED_L = _mm_and_si128 ( L, ED );
__m128i EC_H = _mm_srli_epi64( EC, 32 );
__m128i EC_L = _mm_and_si128 ( L, EC );
__m128i EB_H = _mm_srli_epi64( EB, 32 );
__m128i EB_L = _mm_and_si128 ( L, EB );
__m128i EA_H = _mm_srli_epi64( EA, 32 );
__m128i EA_L = _mm_and_si128 ( L, EA );
__m128i SUM_FC_L_FD_H = _mm_add_epi64( FC_L, FD_H );
__m128i SUM_FB_L_FC_H = _mm_add_epi64( FB_L, FC_H );
__m128i SUM_FA_L_FB_H = _mm_add_epi64( FA_L, FB_H );
__m128i SUM_EC_L_ED_H = _mm_add_epi64( EC_L, ED_H );
__m128i SUM_EB_L_EC_H = _mm_add_epi64( EB_L, EC_H );
__m128i SUM_EA_L_EB_H = _mm_add_epi64( EA_L, EB_H );
__m128i SUM_FC_L_FD_H_ED_L = _mm_add_epi64( SUM_FC_L_FD_H, ED_L );
__m128i SUM_FB_L_FC_H_EC_L_ED_H = _mm_add_epi64( SUM_FB_L_FC_H, SUM_EC_L_ED_H );
__m128i SUM_FA_L_FB_H_EB_L_EC_H = _mm_add_epi64( SUM_FA_L_FB_H, SUM_EB_L_EC_H );
__m128i SUM_FA_H_EA_L_EB_H = _mm_add_epi64( FA_H, SUM_EA_L_EB_H );
__m128i SUM_FC_L_FD_H_ED_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L, 32 );
SUM_FC_L_FD_H_ED_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L, SUM_FB_L_FC_H_EC_L_ED_H );
__m128i SUM_FC_L_FD_H_ED_L_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L_L, 32 );
SUM_FC_L_FD_H_ED_L_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L_L, SUM_FA_L_FB_H_EB_L_EC_H );
__m128i SUM_FC_L_FD_H_ED_L_L_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L_L_L, 32 );
SUM_FC_L_FD_H_ED_L_L_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L_L_L, SUM_FA_H_EA_L_EB_H );
__m128i SUM_FC_L_FD_H_ED_L_L_L_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L_L_L_L, 32 );
SUM_FC_L_FD_H_ED_L_L_L_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L_L_L_L, EA_H );
OUT[0][0] = SUM_FC_L_FD_H_ED_L_L_L_L_L.m128i_u32[0];
OUT[0][1] = SUM_FC_L_FD_H_ED_L_L_L_L.m128i_u32[0];
OUT[0][2] = SUM_FC_L_FD_H_ED_L_L_L.m128i_u32[0];
OUT[0][3] = SUM_FC_L_FD_H_ED_L_L.m128i_u32[0];
OUT[1][0] = SUM_FC_L_FD_H_ED_L_L_L_L_L.m128i_u32[2];
OUT[1][1] = SUM_FC_L_FD_H_ED_L_L_L_L.m128i_u32[2];
OUT[1][2] = SUM_FC_L_FD_H_ED_L_L_L.m128i_u32[2];
OUT[1][3] = SUM_FC_L_FD_H_ED_L_L.m128i_u32[2];
}
__forceinline void multiply128x128_1( const uint32_t ABCD[4], const uint32_t EFGH[4], uint32_t OUT[2][4] ) noexcept
{
uint64_t ED = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t EC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t EB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t EA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t FD = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t FC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t FB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t FA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t GD = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t GC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t GB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t GA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t HD = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t HC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t HB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t HA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t SUM_FC_L_FD_H = ( FC & 0xFFFFFFFF ) + ( FD >> 32u );
uint64_t SUM_FB_L_FC_H = ( FB & 0xFFFFFFFF ) + ( FC >> 32u );
uint64_t SUM_FA_L_FB_H = ( FA & 0xFFFFFFFF ) + ( FB >> 32u );
uint64_t SUM_EC_L_ED_H = ( EC & 0xFFFFFFFF ) + ( ED >> 32u );
uint64_t SUM_EB_L_EC_H = ( EB & 0xFFFFFFFF ) + ( EC >> 32u );
uint64_t SUM_EA_L_EB_H = ( EA & 0xFFFFFFFF ) + ( EB >> 32u );
uint64_t SUM_HC_L_HD_H = ( HC & 0xFFFFFFFF ) + ( HD >> 32u );
uint64_t SUM_HB_L_HC_H = ( HB & 0xFFFFFFFF ) + ( HC >> 32u );
uint64_t SUM_HA_L_HB_H = ( HA & 0xFFFFFFFF ) + ( HB >> 32u );
uint64_t SUM_GC_L_GD_H = ( GC & 0xFFFFFFFF ) + ( GD >> 32u );
uint64_t SUM_GB_L_GC_H = ( GB & 0xFFFFFFFF ) + ( GC >> 32u );
uint64_t SUM_GA_L_GB_H = ( GA & 0xFFFFFFFF ) + ( GB >> 32u );
uint64_t SUM_FC_L_FD_H_ED_L = SUM_FC_L_FD_H + ( ED & 0xFFFFFFFF );
uint64_t SUM_FB_L_FC_H_EC_L_ED_H = SUM_FB_L_FC_H + SUM_EC_L_ED_H;
uint64_t SUM_FA_L_FB_H_EB_L_EC_H = SUM_FA_L_FB_H + SUM_EB_L_EC_H;
uint64_t SUM_FA_H_EA_L_EB_H = SUM_EA_L_EB_H + ( FA >> 32u );
uint64_t SUM_FC_L_FD_H_ED_L_L = ( SUM_FC_L_FD_H_ED_L >> 32u ) + SUM_FB_L_FC_H_EC_L_ED_H;
uint64_t SUM_FC_L_FD_H_ED_L_L_L = ( SUM_FC_L_FD_H_ED_L_L >> 32u ) + SUM_FA_L_FB_H_EB_L_EC_H;
uint64_t SUM_FC_L_FD_H_ED_L_L_L_L = ( SUM_FC_L_FD_H_ED_L_L_L >> 32u ) + SUM_FA_H_EA_L_EB_H;
uint64_t SUM_FC_L_FD_H_ED_L_L_L_L_L = ( SUM_FC_L_FD_H_ED_L_L_L_L >> 32u ) + ( EA >> 32u );
uint64_t SUM_HC_L_HD_H_GD_L = SUM_HC_L_HD_H + ( GD & 0xFFFFFFFF );
uint64_t SUM_HB_L_HC_H_GC_L_GD_H = SUM_HB_L_HC_H + SUM_GC_L_GD_H;
uint64_t SUM_HA_L_HB_H_GB_L_GC_H = SUM_HA_L_HB_H + SUM_GB_L_GC_H;
uint64_t SUM_HA_H_GA_L_GB_H = SUM_GA_L_GB_H + ( HA >> 32u );
uint64_t SUM_HC_L_HD_H_GD_L_L = ( SUM_HC_L_HD_H_GD_L >> 32u ) + SUM_HB_L_HC_H_GC_L_GD_H;
uint64_t SUM_HC_L_HD_H_GD_L_L_L = ( SUM_HC_L_HD_H_GD_L_L >> 32u ) + SUM_HA_L_HB_H_GB_L_GC_H;
uint64_t SUM_HC_L_HD_H_GD_L_L_L_L = ( SUM_HC_L_HD_H_GD_L_L_L >> 32u ) + SUM_HA_H_GA_L_GB_H;
uint64_t SUM_HC_L_HD_H_GD_L_L_L_L_L = ( SUM_HC_L_HD_H_GD_L_L_L_L >> 32u ) + ( GA >> 32u );
OUT[0][0] = SUM_FC_L_FD_H_ED_L_L_L_L_L;
OUT[0][1] = SUM_FC_L_FD_H_ED_L_L_L_L;
OUT[0][2] = SUM_FC_L_FD_H_ED_L_L_L;
OUT[0][3] = SUM_FC_L_FD_H_ED_L_L;
OUT[1][0] = SUM_HC_L_HD_H_GD_L_L_L_L_L;
OUT[1][1] = SUM_HC_L_HD_H_GD_L_L_L_L;
OUT[1][2] = SUM_HC_L_HD_H_GD_L_L_L;
OUT[1][3] = SUM_HC_L_HD_H_GD_L_L;
}
__forceinline void mulShift( const uint64_t* const m, const uint64_t* const mul , uint32_t OUT[2][4]) noexcept
{
uint64_t B0[2];
uint64_t B2[2];
{
B0[0] = _umul128( m[1], mul[0], &B0[1] );
B2[0] = _umul128( m[0], mul[0], &B2[1] );
uint64_t S = B0[1] + B2[0];
OUT[0][2] = S >> 32;
OUT[0][3] = S & 0xFFFFFFFF;
uint64_t M = B2[1] + ( S < B2[0] );
OUT[0][1] = M & 0xFFFFFFFF;
OUT[0][0] = M >> 32;
}
{
B0[0] = _umul128( m[1], mul[1], &B0[1] );
B2[0] = _umul128( m[0], mul[1], &B2[1] );
uint64_t S = B0[1] + B2[0];
OUT[1][2] = S >> 32;
OUT[1][3] = S & 0xFFFFFFFF;
uint64_t M = B2[1] + ( S < B2[0] );
OUT[1][1] = M & 0xFFFFFFFF;
OUT[1][0] = M >> 32;
}
}
constexpr uint32_t N = 1 << 28;
int main()
{
uint32_t OUT[2][4];
uint32_t ABCD[4] = { 4294967295u, 4294967295u, 4294967295u, 4294967295u };
uint32_t EFGH[4] = { 4294967295u, 4294967295u, 4294967295u, 4294967295u };
multiply128x128_1( ABCD, EFGH, OUT );
uint64_t S_1 = 0u;
uint64_t S_2 = 0u;
uint64_t S_3 = 0u;
auto start_1 = std::chrono::high_resolution_clock::now();
for ( uint32_t i = 0; i < N; ++i )
{
EFGH[0] = i;
EFGH[1] = i;
EFGH[2] = i + 1;
EFGH[3] = i + 1;
ABCD[0] = i;
ABCD[1] = i;
ABCD[2] = i + 1;
ABCD[3] = i + 1;
multiply128x128( ABCD, EFGH, OUT );
S_1 += OUT[0][0] + OUT[0][1] + OUT[0][2] + OUT[0][3];
S_1 += OUT[1][0] + OUT[1][1] + OUT[1][2] + OUT[1][3];
}
auto stop_1 = std::chrono::high_resolution_clock::now();
std::cout << "Test A: " << std::chrono::duration_cast<std::chrono::milliseconds>( stop_1 - start_1 ).count() << '\n';
auto start_2 = std::chrono::high_resolution_clock::now();
for ( uint32_t i = 0; i < N; ++i )
{
EFGH[0] = i;
EFGH[1] = i;
EFGH[2] = i + 1;
EFGH[3] = i + 1;
ABCD[0] = i;
ABCD[1] = i;
ABCD[2] = i + 1;
ABCD[3] = i + 1;
mulShift( reinterpret_cast<const uint64_t*>( ABCD ), reinterpret_cast<const uint64_t*>( EFGH ), OUT );
S_2 += OUT[0][0] + OUT[0][1] + OUT[0][2] + OUT[0][3];
S_2 += OUT[1][0] + OUT[1][1] + OUT[1][2] + OUT[1][3];
}
auto stop_2 = std::chrono::high_resolution_clock::now();
std::cout << "Test B: " << std::chrono::duration_cast<std::chrono::milliseconds>( stop_2 - start_2 ).count() << '\n';
auto start_3 = std::chrono::high_resolution_clock::now();
for ( uint32_t i = 0; i < N; ++i )
{
EFGH[0] = i;
EFGH[1] = i;
EFGH[2] = i + 1;
EFGH[3] = i + 1;
ABCD[0] = i;
ABCD[1] = i;
ABCD[2] = i + 1;
ABCD[3] = i + 1;
multiply128x128_1( ABCD, EFGH, OUT );
S_3 += OUT[0][0] + OUT[0][1] + OUT[0][2] + OUT[0][3];
S_3 += OUT[1][0] + OUT[1][1] + OUT[1][2] + OUT[1][3];
}
auto stop_3 = std::chrono::high_resolution_clock::now();
std::cout << "Test C: " << std::chrono::duration_cast<std::chrono::milliseconds>( stop_3 - start_3 ).count() << '\n';
std::cout << S_1 << " " << S_2 << " " << S_3 << '\n';
}
为什么 _umul128 这么慢？也许我在上面的测试代码中犯了一些错误？
我的结果：
测试 A (simd)：4546 毫秒。
测试 B (_umul128)：6637 毫秒。
测试 C（标量）：2333 毫秒。
在 Windows 10、x64、MSVC 2019 上测试
_umul128 版本并没有那么慢，但你通过摆弄 32 位数组使 MSVC发出可怕的 asm.
优化正在击败您的基准；纯 C 版本并没有那么快。
尤其是简单的输入数据：
ABCD[0] = EFGH[0] = i;
ABCD[1] = EFGH[1] = i;
ABCD[2] = EFGH[2] = i + 1;
ABCD[3] = EFGH[3] = i + 1;
像这样初始化两个输入在内联纯 C 版本后创造了大量的优化机会。 i*i 4 次，i*(i+1) = i*i + i 8 次，(i+1)*(i+1) 4 次。 MSVC 并不愚蠢并且注意到了这一点。 这叫做 Common Subexpression Elimination (CSE)。
如果你想看看纯 C 到底有多慢，你需要想出一种更复杂的方法来伪造输入。也许提前生成然后遍历包含输入的内存？从循环计数器设置输入的成本几乎与乘法一样多。
MSVC 的 asm 输出证实大部分工作针对纯 C 版本进行了优化。 (Godbolt with MSVC 19.22 for x64)
...
$LL10@main:
lea r15, QWORD PTR [rax+1]
mov rcx, r15
mov r9, r15
imul rcx, rax # only 3, not 16, imul instructions.
imul rax, rax # (None appear later in this loop in the ... part)
imul r9, r15
mov edi, ecx
mov r14, rcx
mov r8d, eax
shr r14, 32 ; 00000020H
shr rax, 32 ; 00000020H
...
sub r13, 1
jne $LL10@main
MSVC 不善于优化内在函数 并且执行所有 4 个 mul m64 指令而不是注意到 ii * i1i1 执行了两次。
更重要的是，_umul128 循环受到存储转发停顿的伤害 因为它实际上 将数组存储到内存中使用 32 位存储，然后使用 64 位加载来提供 mul m64.
此外，处理 32 位块中的输出只会搬起石头砸自己的脚，引入额外的移位和 mov 操作。
这并不复杂，字面上只需要 3 条指令，mul r64 和 imul r64, r64 加上用于高半部分的 add，就足够了。 GCC/clang 轻松发出正确的东西，x86-64 System V 调用约定可以 return 寄存器中的 128 位 int。
关于神马：https://godbolt.org/z/DcZhSl
#include <stdint.h>
#ifdef __GNUC__
typedef unsigned __int128 u128;
u128 mul128x64( u128 a, uint64_t b) {
return a * b;
}
#endif
# clang -O3 for the x86-64 System V ABI (Linux)
mul128x64(unsigned __int128, unsigned long): #
mov rax, rdi
imul rsi, rdx
mul rdx
add rdx, rsi
ret
对于 MSVC，我们必须自己做，调用约定意味着结果 returned 在内存中。
#ifdef _MSC_VER
#include <intrin.h>
struct u128 { uint64_t u64[2]; };
u128 mul128x64( uint64_t a_lo, uint64_t a_hi, uint64_t b)
{
uint64_t lolo_high;
uint64_t lolo = _umul128( a_lo, b, &lolo_high );
uint64_t lohi = a_hi * b;
return {{lolo, lohi + lolo_high}};
}
#endif
# MSVC x64 -O2
u128 mul128x64(unsigned __int64,unsigned __int64,unsigned __int64) PROC
mov rax, r9
mul rdx
imul r8, r9
mov QWORD PTR [rcx], rax # store the retval into hidden pointer
mov rax, rcx
add r8, rdx
mov QWORD PTR [rcx+8], r8
ret 0
您的 __m128i 内在函数版本不太可能成功 。现代 x86（主流 Intel SnB 系列，AMD Ryzen）在 mul 和 imul 上具有 1/clock 吞吐量。 （除了 Ryzen，其中加宽 i/mul r64 的吞吐量为 2c，但对于 imul r64,r64 仍然是 1/clock。）
因此，如果您在 C 中实现并像这样编译成 asm，则 Sandybridge 系列上 64 x 128 位乘法的总吞吐量是每 2 个周期一次（在端口 1 上出现瓶颈）。
鉴于您需要超过 4 pmuludq 条指令来实现乘法运算，AVX1 是行不通的。 （Skylake 的 pmuludq 吞吐量为 0.5c。Sandybridge 的吞吐量为 1c，因此您需要在每次乘法（平均）2 pmuludq insns 内完成工作才能与标量竞争。而且这还没有考虑所有需要做的轮班/洗牌/添加工作。
可能值得考虑 Bulldozer 系列，其中 64 位标量乘法是 4c 吞吐量，但 pmuludq 是 1c。 (https://agner.org/optimize/) 每个周期产生 128 个产品位（两个 32x32 => 64 位产品）比每 4 个周期产生 128 个产品位要好，如果你可以在不消耗太多额外周期的情况下移动和添加它们。
同样，MSVC 不擅长通过内在函数进行持续传播或 CSE 优化，因此您的内在函数版本没有任何好处。
您的测试代码还使用来自标量整数循环变量的 _mm_set1_epi32( )，需要 vmovd 和 vpshufd 指令。
并且您为这些数组上的 lddqu 内在函数重新加载标量存储/向量，因此您再次遇到存储转发停顿。
SSE2 或 AVX1 的唯一希望是你的数据来自内存，而不是寄存器。或者，如果您可以将数据长期保存在矢量寄存器中，而不是不断地来回移动数据。特别是在 int <-> SIMD 具有高延迟的 Bulldozer 系列上。
我第二次尝试实现快速 mul128x64x2 函数。 First time I ask the question 不与 _umul128 MSVC 版本进行比较。现在我做了这样的比较，我得到的结果表明 _umul128 函数比本地标量和手工制作的 simd AVX 1.0 代码慢。
下面是我的测试代码：
#include <iostream>
#include <chrono>
#include <intrin.h>
#include <emmintrin.h>
#include <immintrin.h>
#pragma intrinsic(_umul128)
constexpr uint32_t LOW[4] = { 4294967295u, 0u, 4294967295u, 0u };
__forceinline void multiply128x128( const uint32_t ABCD[4], const uint32_t EFGH[4], uint32_t OUT[2][4] ) noexcept
{
__m128i L = _mm_lddqu_si128( reinterpret_cast< __m128i const* >( LOW ) );
__m128i IN = _mm_lddqu_si128( reinterpret_cast< __m128i const* >( EFGH ) );
__m128i A = _mm_set1_epi32( ABCD[0] );
__m128i B = _mm_set1_epi32( ABCD[1] );
__m128i C = _mm_set1_epi32( ABCD[2] );
__m128i D = _mm_set1_epi32( ABCD[3] );
__m128i ED = _mm_mul_epu32( IN, D );
__m128i EC = _mm_mul_epu32( IN, C );
__m128i EB = _mm_mul_epu32( IN, B );
__m128i EA = _mm_mul_epu32( IN, A );
IN = _mm_srli_epi64( IN, 32 );
__m128i FD = _mm_mul_epu32( IN, D );
__m128i FC = _mm_mul_epu32( IN, C );
__m128i FB = _mm_mul_epu32( IN, B );
__m128i FA = _mm_mul_epu32( IN, A );
__m128i FD_H = _mm_srli_epi64( FD, 32 );
__m128i FD_L = _mm_and_si128 ( L, FD );
__m128i FC_H = _mm_srli_epi64( FC, 32 );
__m128i FC_L = _mm_and_si128 ( L, FC );
__m128i FB_H = _mm_srli_epi64( FB, 32 );
__m128i FB_L = _mm_and_si128 ( L, FB );
__m128i FA_H = _mm_srli_epi64( FA, 32 );
__m128i FA_L = _mm_and_si128 ( L, FA );
__m128i ED_H = _mm_srli_epi64( ED, 32 );
__m128i ED_L = _mm_and_si128 ( L, ED );
__m128i EC_H = _mm_srli_epi64( EC, 32 );
__m128i EC_L = _mm_and_si128 ( L, EC );
__m128i EB_H = _mm_srli_epi64( EB, 32 );
__m128i EB_L = _mm_and_si128 ( L, EB );
__m128i EA_H = _mm_srli_epi64( EA, 32 );
__m128i EA_L = _mm_and_si128 ( L, EA );
__m128i SUM_FC_L_FD_H = _mm_add_epi64( FC_L, FD_H );
__m128i SUM_FB_L_FC_H = _mm_add_epi64( FB_L, FC_H );
__m128i SUM_FA_L_FB_H = _mm_add_epi64( FA_L, FB_H );
__m128i SUM_EC_L_ED_H = _mm_add_epi64( EC_L, ED_H );
__m128i SUM_EB_L_EC_H = _mm_add_epi64( EB_L, EC_H );
__m128i SUM_EA_L_EB_H = _mm_add_epi64( EA_L, EB_H );
__m128i SUM_FC_L_FD_H_ED_L = _mm_add_epi64( SUM_FC_L_FD_H, ED_L );
__m128i SUM_FB_L_FC_H_EC_L_ED_H = _mm_add_epi64( SUM_FB_L_FC_H, SUM_EC_L_ED_H );
__m128i SUM_FA_L_FB_H_EB_L_EC_H = _mm_add_epi64( SUM_FA_L_FB_H, SUM_EB_L_EC_H );
__m128i SUM_FA_H_EA_L_EB_H = _mm_add_epi64( FA_H, SUM_EA_L_EB_H );
__m128i SUM_FC_L_FD_H_ED_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L, 32 );
SUM_FC_L_FD_H_ED_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L, SUM_FB_L_FC_H_EC_L_ED_H );
__m128i SUM_FC_L_FD_H_ED_L_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L_L, 32 );
SUM_FC_L_FD_H_ED_L_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L_L, SUM_FA_L_FB_H_EB_L_EC_H );
__m128i SUM_FC_L_FD_H_ED_L_L_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L_L_L, 32 );
SUM_FC_L_FD_H_ED_L_L_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L_L_L, SUM_FA_H_EA_L_EB_H );
__m128i SUM_FC_L_FD_H_ED_L_L_L_L_L = _mm_srli_epi64( SUM_FC_L_FD_H_ED_L_L_L_L, 32 );
SUM_FC_L_FD_H_ED_L_L_L_L_L = _mm_add_epi64 ( SUM_FC_L_FD_H_ED_L_L_L_L_L, EA_H );
OUT[0][0] = SUM_FC_L_FD_H_ED_L_L_L_L_L.m128i_u32[0];
OUT[0][1] = SUM_FC_L_FD_H_ED_L_L_L_L.m128i_u32[0];
OUT[0][2] = SUM_FC_L_FD_H_ED_L_L_L.m128i_u32[0];
OUT[0][3] = SUM_FC_L_FD_H_ED_L_L.m128i_u32[0];
OUT[1][0] = SUM_FC_L_FD_H_ED_L_L_L_L_L.m128i_u32[2];
OUT[1][1] = SUM_FC_L_FD_H_ED_L_L_L_L.m128i_u32[2];
OUT[1][2] = SUM_FC_L_FD_H_ED_L_L_L.m128i_u32[2];
OUT[1][3] = SUM_FC_L_FD_H_ED_L_L.m128i_u32[2];
}
__forceinline void multiply128x128_1( const uint32_t ABCD[4], const uint32_t EFGH[4], uint32_t OUT[2][4] ) noexcept
{
uint64_t ED = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t EC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t EB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t EA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[0] );
uint64_t FD = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t FC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t FB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t FA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[1] );
uint64_t GD = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t GC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t GB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t GA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[2] );
uint64_t HD = static_cast<uint64_t>( ABCD[3] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t HC = static_cast<uint64_t>( ABCD[2] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t HB = static_cast<uint64_t>( ABCD[1] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t HA = static_cast<uint64_t>( ABCD[0] ) * static_cast<uint64_t>( EFGH[3] );
uint64_t SUM_FC_L_FD_H = ( FC & 0xFFFFFFFF ) + ( FD >> 32u );
uint64_t SUM_FB_L_FC_H = ( FB & 0xFFFFFFFF ) + ( FC >> 32u );
uint64_t SUM_FA_L_FB_H = ( FA & 0xFFFFFFFF ) + ( FB >> 32u );
uint64_t SUM_EC_L_ED_H = ( EC & 0xFFFFFFFF ) + ( ED >> 32u );
uint64_t SUM_EB_L_EC_H = ( EB & 0xFFFFFFFF ) + ( EC >> 32u );
uint64_t SUM_EA_L_EB_H = ( EA & 0xFFFFFFFF ) + ( EB >> 32u );
uint64_t SUM_HC_L_HD_H = ( HC & 0xFFFFFFFF ) + ( HD >> 32u );
uint64_t SUM_HB_L_HC_H = ( HB & 0xFFFFFFFF ) + ( HC >> 32u );
uint64_t SUM_HA_L_HB_H = ( HA & 0xFFFFFFFF ) + ( HB >> 32u );
uint64_t SUM_GC_L_GD_H = ( GC & 0xFFFFFFFF ) + ( GD >> 32u );
uint64_t SUM_GB_L_GC_H = ( GB & 0xFFFFFFFF ) + ( GC >> 32u );
uint64_t SUM_GA_L_GB_H = ( GA & 0xFFFFFFFF ) + ( GB >> 32u );
uint64_t SUM_FC_L_FD_H_ED_L = SUM_FC_L_FD_H + ( ED & 0xFFFFFFFF );
uint64_t SUM_FB_L_FC_H_EC_L_ED_H = SUM_FB_L_FC_H + SUM_EC_L_ED_H;
uint64_t SUM_FA_L_FB_H_EB_L_EC_H = SUM_FA_L_FB_H + SUM_EB_L_EC_H;
uint64_t SUM_FA_H_EA_L_EB_H = SUM_EA_L_EB_H + ( FA >> 32u );
uint64_t SUM_FC_L_FD_H_ED_L_L = ( SUM_FC_L_FD_H_ED_L >> 32u ) + SUM_FB_L_FC_H_EC_L_ED_H;
uint64_t SUM_FC_L_FD_H_ED_L_L_L = ( SUM_FC_L_FD_H_ED_L_L >> 32u ) + SUM_FA_L_FB_H_EB_L_EC_H;
uint64_t SUM_FC_L_FD_H_ED_L_L_L_L = ( SUM_FC_L_FD_H_ED_L_L_L >> 32u ) + SUM_FA_H_EA_L_EB_H;
uint64_t SUM_FC_L_FD_H_ED_L_L_L_L_L = ( SUM_FC_L_FD_H_ED_L_L_L_L >> 32u ) + ( EA >> 32u );
uint64_t SUM_HC_L_HD_H_GD_L = SUM_HC_L_HD_H + ( GD & 0xFFFFFFFF );
uint64_t SUM_HB_L_HC_H_GC_L_GD_H = SUM_HB_L_HC_H + SUM_GC_L_GD_H;
uint64_t SUM_HA_L_HB_H_GB_L_GC_H = SUM_HA_L_HB_H + SUM_GB_L_GC_H;
uint64_t SUM_HA_H_GA_L_GB_H = SUM_GA_L_GB_H + ( HA >> 32u );
uint64_t SUM_HC_L_HD_H_GD_L_L = ( SUM_HC_L_HD_H_GD_L >> 32u ) + SUM_HB_L_HC_H_GC_L_GD_H;
uint64_t SUM_HC_L_HD_H_GD_L_L_L = ( SUM_HC_L_HD_H_GD_L_L >> 32u ) + SUM_HA_L_HB_H_GB_L_GC_H;
uint64_t SUM_HC_L_HD_H_GD_L_L_L_L = ( SUM_HC_L_HD_H_GD_L_L_L >> 32u ) + SUM_HA_H_GA_L_GB_H;
uint64_t SUM_HC_L_HD_H_GD_L_L_L_L_L = ( SUM_HC_L_HD_H_GD_L_L_L_L >> 32u ) + ( GA >> 32u );
OUT[0][0] = SUM_FC_L_FD_H_ED_L_L_L_L_L;
OUT[0][1] = SUM_FC_L_FD_H_ED_L_L_L_L;
OUT[0][2] = SUM_FC_L_FD_H_ED_L_L_L;
OUT[0][3] = SUM_FC_L_FD_H_ED_L_L;
OUT[1][0] = SUM_HC_L_HD_H_GD_L_L_L_L_L;
OUT[1][1] = SUM_HC_L_HD_H_GD_L_L_L_L;
OUT[1][2] = SUM_HC_L_HD_H_GD_L_L_L;
OUT[1][3] = SUM_HC_L_HD_H_GD_L_L;
}
__forceinline void mulShift( const uint64_t* const m, const uint64_t* const mul , uint32_t OUT[2][4]) noexcept
{
uint64_t B0[2];
uint64_t B2[2];
{
B0[0] = _umul128( m[1], mul[0], &B0[1] );
B2[0] = _umul128( m[0], mul[0], &B2[1] );
uint64_t S = B0[1] + B2[0];
OUT[0][2] = S >> 32;
OUT[0][3] = S & 0xFFFFFFFF;
uint64_t M = B2[1] + ( S < B2[0] );
OUT[0][1] = M & 0xFFFFFFFF;
OUT[0][0] = M >> 32;
}
{
B0[0] = _umul128( m[1], mul[1], &B0[1] );
B2[0] = _umul128( m[0], mul[1], &B2[1] );
uint64_t S = B0[1] + B2[0];
OUT[1][2] = S >> 32;
OUT[1][3] = S & 0xFFFFFFFF;
uint64_t M = B2[1] + ( S < B2[0] );
OUT[1][1] = M & 0xFFFFFFFF;
OUT[1][0] = M >> 32;
}
}
constexpr uint32_t N = 1 << 28;
int main()
{
uint32_t OUT[2][4];
uint32_t ABCD[4] = { 4294967295u, 4294967295u, 4294967295u, 4294967295u };
uint32_t EFGH[4] = { 4294967295u, 4294967295u, 4294967295u, 4294967295u };
multiply128x128_1( ABCD, EFGH, OUT );
uint64_t S_1 = 0u;
uint64_t S_2 = 0u;
uint64_t S_3 = 0u;
auto start_1 = std::chrono::high_resolution_clock::now();
for ( uint32_t i = 0; i < N; ++i )
{
EFGH[0] = i;
EFGH[1] = i;
EFGH[2] = i + 1;
EFGH[3] = i + 1;
ABCD[0] = i;
ABCD[1] = i;
ABCD[2] = i + 1;
ABCD[3] = i + 1;
multiply128x128( ABCD, EFGH, OUT );
S_1 += OUT[0][0] + OUT[0][1] + OUT[0][2] + OUT[0][3];
S_1 += OUT[1][0] + OUT[1][1] + OUT[1][2] + OUT[1][3];
}
auto stop_1 = std::chrono::high_resolution_clock::now();
std::cout << "Test A: " << std::chrono::duration_cast<std::chrono::milliseconds>( stop_1 - start_1 ).count() << '\n';
auto start_2 = std::chrono::high_resolution_clock::now();
for ( uint32_t i = 0; i < N; ++i )
{
EFGH[0] = i;
EFGH[1] = i;
EFGH[2] = i + 1;
EFGH[3] = i + 1;
ABCD[0] = i;
ABCD[1] = i;
ABCD[2] = i + 1;
ABCD[3] = i + 1;
mulShift( reinterpret_cast<const uint64_t*>( ABCD ), reinterpret_cast<const uint64_t*>( EFGH ), OUT );
S_2 += OUT[0][0] + OUT[0][1] + OUT[0][2] + OUT[0][3];
S_2 += OUT[1][0] + OUT[1][1] + OUT[1][2] + OUT[1][3];
}
auto stop_2 = std::chrono::high_resolution_clock::now();
std::cout << "Test B: " << std::chrono::duration_cast<std::chrono::milliseconds>( stop_2 - start_2 ).count() << '\n';
auto start_3 = std::chrono::high_resolution_clock::now();
for ( uint32_t i = 0; i < N; ++i )
{
EFGH[0] = i;
EFGH[1] = i;
EFGH[2] = i + 1;
EFGH[3] = i + 1;
ABCD[0] = i;
ABCD[1] = i;
ABCD[2] = i + 1;
ABCD[3] = i + 1;
multiply128x128_1( ABCD, EFGH, OUT );
S_3 += OUT[0][0] + OUT[0][1] + OUT[0][2] + OUT[0][3];
S_3 += OUT[1][0] + OUT[1][1] + OUT[1][2] + OUT[1][3];
}
auto stop_3 = std::chrono::high_resolution_clock::now();
std::cout << "Test C: " << std::chrono::duration_cast<std::chrono::milliseconds>( stop_3 - start_3 ).count() << '\n';
std::cout << S_1 << " " << S_2 << " " << S_3 << '\n';
}
为什么 _umul128 这么慢？也许我在上面的测试代码中犯了一些错误？
我的结果： 测试 A (simd)：4546 毫秒。 测试 B (_umul128)：6637 毫秒。 测试 C（标量）：2333 毫秒。
在 Windows 10、x64、MSVC 2019 上测试
_umul128 版本并没有那么慢，但你通过摆弄 32 位数组使 MSVC发出可怕的 asm.
优化正在击败您的基准；纯 C 版本并没有那么快。
尤其是简单的输入数据：
ABCD[0] = EFGH[0] = i;
ABCD[1] = EFGH[1] = i;
ABCD[2] = EFGH[2] = i + 1;
ABCD[3] = EFGH[3] = i + 1;
像这样初始化两个输入在内联纯 C 版本后创造了大量的优化机会。 i*i 4 次，i*(i+1) = i*i + i 8 次，(i+1)*(i+1) 4 次。 MSVC 并不愚蠢并且注意到了这一点。 这叫做 Common Subexpression Elimination (CSE)。
如果你想看看纯 C 到底有多慢，你需要想出一种更复杂的方法来伪造输入。也许提前生成然后遍历包含输入的内存？从循环计数器设置输入的成本几乎与乘法一样多。
MSVC 的 asm 输出证实大部分工作针对纯 C 版本进行了优化。 (Godbolt with MSVC 19.22 for x64)
...
$LL10@main:
lea r15, QWORD PTR [rax+1]
mov rcx, r15
mov r9, r15
imul rcx, rax # only 3, not 16, imul instructions.
imul rax, rax # (None appear later in this loop in the ... part)
imul r9, r15
mov edi, ecx
mov r14, rcx
mov r8d, eax
shr r14, 32 ; 00000020H
shr rax, 32 ; 00000020H
...
sub r13, 1
jne $LL10@main
MSVC 不善于优化内在函数 并且执行所有 4 个 mul m64 指令而不是注意到 ii * i1i1 执行了两次。
更重要的是，_umul128 循环受到存储转发停顿的伤害 因为它实际上 将数组存储到内存中使用 32 位存储，然后使用 64 位加载来提供 mul m64.
此外，处理 32 位块中的输出只会搬起石头砸自己的脚，引入额外的移位和 mov 操作。
这并不复杂，字面上只需要 3 条指令，mul r64 和 imul r64, r64 加上用于高半部分的 add，就足够了。 GCC/clang 轻松发出正确的东西，x86-64 System V 调用约定可以 return 寄存器中的 128 位 int。
关于神马：https://godbolt.org/z/DcZhSl
#include <stdint.h>
#ifdef __GNUC__
typedef unsigned __int128 u128;
u128 mul128x64( u128 a, uint64_t b) {
return a * b;
}
#endif
# clang -O3 for the x86-64 System V ABI (Linux)
mul128x64(unsigned __int128, unsigned long): #
mov rax, rdi
imul rsi, rdx
mul rdx
add rdx, rsi
ret
对于 MSVC，我们必须自己做，调用约定意味着结果 returned 在内存中。
#ifdef _MSC_VER
#include <intrin.h>
struct u128 { uint64_t u64[2]; };
u128 mul128x64( uint64_t a_lo, uint64_t a_hi, uint64_t b)
{
uint64_t lolo_high;
uint64_t lolo = _umul128( a_lo, b, &lolo_high );
uint64_t lohi = a_hi * b;
return {{lolo, lohi + lolo_high}};
}
#endif
# MSVC x64 -O2
u128 mul128x64(unsigned __int64,unsigned __int64,unsigned __int64) PROC
mov rax, r9
mul rdx
imul r8, r9
mov QWORD PTR [rcx], rax # store the retval into hidden pointer
mov rax, rcx
add r8, rdx
mov QWORD PTR [rcx+8], r8
ret 0
您的 __m128i 内在函数版本不太可能成功 。现代 x86（主流 Intel SnB 系列，AMD Ryzen）在 mul 和 imul 上具有 1/clock 吞吐量。 （除了 Ryzen，其中加宽 i/mul r64 的吞吐量为 2c，但对于 imul r64,r64 仍然是 1/clock。）
因此，如果您在 C 中实现并像这样编译成 asm，则 Sandybridge 系列上 64 x 128 位乘法的总吞吐量是每 2 个周期一次（在端口 1 上出现瓶颈）。
鉴于您需要超过 4 pmuludq 条指令来实现乘法运算，AVX1 是行不通的。 （Skylake 的 pmuludq 吞吐量为 0.5c。Sandybridge 的吞吐量为 1c，因此您需要在每次乘法（平均）2 pmuludq insns 内完成工作才能与标量竞争。而且这还没有考虑所有需要做的轮班/洗牌/添加工作。
可能值得考虑 Bulldozer 系列，其中 64 位标量乘法是 4c 吞吐量，但 pmuludq 是 1c。 (https://agner.org/optimize/) 每个周期产生 128 个产品位（两个 32x32 => 64 位产品）比每 4 个周期产生 128 个产品位要好，如果你可以在不消耗太多额外周期的情况下移动和添加它们。
同样，MSVC 不擅长通过内在函数进行持续传播或 CSE 优化，因此您的内在函数版本没有任何好处。
您的测试代码还使用来自标量整数循环变量的 _mm_set1_epi32( )，需要 vmovd 和 vpshufd 指令。
并且您为这些数组上的 lddqu 内在函数重新加载标量存储/向量，因此您再次遇到存储转发停顿。
SSE2 或 AVX1 的唯一希望是你的数据来自内存，而不是寄存器。或者，如果您可以将数据长期保存在矢量寄存器中，而不是不断地来回移动数据。特别是在 int <-> SIMD 具有高延迟的 Bulldozer 系列上。