Ответ 1

Для этого вы можете использовать малоизвестные PMADDUBSW (Multiply и Add Packed Signed and 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 раз быстрее, чем простой способ (тестируется с 2 Кбайт случайных данных на 4770 КБ, беря максимум из ста пробегов для каждого).

Использование указателей, предложенных другим ответом, улучшает его дальше, примерно в 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;
}

Ответ 2

Этот цикл,

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

может получить выгоду от перехода на:

char *b = a + top;
a++;
for (; a < b; ) 
{
    sum += ( *a++ ) * ( *b--);
}

Удалив [], вы сохраняете много арифметики для каждого доступа к массиву. Это уменьшает теоретическое число адресных вычислений: 4 с ++m --k и a[m] a[k] до 2 с *a++ *b--

Приращение простого указателя дешевле и обычно быстрее общего, поскольку доступ к массиву не всегда оптимизирован с помощью [].

Надеюсь, что это поможет

Ответ 3

Операция, которую вы хотите выполнить, называется дискретной сверткой, она появляется при умножении больших чисел. Наивный алгоритм, который вы используете, имеет сложность O (n & thinsp; ²), но решение O (n log n) получается с помощью дискретное преобразование Фурье.

Дискретные конволюции

Дискретная свертка c = a & lowast; b двух последовательностей a = a ₀, a ₁, & hellip;, a _{n & minus; 1} и b = b ₀, b ₁, & hellip;, b _{n & minus; 1} с n элементами каждая представляет собой последовательность из 2 n & minus; 1 элементов, определенных для каждого k в виде суммы:

c _i= & sum; _{max (0, я & minus; n + 1) & le; j < min (n, я + 1)} a _{j & thinsp;} b _{я & minus; Jсуб >}

Если предположить, что a _i= b _i= 0 для я & notin; {0, & hellip;, n & minus; 1}, то мы можем упростить это и суммировать по целому числу i:

c _i= & sum; _j a _j & thinsp; b _{я & minus; Jсуб >}

Обратите внимание, как это операция, которую вы хотите выполнить: a = b = a[], а i -я итерация вашего цикла, sum - это просто c _i.

Дискретная свертка хорошо изучена и появляется в ряде математических задач, связанных с обработкой сигналов, комбинаторика и статистикой. К счастью, его можно вычислить в сверхлинейном времени O (n log n) вместо наивного O (n ²).

Циклические дискретные конволюции

Мы можем разложить b на b _N такие, что

b _k= b _{k mod n, & emsp; или
b N= b 0, b 1, & hellip;, b n & minus; 1, b 0, & hellip;, b n & minus; 1, b 0, & hellip;}

Это называется циклическим расширением b до b _N. Дискретная свертка a & lowast; b _N называется дискретной циклической сверткой a и b.

Заметим, что нециклическая дискретная свертка a и b может быть вычислена из циклической свертки путем добавления достаточно большого числа нулей к a и b до свертывания, так что циклическое разложение b не изменяет результат. Подробнее см. в этой статье.

Дискретное преобразование Фурье

Дискретное преобразование Фурье (DFT) преобразует последовательность выборок a в комплексный частотный спектр ℱ (a) дискретизированного сигнала. Преобразование Фурье обратимо и вычислимо в суперлинейном времени O (n log n) с использованием различных алгоритмов быстрого преобразования Фурье (FFT). Это преобразование имеет много применений в обработке сигналов. Он имеет полезное свойство, что

a & lowast; bN= ℱ ^{& minus; 1} (ℱ (a) & times; ℱ (b))

То есть циклическая свертка двух последовательностей a и b является умножением по элементам при преобразовании Фурье. Так как элементарное умножение может быть выполнено в линейном времени, это позволяет вычислить результат, который вам нужен в общем времени O (n log n).

Реализация

Проект FFTW обеспечивает высоко оптимизированные реализации быстрого преобразования Фурье.

Проблемы

Преобразования Фурье работают над комплексными числами и требуют, чтобы блок с плавающей точкой был быстрым. Результат, скорее всего, будет немного неточным. Насколько мне известно, точные методы возможны с помощью некоторой модульной магии (это то, что используется в алгоритме Шенхаге-Штрассена для многоцелевого умножения), но я не уверен в деталях.