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

Question

Я хочу оптимизировать следующую функцию с использованием SIMD (SSE2 и т. Д.):

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

Это похоже на чрезвычайно векторизуемую задачу, за исключением того, что необходимых инструкций просто нет ...

Можно предположить, что N очень велико (от 10 ^ 12 до 10 ^ 18) и имеет размер ~ sqrt (N).Мы также можем предположить, что p может принимать только значения -1, 0 и 1;таким образом, нам не нужно реальное умножение, (N / i) * p [i] может быть выполнено с четырьмя инструкциями (pcmpgt, pxor, psub, pand), если бы мы могли просто как-то вычислить N / i.

Answer 1

Производная от 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. Конечно, если вы не будете осторожны, вы достигнете деления на ноль.

Answer 2

Это как можно ближе к векторизации этого кода. Я не ожидаю, что это будет быстрее. Я просто пытался писать 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;
}

Answer 3

Затраты сосредоточены на вычислении делений. В 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] таким образом, который не наносит ущерб кэшам, поэтому это не кажется очень простым).

Answer 4

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