在 64 位 x 64 位乘法中使用 Karatsuba 算法真的有效吗？答案

【问题标题】：Is it really efficient to use Karatsuba algorithm in 64-bit x 64-bit multiplication?在 64 位 x 64 位乘法中使用 Karatsuba 算法真的有效吗？
【发布时间】：2015-09-13 04:16:32
【问题描述】：

我在 AVX2 上工作，需要计算 64 位 x64 位 -> 128 位加宽乘法，并以最快的方式获得 64 位高位部分。既然AVX2没有这样的指令，那我用Karatsuba算法提高效率和速度是否合理？

【问题讨论】：

这在很大程度上取决于架构。 25 年前，我确实在 32 位 SPARC 处理器上使用 Karatsuba 进行 64x64->128 位乘法，以提高性能。我还没有看过 AVX2，我的机器中缺少 Haswell 级 CPU。您是否搜索了文献（或一般的互联网）以查看其他人发现了什么？使用标准方法，您的基于 AVX2 的 64x64->128 位乘法有多快？
你真的需要 64bx64b 到 128b 吗？或者你可以使用 56bx56b 到 106b 吗？
其实我需要 64bx64b -> 64b 高位，这样codepaste.net/29m5qm `
另外我认为你应该在你的问题中说明你想要 64bx64b 到高 64b。
你的问题会更有趣，我想如果你问是否有一种使用 AVX2 的快速方法来做high(64bx64b)。您的问题可以被视为xy problem，因为您确实想找到一种更快的方法来处理high(64bx64b)。

标签： c++ performance parallel-processing simd avx2

【解决方案1】：

如果不尝试就很难判断，但仅使用 AMD64 MUL 指令可能会更快，它支持 64x64=128，吞吐量与大多数 AVX2 指令相同（但未矢量化）。缺点是如果操作数在 YMM 寄存器中，则需要加载到常规寄存器。这会给单个 64x64=128 类似 LOAD + MUL + STORE 的东西。

如果您可以在 AVX2 中矢量化 Karatsuba，请尝试 AVX2 和 MUL，看看哪个更快。如果你不能矢量化，单个MUL 可能会更快。如果您可以将负载移除并存储到常规寄存器中，那么单个MUL 肯定会更快。

MUL 和 AVX2 指令都可以在内存中具有相同吞吐量的操作数，这可能有助于为 MUL 消除一个负载。

【讨论】：

在零件上花费的时间最多：__uint128_t product = ((__uint128_t)a)*((__uint128_t)b);那么你能对这个答案说些什么呢？ stackoverflow.com/a/24575626/1979163
@user1979163 您是否启用了所有优化？使用 GCC，整个函数只转换为 4 条汇编指令（mov,mov,mul,mov）。这就是禁用内联的情况。并且 4 条指令在 CPU 内部应该具有良好的并行吞吐量。我不确定你能不能做得更好。
是的，我启用了所有优化并编译为 g++ -mavx2 -O3 -ftree-vectorize -ftree-vectorizer-verbose=2 mulhiExample.cpp -std=c++11 -fabi-version=0 。你能提供你编译的代码和你使用的标志吗？
@user1979163 我使用了您提供的代码，删除了静态内联。我用你的选项得到了相同的程序集，所以你可能有相同的结果。此功能是否仅在同一个 .cpp 文件中使用？如果它在其他地方使用，它不会被内联，即使是inline。在这种情况下，它可能会花费很多。要启用内联，请将函数定义放在 .h 中，或使用 -flto 进行编译。如果你只在同一个 .cpp 文件中使用这个函数，你可能不走运。
@user1979163 注意，如果使用-flto，还需要将-O3传递给链接器命令，否则将无用。

【解决方案2】：

It's highly unlikely that AVX2 will beat the mulx instruction 在一条指令中执行 64bx64b 到 128b。我知道有一个例外large multiplications using floating point FFT。

但是，如果您不需要 64bx64b 到 128b，您可以考虑 53bx53b 到 106b 使用double-double arithmetic。

四个53位数字a和b相乘得到四个106位数字只需要两条指令：

__m256 p = _mm256_mul_pd(a,b);
__m256 e = _mm256_fmsub_pd(a,b,p);

与使用mulx 的一条指令中的一个 128 位数字相比，这在两条指令中提供了四个 106 位数字。

【讨论】：

我很确定所有 64+64 位都进入了高位（64bx64b）。只是一个简单的例子（32bx32b）：A = 0xE567 89AB，B = C123 4567，高（AxB）= 0xAD12 AA22。现在，如果我从 A 和 B 中仅删除 1 位：A = 0xE567 89AA，B = 0xC123 4566，high(AxB) = 0xAD12 AA21。 2x32b 中的 2x1b 进入高位（32bx32b）。还是我的推理有问题？
@Zboson：很遗憾，您无法修复它（或者至少，性能不佳）。您确实需要所有输入位来计算正确的乘高（这就是我在提到“操作数的舍入”时所指的）。
@Zboson - 我没有遵循答案中给出的 sn-p。 fmsub(x,y,z) 是 x * y - z，对吧？据我所知，sn-p 是p = a * b; e = a * b - (p) = a * b - (a * b) = ~0，对吧？
@Zboson 有办法让它工作。但从我的测试来看，还不足以克服MULX/ADCX/ADOX。所以我很期待AVX512。但它将是短暂的，因为 AVX512-IFMA 将取代它。
诀窍是将输入限制为 e 的 fma 之前舍入 p = a*b。但要完成这项工作，您需要正确缩放输入，以便 a*b 在中间被小数位分割。顺序是：p = a*b; p = round(p); e = fma(-a,b,p) 但是你的输出p 和e 的缩放比例也不同。最终，设置和校正缩放的开销使它比通常的标量方法慢。有一些技巧可以消除一些缩放，但很难消除足够多的缩放以使该方法可行。

【解决方案3】：

没有。在现代架构中，Karatsuba 击败教科书乘法的交叉点通常介于 8 到 24 个机器字之间（例如，在 x86_64 上介于 512 到 1536 位之间）。对于固定大小，阈值位于该范围的较小端，新的 ADCX/ADOX 指令可能会在一定程度上进一步提高标量代码，但 64x64 仍然太小，无法从 Karatsuba 中受益。

【讨论】：

这回答了 OPs 问题，但我认为如果 OP 询问是否有办法使用 AVX2 击败 mulx 为 high(64bx64b) 会更有趣（这就是我选择解释 OPs 问题的方式; 即作为 xy 问题）。我认为在浮点域中使用双双算法可能是可能的（假设 OP 可以在大多数算法中保持在浮点域中）。
如果提问者真的在做整数运算，并且乘法是所讨论代码的重要部分，则不可能击败 MULX。使用 FMA 很可爱，但如果任务本质上是整数，则会遇到两个痛点：(a) 操作数和高乘积的舍入，以及 (b) 如果输入并不总是全宽，则乘积最终会有不同的指数，这使得后续操作非常痛苦。
你可能是对的。我还没有尝试实现这一点。这是我一直在考虑的一个想法。从某种意义上说，舍入可能是一个前期成本，这就是我在大多数算法中保持浮点域的意思。不同的指数是个问题。也许更大的问题是是否也需要添加。双双加法很慢，回到整数域进行加法是另一个开销，AVX2 比较和加法无论如何都可以用 ADCX 收支平衡。
根据我的经验，Karatsuba 的“有用”范围在现代处理器上正在迅速缩小。像 MULX 和 ADX 这样的东西正在使基本情况算法越来越快，从而推高了基本情况 / karatsuba 阈值。另一方面，浮点 FFT 正在充分利用 SIMD 的优势，它正在蚕食 Karatsuba/FFT 阈值并将其推低。
@Zboson：在 32b 过程中，这是一个更更具吸引力的选择，因为您拥有的最宽标量乘数是 32x32 -> 64。但真正的解决方案是只需在 64b 中编译 =)。