Головоломка по оптимизации SIMD

Это настолько близко, насколько я мог приблизиться к векторизации этого кода. Я действительно не ожидаю, что это будет быстрее. Я просто пробовал свои силы в написании SIMD-кода.

#include <stdint.h>

int64_t fun(int64_t N, int size, const int* p)
{
    int64_t sum = 0;
    int i;
    for(i=1; i<size; i++) {
        sum += (N/i)*p[i];
    }
    return sum;
}

typedef int64_t v2sl __attribute__ ((vector_size (2*sizeof(int64_t))));

int64_t fun_simd(int64_t N, int size, const int* p)
{
    int64_t sum = 0;
    int i;
    v2sl v_2 = { 2, 2 };
    v2sl v_N = { N, N };
    v2sl v_i = { 1, 2 };
    union { v2sl v; int64_t a[2]; } v_sum;

    v_sum.a[0] = 0;
    v_sum.a[1] = 0;
    for(i=1; i<size-1; i+=2) {
        v2sl v_p = { p[i], p[i+1] };
        v_sum.v += (v_N / v_i) * v_p;
        v_i += v_2;
    }
    sum = v_sum.a[0] + v_sum.a[1];
    for(; i<size; i++) {
        sum += (N/i)*p[i];
    }
    return sum;
}

typedef double v2df __attribute__ ((vector_size (2*sizeof(double))));

int64_t fun_simd_double(int64_t N, int size, const int* p)
{
    int64_t sum = 0;
    int i;
    v2df v_2 = { 2, 2 };
    v2df v_N = { N, N };
    v2df v_i = { 1, 2 };
    union { v2df v; double a[2]; } v_sum;

    v_sum.a[0] = 0;
    v_sum.a[1] = 0;
    for(i=1; i<size-1; i+=2) {
        v2df v_p = { p[i], p[i+1] };
        v_sum.v += (v_N / v_i) * v_p;
        v_i += v_2;
    }
    sum = v_sum.a[0] + v_sum.a[1];
    for(; i<size; i++) {
        sum += (N/i)*p[i];
    }
    return sum;
}

#include <stdio.h>

static const int test_array[] = {
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0
};
#define test_array_len (sizeof(test_array)/sizeof(int))

#define big_N (1024 * 1024 * 1024)

int main(int argc, char *argv[]) {
    int64_t res1;
    int64_t res2;
    int64_t res3;
    v2sl a = { 123, 456 };
    v2sl b = { 100, 200 };
    union { v2sl v; int64_t a[2]; } tmp;

    a = a + b;
    tmp.v = a;
    printf("a = { %ld, %ld }\n", tmp.a[0], tmp.a[1]);

    printf("test_array size = %zd\n", test_array_len);

    res1 = fun(big_N, test_array_len, test_array);
    printf("fun() = %ld\n", res1);

    res2 = fun_simd(big_N, test_array_len, test_array);
    printf("fun_simd() = %ld\n", res2);

    res3 = fun_simd_double(big_N, test_array_len, test_array);
    printf("fun_simd_double() = %ld\n", res3);

    return 0;
}

Производная от 1/x равна -1/x^2, что означает, что чем больше x, тем больше N/x==N/(x + 1).

Для известного значения N/x (назовем это значение r) мы можем определить следующее значение x (назовем это значение x' так, чтобы N/x'<r:

x'= N/(r - 1)

А так как мы имеем дело с целыми числами:

x'= ceiling(N/(r - 1))

Таким образом, цикл становится примерно таким:

int64_t sum = 0;   
int i=1; 
int r= N;
while (i<size)
{
    int s= (N + r - 1 - 1)/(r - 1);

    while (i<s && i<size)
    {
        sum += (r)*p[i];
        ++i;
    }

    r= N/s;
}
return sum;   

Для достаточно большого N у вас будет много прогонов одинаковых значений для N/i. Конечно, вы попадете в деление на ноль, если не будете осторожны.

Я предлагаю вам сделать это с операциями SIMD с плавающей запятой - с одинарной или двойной точностью в зависимости от ваших требований к точности. Преобразование из int в float или double выполняется относительно быстро с использованием SSE.

Стоимость сосредоточена на вычислении подразделений. В SSE2 нет кода операции для целочисленного деления, поэтому вам придется реализовать алгоритм деления самостоятельно, побитно. Я не думаю, что это стоило бы усилий: SSE2 позволяет вам выполнять два экземпляра параллельно (вы используете 64-битные числа, а регистры SSE2 128-битные), но я считаю вероятным, что ручной алгоритм деления будет по крайней мере в два раза медленнее, чем код операции ЦП idiv.

(Кстати, вы компилируете в 32-битном или 64-битном режиме? В последнем будет удобнее работать с 64-битными целыми числами.)

Сокращение общего числа дивизий выглядит более перспективным путем. Можно заметить, что для положительных целых чисел x и y, тогда floor(x/(2y)) = floor(floor(x/y)/2)< /эм>. В терминологии C, как только вы вычислили N/i (усеченное деление), вам просто нужно сдвинуть его на один бит вправо, чтобы получить N/(2*i). При правильном использовании это делает половину ваших делений почти бесплатными (что «правильно» также включает доступ к миллиардам значений p[i] таким образом, который не наносит ущерба кешам, поэтому это не кажется очень простым).