Быстрое скалярное произведение для совершенно особого случая

Представляет X , используя связанный список мест, где x [i] = 1 . {{1} } Чтобы найти нужную сумму, вам потребуется O (N) операций, где N - размер вашего списка.

Вы можете сохранить свой битовый вектор как последовательность of int, где каждый int упаковывает пару коэффициентов в виде битов. Тогда покомпонентное умножение эквивалентно битовому и. При этом вам просто нужно подсчитать количество установленных битов, что можно сделать следующим образом:

inline int count(uint32_t x) {
    // see link
}

int dot(uint32_t a, uint32_t b) {
    return count(a & b);
}

Чтобы немного взломать подсчет установленных битов, см. http://graphics.stanford.edu/~seander/bithacks. html # CountBitsSetParallel

Edit: Извините, я только что понял, что только один из векторов содержит элементы {0,1}, а другой - нет. Этот ответ применим только к случаю, когда оба вектора ограничены коэффициентами из набора {0,1}.

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

Попробуйте следующее:

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

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

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

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

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

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

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

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

    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, но, по моему опыту, в настоящее время они довольно хороши в этом ....

Для 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 ().

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

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 быстрых инструкций. Умный компилятор, вероятно, мог бы развернуть другие предложенные ответы, чтобы дать тот же результат.

Но я все еще перепроверяю.

Вполне вероятно, что время, затрачиваемое на загрузку 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;
}

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

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

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

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

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

На этот вопрос, вероятно, нет общего ответа. Вам нужно профилировать свой код под все различные цели. Производительность будет зависеть от оптимизаций компилятора, таких как разворачивание циклов и SIMD-инструкции, которые доступны на большинстве современных процессоров (x86, PPC, ARM имеют свои собственные реализации).

Ну, вы хотите, чтобы все биты прошли, если это 1, и ничего, если это 0. Итак, вы хотите каким-то образом превратить 1 в -1 (т.е. 0xffffffff), а 0 останется прежним. Это просто -X .... так что вы делаете ...

Y & (-X)

для каждого элемента ... работа выполнена?

Edit2: Чтобы дать пример кода, вы можете сделать что-то вроде этого и избежать ветки:

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

Конечно, лучше всего хранить единицы и нули в массиве целых чисел и, следовательно, избегать сдвигов.

Изменить: также стоит отметить, что если значения в Y имеют размер 16 бит, вы можете выполнять 2 из них и операции за операцию (4, если у вас 64-битные регистры). Однако это означает преобразование значений X 1 на 1 в большее целое число.

т.е. YVals = -4, 3 в 16-битном формате = 0xFFFC, 0x3 ... вставить в 1 32-битный, и вы получите 0xFFFC0003. Если у вас есть 1, 0 в качестве значений X, вы формируете битовую маску из 0xFFFF0000 и 2 вместе, и у вас есть 2 результата в 1 побитовой операции и операции.

Другое редактирование:

ЕСЛИ вам нужен код, описывающий, как выполнять второй метод, что-то вроде , это должно сработать (хотя он использует неопределенное поведение, поэтому он может работать не на всех компиляторах .. работает на всех компиляторах, с которыми я сталкивался).

union int1632
{
     int32_t i32;
     int16_t i16[2];
};

int result=0;
for ( int i = 0; i < (L & ~0x1); i += 2 )
{
    int3264 y3264;
    y3264.i16[0] = Y[i + 0];
    y3264.i16[1] = Y[i + 1];

    int3264 x3264;
    x3264.i16[0] = -(int16_t)((X >> (i + 0)) & 1);
    x3264.i16[1] = -(int16_t)((X >> (i + 1)) & 1);

    int3264 res3264;
    res3264.i32  = y3264.i32 & x3264.i32;

    result += res3264.i16[0] + res3264.i16[1];    
}

if ( i < L )
    result+=Y[i] & -(int)((X >> i) & 1);

Будем надеяться, что компилятор оптимизирует присваивания (я не уверен, но идею можно было бы переработать так, чтобы они точно были) и немного ускорится, поскольку теперь вам нужно только сделать 1 побитовое - и вместо 2. Ускорение будет незначительным ...