Próbowałem napisać funkcję mnożenia carry-less, bez dzielenia przez wielomian, samo mnożenie, wykorzystując instrukcje PCLMULQDQ:
__int128 carry_less(__int128 a, __int128 b)
{
__m128i c1_c0 = _mm_clmulepi64_si128(_mm_set_epi64x(0, a >> 64), _mm_set_epi64x(0, b >> 64), 0);
__m128i d1_d0 = _mm_clmulepi64_si128(_mm_set_epi64x(0, (uint64_t)a), _mm_set_epi64x(0, (uint64_t)b), 0);
__m128i e1_e0 = _mm_clmulepi64_si128(_mm_set_epi64x(0, (uint64_t)a ^ (a >> 64)), _mm_set_epi64x(0, (uint64_t)b ^ (b >> 64)), 0);
uint64_t y4 = _mm_extract_epi64(c1_c0, 1);
uint64_t y3 = _mm_cvtsi128_si64(c1_c0) ^ _mm_extract_epi64(c1_c0, 1) ^ _mm_extract_epi64(d1_d0, 1) ^ _mm_extract_epi64(e1_e0, 1);
uint64_t y2 = _mm_extract_epi64(d1_d0, 1) ^ _mm_cvtsi128_si64(c1_c0) ^ _mm_cvtsi128_si64(d1_d0) ^ _mm_cvtsi128_si64(e1_e0);
uint64_t y1 = _mm_cvtsi128_si64(d1_d0);
return y3;
}
Na podstawie publikacji Efficient implementation of the Galois Counter Mode using a carry-less multiplier and a fast reduction algorithm:
https://sci-hub.se/10.1016/j.ipl.2010.04.011
Algorithm 1: Carry-less Karatsuba. Ale nie działa. Mnożenie carry-less dwóch liczb 128-bitowych powinno zwracać wynik 255-bitowy, 64 bitowe y4, y3, y2, y1. Chcę pomnożyć liczby:
__int128 wynik = carry_less(331755698128909410983678231014631813416, 36893488147419103231);
A zatem y3 powinno zwracać jakąś wartość, ale y4 będzie równe 0. Tymczasem y3 też wychodzi mi równe 0, co jest błędem.
Znalazłem implementację Intela takiego mnożenia, ale razem z dzieleniem przez wielomian (a ja chcę samo mnożenie):
__m128i gfmul(__m128i a, __m128i b)
{
__m128i tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, tmp9, tmp10, tmp11, tmp12;
__m128i XMMMASK = _mm_setr_epi32(0xffffffff, 0x0, 0x0, 0x0);
tmp3 = _mm_clmulepi64_si128(a, b, 0x00);
tmp6 = _mm_clmulepi64_si128(a, b, 0x11);
tmp4 = _mm_shuffle_epi32(a, 78);
tmp5 = _mm_shuffle_epi32(b, 78);
tmp4 = _mm_xor_si128(tmp4, a);
tmp5 = _mm_xor_si128(tmp5, b);
tmp4 = _mm_clmulepi64_si128(tmp4, tmp5, 0x00);
tmp4 = _mm_xor_si128(tmp4, tmp3);
tmp4 = _mm_xor_si128(tmp4, tmp6);
tmp5 = _mm_slli_si128(tmp4, 8);
tmp4 = _mm_srli_si128(tmp4, 8);
tmp3 = _mm_xor_si128(tmp3, tmp5);
tmp6 = _mm_xor_si128(tmp6, tmp4);
tmp7 = _mm_srli_epi32(tmp6, 31);
tmp8 = _mm_srli_epi32(tmp6, 30);
tmp9 = _mm_srli_epi32(tmp6, 25);
tmp7 = _mm_xor_si128(tmp7, tmp8);
tmp7 = _mm_xor_si128(tmp7, tmp9);
tmp8 = _mm_shuffle_epi32(tmp7, 147);
tmp7 = _mm_and_si128(XMMMASK, tmp8);
tmp8 = _mm_andnot_si128(XMMMASK, tmp8);
tmp3 = _mm_xor_si128(tmp3 , tmp8);
tmp6 = _mm_xor_si128(tmp6, tmp7);
tmp10 = _mm_slli_epi32(tmp6, 1);
tmp3 = _mm_xor_si128(tmp3, tmp10);
tmp11 = _mm_slli_epi32(tmp6, 2);
tmp3 = _mm_xor_si128(tmp3, tmp11);
tmp12 = _mm_slli_epi32(tmp6, 7);
tmp3 = _mm_xor_si128(tmp3, tmp12);
__m128i result = _mm_xor_si128(tmp3, tmp6);
return result;
}
To chyba działa poprawnie.
https://www.intel.com/content/dam/develop/external/us/en/documents/clmul-wp-rev-2-02-2014-04-20.pdf
Figure 7, algorytm 2 i 4. Tymczasem ja chcę tylko algorytm 2 (jest na stronie 13). Nie rozumiem tego kodu i nie potrafię wyodrębnić z niego algorytmu 2. Instrukcji mnożenia PCLMULQDQ nauczyłem się stąd:
https://wunkolo.github.io/post/2020/05/pclmulqdq-tricks/
Ale tak jak pisałem, wygląda na to, że źle ich użyłem w moim kodzie. Nie wiem co to jest _mm_shuffle_epi32(a, 78) i ciężko znaleźć coś klarownego na ten temat. Czy w tym kodzie Intela można łatwo oddzielić mnożenie od algorytmu 4, który wykonuje dziele przez wielomian? A może ktoś potrafi poprawić mój kod?
PS Tu są pewne wyjaśnienia SIMD intrinsics and instructions:
https://arxiv.org/pdf/1503.03465.pdf
Czyli mm shuffle_epi32 jakoś miesza bajty, ale po co to jest w ogóle robione w tym algorytmie i co oznacza 78 nadal nie wiem.
PS 2 Wydaje mi się, że algorytm mnożenia kończy się tutaj:
tmp3 = _mm_clmulepi64_si128(a, b, 0x00);
tmp6 = _mm_clmulepi64_si128(a, b, 0x11);
tmp4 = _mm_shuffle_epi32(a, 78);
tmp5 = _mm_shuffle_epi32(b, 78);
tmp4 = _mm_xor_si128(tmp4, a);
tmp5 = _mm_xor_si128(tmp5, b);
tmp4 = _mm_clmulepi64_si128(tmp4, tmp5, 0x00);
tmp4 = _mm_xor_si128(tmp4, tmp3);
tmp4 = _mm_xor_si128(tmp4, tmp6);
tmp5 = _mm_slli_si128(tmp4, 8);
tmp4 = _mm_srli_si128(tmp4, 8);
tmp3 = _mm_xor_si128(tmp3, tmp5);
tmp6 = _mm_xor_si128(tmp6, tmp4);
Bo dalej są już przesunięcia, które są charakterystyczne dla algorytmu 4, czyli dzielenia przez wielomian. Nie wiem tylko, które tmp są ostatecznym wynikiem. Chyba dwa ostatnie.
PS 3 Zrobiłem to tak:
__int128 carry_less(__m128i a, __m128i b)
{
__m128i tmp3, tmp4, tmp5, tmp6;
tmp3 = _mm_clmulepi64_si128(a, b, 0x00);
tmp6 = _mm_clmulepi64_si128(a, b, 0x11);
tmp4 = _mm_shuffle_epi32(a, 78);
tmp5 = _mm_shuffle_epi32(b, 78);
tmp4 = _mm_xor_si128(tmp4, a);
tmp5 = _mm_xor_si128(tmp5, b);
tmp4 = _mm_clmulepi64_si128(tmp4, tmp5, 0x00);
tmp4 = _mm_xor_si128(tmp4, tmp3);
tmp4 = _mm_xor_si128(tmp4, tmp6);
tmp5 = _mm_slli_si128(tmp4, 8);
tmp4 = _mm_srli_si128(tmp4, 8);
tmp3 = _mm_xor_si128(tmp3, tmp5);
tmp6 = _mm_xor_si128(tmp6, tmp4);
uint64_t Lo_first = _mm_extract_epi64(tmp3, 1);
uint64_t Hi_first = _mm_cvtsi128_si64(tmp3);
uint64_t Lo_second = _mm_cvtsi128_si64(tmp6);
uint64_t Hi_second = _mm_extract_epi64(tmp6, 1);
return Lo_second;
}
Lo_first to niskie 64-bity pierwszej 128-bitowej połówki (wynik jest 256 bitowy). Hi_second to najwyższe 64-bity, drugiej 128-bitowej połówki.
Przy próbie mnożenia losowej liczby z 2^65-1:
wynik = carry_less(179854719828577785056322020520430909555, _mm_setr_epi32(4294967295, 4294967295, 1, 0));
Najwyższe bity są zgodnie z oczekiwaniami równe zero. Pytanie, czy dobrze inicjuję __m128i jako _mm_setr_epi32(4294967295, 4294967295, 1, 0), czyli wpisuję liczby w odwróconym porządku, pierwsza liczba 32-bitowa stanowi najniższe 32 bity.
Widzę, że można to normalnie zainicjować unikając "r":
_mm_set_epi32(0, 1, 4294967295, 4294967295)
A nawet używając liczb 64-bitowych:
_mm_set_epi64x(1, 18446744073709551615)