uint32_t ADLER32_update_(uint32_t adler, uint8_t *p, size_t len){
const uint32_t BASE = 65521;
uint32_t s1 = (adler ) & 0xFFFF;
uint32_t s2 = (adler >> 16);
if (len) do{
s1 += *p++;
if (s1 >= BASE) s1 -= BASE;
s2 += s1;
if (s2 >= BASE) s2 -= BASE;
} while (--len);
return (s2<<16) + s1;
}
Статья описывает процесс оптимизации алгоритма для векторного-паралельного вычисления контрольной суммы.
Оптимизация - это последовательные тождественные преобразования алгоритма.
const int Nmax = 5552;
do {
int n = len>Nmax? Nmax: len;
len-=n;
do {
s1 += *p++;
s2 += s1;
}while (--n);
// не полное редуцирование
s1-= (s1>>16)*0xFFF1u;
s2-= (s2>>16)*0xFFF1u;
} while (len);
s1 = mod65521(s1);// s1%=BASE;
s2 = mod65521(s2);// s2%=BASE;
s2 = s2 + s1*n + p0*(n) + p1*(n-1) ... MAX(s2) = (n+1)*(0xFFF0 + 0xFF*n/2)
В оригинальном алгоритме Марка Адлера используется операция остаток от деления, на простое число 65521. Вычисление остатка можно упростить. В предложенном алгоритме производится не полное редуцирование старших 16 бит, это недоделанная операция остатка от деления - не полный остаток.
Финальное редуцирование можно оптимизировать за счет использования обратной операции к делению. Вместо деления выполняем умножение на обратное число со сдвигом. Способ вычисления обратного числа рассматривал в одной из статей.
#define mod65521(x) ((x) - 0xFFF1*(unsigned long)(((x)*0x80078071ULL)>>47))
// Оптимизация вложенного цикла
s2+=s1*n;
do {
s1 += p[0];
s2 += p[0]*n;
p++;
} while (--n);
В таком варианте компилятор (GCC 10+ и clang) векторизует цикл и выполняет вычисления на векторе 16 байт. Вопреки ожиданиям, оптимизированный компилятором код довольно громоздкий и не дает ожидаемого ускорения в 16 раз, поскольку выполняется перепаковка входных дынных в вектор uint32х4. По этой причине дальнейшая оптимизация выполняется вручную за счет использования векторных инструкций горизонтального сложения для суммирования элементов (s1). Горизонтальное сложение и умножение на векторе можно выполнить с использованием инструкции _mm256_maddubs_epi16, а горизонтальное суммирование по вектору можно выполнить на инструкции _mm_sad_epu8
Ниже показываю, как разворачивается и сворачивается цикл. Начнем с шага 2.
// Развертывание вложенного цикла
do {
s3 += s1;
s1 += p[0] + p[1];
s2 += p[0]*2 + p[1]*1;
p+=2;
} while ((n-=2)>0);
s2+=s3*2;
// Свертка по L-элементов
do {
s3 += s1;
for (i=0;i<L; i++) {
s1 += p[i];
s2 += p[i]*(L-i);
}
p+=L;
} while ((n-=L)>0);
s2+=s3*L;
Тоже на векторных инструкциях Intel AVX2.
do{
vs3 = _mm256_add_epi32(vs3, vs1);
__m256i v = _mm256_lddqu_si256((void*)p); p+=32;
__m256i v1 = _mm256_maddubs_epi16(v, E);
__m256i v2 = _mm256_maddubs_epi16(v, M);
v1 = _mm256_madd_epi16 (v1, _mm256_set1_epi16(1));
v2 = _mm256_madd_epi16 (v2, _mm256_set1_epi16(1));
vs1 = _mm256_add_epi32(vs1,v1);
vs2 = _mm256_add_epi32(vs2,v2);
} while(--n);
Тоже на инструкциях Intel AVX512_VNNI (инструкции для сверточных нейро сетей подходят для наших целей).
do {
__m256i v = _mm256_lddqu_si256((void*)p); p+=32;
vs3 = _mm256_add_epi32(vs3, vs1);
vs1 = _mm256_dpbusd_epi32(vs1, v, E);
vs2 = _mm256_dpbusd_epi32(vs2, v, M);
} while((n-=32)>0);
Я применил инструкции разрядностью 256 бит, поскольку производительность одной инструкции 512 бит такая же, как двух инструкций 256бит. Инструкции _mm256_dpbusd_epi32 создают задержку в пять тактов, поэтому есть дополнительная возможность ускорения, если разложить расчет на пять-десять независимых переменных.
Быстрое редуцирование на векторе удалость описать таким образом:
vs2 = _mm256_sub_epi32 (vs2, _mm256_mullo_epi32(_mm256_srli_epi32(vs2, 16), P));
Мы показали порядок вывода паралельного алгоритма из последовательного. Теперь перейдем к анализу быстродействия. Анализ выполняется с использованием анализатора машинного кода LLVM-MCA.
# gcc -O3 -march=skylake -S -o test.s adler.c # llvm-mca.exe -mcpu=skylake -timeline test.s Фрагмент отчета (временная диаграмма), начиная со второго цикла: [2,0] . .DeeeeeeeE----------R . vlddqu (%rax), %ymm0 [2,1] . .D================eER . vpaddd %ymm5, %ymm2, %ymm3 [2,2] . .D=eeeeeeeeeeeeE----R . vpmaddubsw .LC10(%rip), %ymm0, %ymm1 [2,3] . .D=eeeeeeeeeeeeE----R . vpmaddubsw .LC11(%rip), %ymm0, %ymm0 [2,4] . . D=====eeeeeeeeeeeeER . vpmaddwd .LC12(%rip), %ymm1, %ymm1 [2,5] . . D=====eeeeeeeeeeeeER . vpmaddwd .LC12(%rip), %ymm0, %ymm0 [2,6] . . D=================eER . vpaddd %ymm2, %ymm1, %ymm1 [2,7] . . D=================eER . vpaddd %ymm4, %ymm0, %ymm0 [2,8] . . DeE----------------R . addq $32, %rax [2,9] . . D===============eE-R . vmovdqa %ymm3, %ymm5 [2,10] . . D=================eER . vmovdqa %ymm1, %ymm2 [2,11] . . D=================eER . vmovdqa %ymm0, %ymm4 [2,12] . . D=eE----------------R . cmpq %rcx, %rax [2,13] . . D==eE---------------R . jne .L21
Из отчета видно, что скорость реализиации алгоритма на инструкциях AVX2 равняется 3 такта на цикл, на 32 байта. Раз в 10 ускорились! И это далеко не предел.
Полностью код доступен на Github
