字符/数字的乘法可以更高效吗？

我有以下代码,其中计算总和,基于一个非常大的系列.

该系列char *a是一个char数组,仅包含数字(0..9).

我想问一下是否有可能让代码更快.它目前是分布式计算应用程序的瓶颈.

一个小的复制代码.不是实际的代码,而是更简化.

int top = 999999999;

char *a;
a = (char*) calloc(top+1, sizeof(char));

// ... fill a with initial values ...

for (int i=0; i<10; ++i) {
    unsigned long long int sum = 0;

    for (m = 1, k = top; m < k; ++m, --k) {
        // Here is the bottle neck!!
        sum += a[m]*a[k];
    }

    printf("%d\n", sum);

    // ... Add something at the end of a, and increase top ...
}

我已经尝试过以下内容:

使用-O3(gcc编译器)优化代码.编译器行现在是:

gcc -c -Wall -fopenmp -Wno-unused-function -O3 -std=c99 -g0 -march=native -pipe -D_FILE_OFFSET_BITS=64 -m64 -fwhole-program -fprefetch-loop-arrays -funsafe-loop-optimizations -Wunsafe-loop-optimizations -fselective-scheduling -fselective-scheduling2 -fsel-sched-pipelining -fsel-sched-pipelining-outer-loops -fgcse-sm -fgcse-lm -fgcse-las -fmodulo-sched -fgcse-after-reload -fsee -DLIBDIVIDE_USE_SSE2 -DLIBDIVIDE_USE_SSE4_1 xxx.c -o xxx.o

使用GNU openMP将for循环拆分为多个核心

unsigned long long int halfway = (top>>1) + 1; // = top/2 + 1
// digits is defined as top+1

#pragma omp parallel // firstprivate/*shared*/(a, digits, halfway)
for (unsigned long long int m = 1; m < halfway; ++m) {
    sum += a[m] * a[digits-m];
}

结果:更快,更快,但需要更多核心,我仍然希望加快速度.

铸造a[m]到unsigned long long int乘法之前

sum += (unsigned long long int)a[m] * a[k];

结果:性能提升很小.

使用乘法查找表,因为数组查找比实际乘法更快.

sum += multiply_lookup[a[m]][a[k]]; // a[m]*a[k];

结果:性能提升很小.

我试图找到一种减少操作的数学解决方案,但似乎没有任何东西可以通过数学方式进行优化.

我有以下想法进行优化:

我已经读过浮点数(asm fmul)的乘法比整数乘法(asm mul)快得多.只是更改int为float没有帮助 - 但我认为如果使用MMX或SSE指令集完成工作,或者如果工作由FPU完成,则代码可能会变得更加高效.虽然我有一些汇编知识,但我对这些主题一无所知.

但是,如果您有其他想法如何优化它,我很高兴听到它们.

更新一些其他信息:

该系列在每个循环后增长1个元素.

随着系列的增长,top增加.

当top达到数组限制时,a将使用增加100000字节realloc().

平台:Debian Linux Jessie x64,在Intel(R)Xeon(R)CPU X3440 @ 2.53GHz上

额外的题外话题:你知道这个总和的数学名称,该系列的元素对从外到内成倍增加吗？

您可以使用鲜为人知的PMADDUBSW(Multiply和Add Packed Signed和Unsigned Bytes).签名/未签名的业务在这里无关紧要,无论如何,一切都在[0 .. 9]的区间内.添加是饱和的,但这并不重要,因为9*9只有81.内在函数是_mm_maddubs_epi16.因为k索引会下降,所以你必须对它进行字节反转,这可以用PSHUFB(_mm_shuffle_epi8)来完成.当索引在中间"相遇"时会发生令人讨厌的事情,你可以逐一完成这一部分.

这是一个尝试,只有轻微测试:

__m128i sum = _mm_setzero_si128();
int m, k;
for (m = 1, k = top - 15; m + 15 < k; m += 16, k -= 16) {
   __m128i am = _mm_loadu_si128((__m128i*)(a + m));
   __m128i ak = _mm_loadu_si128((__m128i*)(a + k));
   ak = _mm_shuffle_epi8(ak, _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ,15));
   sum = _mm_add_epi16(sum, _mm_maddubs_epi16(am, ak));
}
// could use phaddw, but I do this the long way to avoid overflow slightly longer
sum = _mm_add_epi32(_mm_unpacklo_epi16(sum, _mm_setzero_si128()),
                    _mm_unpackhi_epi16(sum, _mm_setzero_si128()));
sum = _mm_hadd_epi32(sum, sum);
sum = _mm_hadd_epi32(sum, sum);
int s = _mm_cvtsi128_si32(sum);
// this is for the "tail"
k += 15;
for (; m < k; ++m, --k)
    s += a[m] * a[k];

我也忽略溢出.您可以为(2 ¹⁶ -1)/(2*81)= 404迭代执行此操作,但仍然绝对没有溢出.如果您需要更多,请定期将其添加到32位结果中.

在快速基准测试中,这个速度大约是简单方法的7倍(在4770K上使用2KB的随机数据进行测试,每次运行最多可以进行100次运行).

使用其他答案所建议的指针可以进一步提高它,大约是简单方法的9倍.有了索引,有一些奇怪的符号扩展正在进行.

int foobar(char* a, int top)
{
    __m128i sum = _mm_setzero_si128();

    char *m, *k;
    for (m = a + 1, k = a + top - 15; m + 15 < k; m += 16, k -= 16) {
       __m128i am = _mm_loadu_si128((__m128i*)(m));
       __m128i ak = _mm_loadu_si128((__m128i*)(k));
       ak = _mm_shuffle_epi8(ak, _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15));
       sum = _mm_add_epi16(sum, _mm_maddubs_epi16(am, ak));
    }

    sum = _mm_add_epi32(_mm_unpacklo_epi16(sum, _mm_setzero_si128()),
                        _mm_unpackhi_epi16(sum, _mm_setzero_si128()));
    sum = _mm_hadd_epi32(sum, sum);
    sum = _mm_hadd_epi32(sum, sum);
    int s = _mm_cvtsi128_si32(sum);

    k += 15;
    for (; m < k; ++m, --k)
        s += *m * *k;

    return s;
}

除了额外的逻辑之外,分成几部分,仍然是原始速度的9倍:

int foobar(char* a, int top)
{
    int s = 0;
    char *m, *k;
    for (m = a + 1, k = a + top - 15; m + 15 < k;) {
        __m128i sum = _mm_setzero_si128();
        for (int i = 0; i < 404 && m + 15 < k; m += 16, k -= 16, ++i) {
           __m128i am = _mm_loadu_si128((__m128i*)(m));
           __m128i ak = _mm_loadu_si128((__m128i*)(k));
           ak = _mm_shuffle_epi8(ak, _mm_set_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ,15));
           sum = _mm_add_epi16(sum, _mm_maddubs_epi16(am, ak));
        }
        sum = _mm_add_epi32(_mm_unpacklo_epi16(sum, _mm_setzero_si128()),
                            _mm_unpackhi_epi16(sum, _mm_setzero_si128()));
        sum = _mm_hadd_epi32(sum, sum);
        sum = _mm_hadd_epi32(sum, sum);
        s += _mm_cvtsi128_si32(sum);
    }

    k += 15;
    for (; m < k; ++m, --k)
        s += *m * *k;

    return s;
}