SSE 等矢量编程 (SIMD)

Question

我对 SSE 编程完全陌生，但我拥有 Intel Core i7 处理器。

基本上，我想一次取 4 个 32 位无符号整数并将它们全部立方（提高到 3 的幂）。据我了解，SSE 及其后继产品的 SIMD 功能使这成为可能，但我到底该怎么做呢？最好在 C 中，但如果需要，我可以管理程序集。

编辑以明确我的最终目标：

然后，我想把所有的立方体加在一起得到一个数字。

背景：我只是想使用 SSE 来优化确定一个数字是否是 Armstrong 数字（一个三位数字，其每个数字的立方和与数字本身相同）。一个例子是 153。除了蛮力之外似乎没有办法做到这一点。这些是自恋数字的子集，其所有数字的总和到十进制数字的长度的幂等于数字本身。希望我最终将其扩展为更灵活，开始我只是在做阿姆斯壮的数字。正如您可能想象的那样，这出现在另一个网站上，我们中的一些人正在努力优化它。通过您的想法和我自己的研究，我想出了以下代码：

#include <stdio.h>
#include <smmintrin.h>  // SSE 4.1

__m128i vcube(const __m128i v)
{
    return _mm_mullo_epi32(v, _mm_mullo_epi32(v, v));
}


int main(int argc, const char * argv[]) {
    for (unsigned int i = 1; i <= 500; i++) {
        unsigned int firstDigit = i / 100;
        unsigned int secondDigit = (i - firstDigit * 100) / 10;
        unsigned int thirdDigit = (i - firstDigit * 100 - secondDigit * 10);

        __m128i v = _mm_setr_epi32(0, firstDigit, secondDigit, thirdDigit);
        __m128 v3 = (__m128) vcube(v);

        v3 = _mm_hadd_ps(v3, v3);
        v3 = _mm_hadd_ps(v3, v3);

        if (_mm_extract_epi32((__m128i) v3, 0) == i)    
            printf ("%03d is an Armstrong number\n", i);
        }
    return 0;
}

注意：我必须进行一些类型强制才能使其在某些系统（Solaris，至少是某些 Linux）中编译。

所以这行得通，但也许可以简化。抱歉，我没有发布整个任务，但我试图将其分解为多个步骤，并且我想确保每个数字都正确立方。

（结束编辑）

谢谢！

编辑：我想我应该添加我正在运行 Mac OS X Sierra。

再次编辑：

所以，假设我制作了所有这些无符号短裤而不是无符号整数并添加更多数字，当短裤可能无法容纳所有数字的总和时，我如何将它们加在一起？如果您知道我的意思，有没有办法将它们添加并存储在更大变量的向量中，或者像 UInt64 这样的普通更大的数字？

抱歉所有问题，但就像我说的，我在矢量处理方面完全是新手，尽管我从我的第一台 Mac G4 开始就可以使用它。

Answer 1

如果您的输入值在 0..1625 范围内（以便结果适合 32 位），那么您可以使用 _mm_mullo_epi32:

__m128i vcube(const __m128i v)
{
    return _mm_mullo_epi32(v, _mm_mullo_epi32(v, v));
}

---

演示：

#include <stdio.h>
#include <smmintrin.h>  // SSE 4.1

__m128i vcube(const __m128i v)
{
    return _mm_mullo_epi32(v, _mm_mullo_epi32(v, v));
}

int main()
{
    __m128i v = _mm_setr_epi32(0, 1, 1000, 1625);
    __m128i v3 = vcube(v);

    printf("%vlu => %vlu\n", v, v3);

    return 0;
}

编译测试：

$ gcc -Wall -Wno-format-invalid-specifier -Wno-format-extra-args -msse4 vcube.c && ./a.out
0 1 1000 1625 => 0 1 1000000000 4291015625

Answer 2

对于x<=2642245，您可以使用下面的foo_SSE 函数使用SSE4.1 执行x*x*x。这将两个 32 位无符号整数作为输入打包到 SSE 寄存器的高 64 位和低 64 位，并输出两个 64 位整数。

#include <stdio.h>
#include <x86intrin.h>
#include <inttypes.h>

__m128i foo_SSE(__m128i x) {
  __m128i mask = _mm_set_epi32(-1, 0, -1, 0);
  __m128i x2 =_mm_shuffle_epi32(x, 0x80);
  __m128i t0 = _mm_mul_epu32(x,x);
  __m128i t1 = _mm_mul_epu32(t0,x);
  __m128i t2 = _mm_mullo_epi32(t0,x2);
  __m128i t3 = _mm_and_si128(t2, mask);
  __m128i t4 = _mm_add_epi32(t3, t1);
  return t4;
}

int main(void) {
  uint64_t k1 = 100000;
  uint64_t k2 = 2642245;                                                                                                                                                             
  __m128i x = _mm_setr_epi32(k1, 0, k2, 0);
  uint64_t t[2];
  _mm_store_si128((__m128i*)t, foo_SSE(x));
  printf("%20" PRIu64 " ",  t[0]);
  printf("%20" PRIu64 "\n", t[1]);
  printf("%20" PRIu64 " ",  k1*k1*k1);
  printf("%20" PRIu64 "\n", k2*k2*k2);    
}

这可能会有所改进。我有点不习惯。

Answer 3

要快速了解 3 个主要阶段（加载、操作、存储），请参阅以下 sn-p。对于整数e0 和e1：

#include "emmintrin.h"
__m128i result __attribute__((aligned(16)));
__m128i x = _mm_setr_epi32(0, e1, 0, e0);
__m128i cube = _mm_mul_epu32(x, _mm_mul_epu32(x, x));
_mm_store_si128(&result, cube);

_mm_mul_epu32 取两个 _m128i 寄存器的 32 位的偶数倍，将它们相乘并将结果作为 64 位的 2 元组放入结果寄存器中。

要让它们离开那里，要么通过强制转换访问它们，要么使用编译器对 __m128i 的方便定义，例如国际商会：

printf("%llu %llu\n", result.m128i_i64[0], result.m128i_i64[1]); /* msc style */

注意：我将 Intel Intrinsics guide 用于 SSE 原语。

为了清楚代码的实际作用而进行了编辑。