x86 - 2 个 64 位整数的 SSE 乘法

Question

如何将两个 64 位整数乘以另外 2 个 64 位整数？我没有找到任何可以做到这一点的指令。

Answer 1

迟到的答案，但这是巴拉巴斯发布的更好版本。

如果您曾经使用过 GCC 或 Clang 的向量扩展，这就是他们使用的例程。

这使用与长乘法和网格乘法相同的方法。

    65
  * 73
  ----
    15 //   (5 * 3)
   180 //   (6 * 3) * 10
   350 //   (5 * 7) * 10
+ 4200 // + (6 * 7) * 100
------
  4745

但是，它不是使用 10 的每个单位，而是使用 32 位的每个单位，并且它省略了最后一个乘法，因为它总是会移动到第 64 位之后，就像你不会乘以 6*7 一样截断大于 99 的值。

#include <emmintrin.h>

/*
 * Grid/long multiply two 64-bit SSE lanes.
 * Works for both signed and unsigned.
 *   ----------------.--------------.----------------.
 *  |                |   b >> 32    | a & 0xFFFFFFFF |
 *  |----------------|--------------|----------------|  
 *  | d >> 32        |   b*d << 64  |    a*d << 32   |
 *  |----------------|--------------|----------------|
 *  | c & 0xFFFFFFFF |   b*c << 32  |       a*c      |
 *  '----------------'--------------'----------------'
 *  Add all of them together to get the product.
 *
 *  Because we truncate the value to 64 bits, b*d << 64 will be zero,
 *  so we can leave it out.
 *
 *  We also can add a*d and b*c first and then shift because of the
 *  distributive property: (a << 32) + (b << 32) == (a + b) << 32.
 */

__m128i Multiply64Bit(__m128i ab, __m128i cd)
{
    /* ac = (ab & 0xFFFFFFFF) * (cd & 0xFFFFFFFF); */
    __m128i ac = _mm_mul_epu32(ab, cd);

    /* b = ab >> 32; */
    __m128i b = _mm_srli_epi64(ab, 32);

    /* bc = b * (cd & 0xFFFFFFFF); */
    __m128i bc = _mm_mul_epu32(b, cd);

    /* d = cd >> 32; */
    __m128i d = _mm_srli_epi64(cd, 32);

    /* ad = (ab & 0xFFFFFFFF) * d; */
    __m128i ad = _mm_mul_epu32(ab, d);

    /* high = bc + ad; */
    __m128i high = _mm_add_epi64(bc, ad);

    /* high <<= 32; */
    high = _mm_slli_epi64(high, 32);

    /* return ac + high; */
    return _mm_add_epi64(high, ac);
}

Compiler Explorer注意：下面还包含了 GCC 矢量扩展版本以供比较。

Answer 2

我知道这是一个老问题，但我实际上正在寻找这个。由于仍然没有关于它的指令，我实现了 64 位乘法自己与 Paul R 提到的 pmuldq。这就是我想出的：

// requires g++ -msse4.1 ...

#include <emmintrin.h>
#include <smmintrin.h>

__m128i Multiply64Bit(__m128i a, __m128i b)
{
    auto ax0_ax1_ay0_ay1 = a;
    auto bx0_bx1_by0_by1 = b;

    // i means ignored

    auto ax1_i_ay1_i = _mm_shuffle_epi32(ax0_ax1_ay0_ay1, _MM_SHUFFLE(3, 3, 1, 1));
    auto bx1_i_by1_i = _mm_shuffle_epi32(bx0_bx1_by0_by1, _MM_SHUFFLE(3, 3, 1, 1));

    auto ax0bx0_ay0by0 = _mm_mul_epi32(ax0_ax1_ay0_ay1, bx0_bx1_by0_by1);
    auto ax0bx1_ay0by1 = _mm_mul_epi32(ax0_ax1_ay0_ay1, bx1_i_by1_i);
    auto ax1bx0_ay1by0 = _mm_mul_epi32(ax1_i_ay1_i, bx0_bx1_by0_by1);

    auto ax0bx1_ay0by1_32 = _mm_slli_epi64(ax0bx1_ay0by1, 32);
    auto ax1bx0_ay1by0_32 = _mm_slli_epi64(ax1bx0_ay1by0, 32);

    return _mm_add_epi64(ax0bx0_ay0by0, _mm_add_epi64(ax0bx1_ay0by1_32, ax1bx0_ay1by0_32));
}

SSE Multiply64Bit的 Godbolt 。

Answer 3

您需要使用 32 位乘法运算来实现自己的 64 位乘法例程。不过，它可能不会比仅使用标量代码更有效，特别是因为要对向量进行大量改组才能获得所有必需的操作。