Продукт Fast Dot для особого случая

Question

Учитывая вектор X размера L, где каждый скалярный элемент X взят из двоичного набора {0,1}, найти произведение точки z = точка (X, Y), если вектор Y размера L состоит целочисленных элементов. Я полагаю, должен существовать очень быстрый способ сделать это.

Допустим, у нас есть L=4; X[L]={1, 0, 0, 1}; Y[L]={-4, 2, 1, 0}, и мы должны найти z=X[0]*Y[0] + X[1]*Y[1] + X[2]*Y[2] + X[3]*Y[3] (что в этом случае даст нам -4).

Очевидно, что X можно представить с помощью двоичных цифр, например, целочисленный тип int32 для L = 32. Затем все, что нам нужно сделать, это найти точечное произведение этого целого числа с массивом из 32 целых чисел. У вас есть идеи или предложения, как это сделать очень быстро?

Answer 1

Это действительно потребует профилирования, но вы можете рассмотреть альтернативу:

int result=0;
int mask=1;
for ( int i = 0; i < L; i++ ){
    if ( X & mask ){
        result+=Y[i];
    }
    mask <<= 1;
}

Обычно операции сдвига битов и побитовые операции выполняются быстрее, чем умножение, однако оператор if может быть медленнее, чем умножение, хотя при прогнозировании ветвлений и больших значениях L я думаю, что это может быть быстрее. Вы действительно должны были бы профилировать это, однако, чтобы определить, привело ли это к любому ускорению.

Как было указано в комментариях ниже, развертывание цикла либо вручную, либо с помощью флага компилятора (например, "-funroll-loops" в GCC) также может ускорить это (исключая условие цикла).

Редактировать
В комментариях ниже был предложен следующий хороший твик:

int result=0;
for ( int i = 0; i < L; i++ ){
    if ( X & 1 ){
        result+=Y[i];
    }
    X >>= 1;
}

Answer 2

Попробуйте это:

int result=0;
for ( int i = 0; i < L; i++ ){
    result+=Y[i] & (~(((X>>i)&1)-1));
}

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

Answer 3

Является ли предложение изучить SSE2 полезным? У него уже есть операции с типом точечного произведения, плюс вы можете тривиально выполнить 4 (или, может быть, 8, я забыл размер регистра) простые итерации вашего наивного цикла параллельно. В SSE также есть несколько простых операций логического типа, поэтому он может выполнять сложения, а не умножения без каких-либо условных операций ... опять же, вам нужно посмотреть, какие операции доступны.

Answer 4

Поскольку размер явно не имеет значения , я думаю, что следующий код, вероятно, является наиболее эффективным кодом общего назначения:

int result = 0;
for (size_t i = 0; i < 32; ++i)
    result += Y[i] & -X[i];

Битовое кодирование X просто ничего не вносит в таблицу (даже если цикл потенциально может завершиться раньше, как правильно заметил @Mathieu) Но пропуск if внутри цикла делает .

Конечно, развертывание цикла может значительно ускорить это, как отмечали другие.

Answer 5

Как насчет комбинирования циклического сдвига с небольшой справочной таблицей?

    int result=0;

    for ( int x=X; x!=0; x>>=4 ){
        switch (x&15) {
            case 0: break;
            case 1: result+=Y[0]; break;
            case 2: result+=Y[1]; break;
            case 3: result+=Y[0]+Y[1]; break;
            case 4: result+=Y[2]; break;
            case 5: result+=Y[0]+Y[2]; break;
            case 6: result+=Y[1]+Y[2]; break;
            case 7: result+=Y[0]+Y[1]+Y[2]; break;
            case 8: result+=Y[3]; break;
            case 9: result+=Y[0]+Y[3]; break;
            case 10: result+=Y[1]+Y[3]; break;
            case 11: result+=Y[0]+Y[1]+Y[3]; break;
            case 12: result+=Y[2]+Y[3]; break;
            case 13: result+=Y[0]+Y[2]+Y[3]; break;
            case 14: result+=Y[1]+Y[2]+Y[3]; break;
            case 15: result+=Y[0]+Y[1]+Y[2]+Y[3]; break;
        }
        Y+=4;
    }

Производительность этого будет зависеть от того, насколько хорош компилятор в оптимизации оператора switch, но, по моему опыту, в настоящее время они довольно хороши в этом ....

Answer 6

Я видел несколько ответов с хитростью (чтобы избежать ветвления), но ни один из них не получил правильного цикла imho: /

Оптимизация @Goz ответ:

int result=0;
for (int i = 0, x = X; x > 0; ++i, x>>= 1 )
{
   result += Y[i] & -(int)(x & 1);
}

Преимущества:

• нет необходимости делать i операции сдвига битов каждый раз (X>>i)
• цикл останавливается раньше, если X содержит 0 в старших битах

Теперь мне интересно, работает ли он быстрее, тем более что преждевременная остановка цикла for может быть не такой легкой для развертывания цикла (по сравнению с константой времени компиляции).

Answer 7

Это решение идентично, но немного быстрее (по моим тестам), чем решение Майкла Аарона:

long Lev=1;
long Result=0
for (int i=0;i<L;i++) {
  if (X & Lev)
     Result+=Y[i];
  Lev*=2;
}

Я думал, что существует численный способ быстро установить следующий установленный бит в слове, который должен улучшить производительность, если ваши данные X очень скудны, но в настоящее время не могут найти указанную числовую формулировку в настоящее время.

Answer 8

Вполне вероятно, что время, потраченное на загрузку X и Y из основной памяти, будет доминировать. Если это так для архитектуры вашего процессора, алгоритм загружается быстрее при меньшей загрузке. Это означает, что хранение X в качестве битовой маски и расширение ее в кэш-память L1 ускорит алгоритм в целом.

Другой важный вопрос - будет ли ваш компилятор генерировать оптимальные нагрузки для Y. Это сильно зависит от процессора и компилятора. Но в целом полезно, если компилятор может точно видеть, какие значения нужны когда. Вы можете вручную развернуть цикл. Однако, если L является константой, оставьте это компилятору:

template<int I> inline void calcZ(int (&X)[L], int(&Y)[L], int &Z) {
  Z += X[I] * Y[I]; // Essentially free, as it operates in parallel with loads.
  calcZ<I-1>(X,Y,Z);
}
template< > inline void calcZ<0>(int (&X)[L], int(&Y)[L], int &Z) {
  Z += X[0] * Y[0];
}
inline int calcZ(int (&X)[L], int(&Y)[L]) {
    int Z = 0;
    calcZ<L-1>(X,Y,Z);
    return Z;
}

(Конрад Рудольф подверг сомнению это в комментарии, задаваясь вопросом о памяти использование . Это не реальное узкое место в современных компьютерных архитектурах, пропускная способность между памятью и процессором равна. Этот ответ почти не имеет значения, если Y как-то уже в кеше.)

Answer 9

Для small L вы можете использовать оператор switch вместо цикла. Например, если L = 8, вы можете иметь:

int dot8(unsigned int X, const int Y[])
{
    switch (X)
    {
       case 0: return 0;
       case 1: return Y[0];
       case 2: return Y[1];
       case 3: return Y[0]+Y[1];
       // ...
       case 255: return Y[0]+Y[1]+Y[2]+Y[3]+Y[4]+Y[5]+Y[6]+Y[7];
    }
    assert(0 && "X too big");
}   

И если L = 32, вы можете написать функцию dot32 (), которая вызывает dot8 () четыре раз, если возможно, встроенная. (Если ваш компилятор отказывается встроить dot8 (), вы можете переписать dot8 () как макрос для принудительного встраивания.) Добавлено :

int dot32(unsigned int X, const int Y[])
{
    return dot8(X >> 0  & 255, Y + 0)  +
           dot8(X >> 8  & 255, Y + 8)  +
           dot8(X >> 16 & 255, Y + 16) +
           dot8(X >> 24 & 255, Y + 24);
}

Это решение, как указывает Микера, может иметь стоимость кэша команд; в этом случае может помочь использование функции dot4 ().

Дальнейшее обновление : это можно комбинировать с решением Микеры:

static int dot4(unsigned int X, const int Y[])
{
    switch (X)
    {
        case 0: return 0;
        case 1: return Y[0];
        case 2: return Y[1];
        case 3: return Y[0]+Y[1];
        //...
        case 15: return Y[0]+Y[1]+Y[2]+Y[3];
    }
}

Глядя на полученный ассемблерный код с опциями -S -O3 с gcc 4.3.4 на CYGWIN, я немного удивлен, увидев, что это автоматически встроено в dot32 (), с восемь таблицы переходов по 16 записей.

Но добавление __attribute __ ((__ noinline__)), кажется, производит более привлекательный ассемблер.

Другим вариантом является использование сквозных операторов в операторе switch, но gcc добавляет инструкции jmp и выглядит не так быстро.

Редактировать - Совершенно новый ответ: После размышления о 100-процентном штрафе, упомянутом Антсом Аасмой, и другими ответами, вышеупомянутое, вероятно, не оптимально. Вместо этого вы можете вручную развернуть цикл как в:

int dot(unsigned int X, const int Y[])
{
    return (Y[0] & -!!(X & 1<<0)) +
           (Y[1] & -!!(X & 1<<1)) +
           (Y[2] & -!!(X & 1<<2)) +
           (Y[3] & -!!(X & 1<<3)) +
           //...
           (Y[31] & -!!(X & 1<<31));
}

Это, на моей машине, генерирует 32 x 5 = 160 быстрых инструкций. Умный компилятор мог предположительно развернуть другие предложенные ответы, чтобы дать тот же результат.

Но я все еще проверяю дважды.

Answer 10

result = 0;
for(int i = 0; i < L ; i++)
    if(X[i]!=0)
      result += Y[i];