Application à Karatsuba¶
Dans cette partie, nous allons voir l’intérêt du jeu d’instructions AVX en s’intéressant à la multiplication polynomiale dans l’anneau \({\mathbb F}_2[x]\).
Rappelons que si \(A(x)\) et \(B(x)\) sont deux polynômes de degré au plus \(n-1\) de \({\mathbb R}[x]\), on a alors
\(A(X)B(X) = \displaystyle\sum_{k=0}^{2n-2}\bigg (\sum_{i=0}^k a_ib_{k-i}\bigg )x^k\).
Cette formule fait intervenir des calculs dans lesquels des propagations de retenues doivent être gérées (somme interne). En travaillant dans \({\mathbb F}_2[x]\), les coefficients des polynômes ne peuvent être égaux qu’à 0 ou 1 et les additions de la formule ci-dessus se transforment en addition modulo 2, ce qui correspond à l’opérateur logique nommé ou-exclusif, noté \(\oplus\). La formule devient
\(A(X)B(X) = \displaystyle\bigoplus_{k=0}^{2n-2}\bigg (\bigoplus_{i=0}^k a_ib_{k-i}\bigg )x^k\),
dans laquelle aucune retenue (Carry Less) n’est à gérer puisque 1+1=0 !
En 2010, Intel introduit dans ses processeurs l’instruction PCLMULQDQ (Perform Carry Less MULtiplication of one QuadworD bye one Quadword) permettant de réaliser le produit dans \({\mathbb F}_2[x]\) de 2 polynômes de degré au plus 63. Chaque polynôme est représenté par un entier de 64 bits (le bit \(i\) vaut 1 si \(x^i\) est un monôme du polynôme). Ces 2 entiers doivent être stockés dans une variable de type __m128i, soit en partie haute, soit en partie basse. Un paramètre de l’instruction de multiplication permet de spécifier quelle partie contient les polynômes. Le résultat est stocké dans un registre de 128 bits :
__m128i _mm_clmulepi64_si128 (__m128i A, __m128i B, const int imm8)
Les registres A et B contiennent 2 mots de 64 bits :
127 63 0 127 63 0
|--------------------|--------------------| |--------------------|--------------------|
A: |________A_1_________|________A_0_________| B: |________B_1_________|________B_0_________|
64 64
Les bits 0 et 4 de la variable imm8 indiquent à l’instruction quelle multiplication effectuer:
- si imm8=0, l’opération effectuée est \(A_0(x)B_0(x)\),
- si imm8=1, l’opération effectuée est \(A_1(x)B_0(x)\),
- si imm8=16, l’opération effectuée est \(A_0(x)B_1(x)\),
- si imm8=17, l’opération effectuée est \(A_1(x)B_1(x)\).
La multiplication de Karatsuba (version 128 bits)¶
Nous allons utiliser cette instruction afin d’étendre la multiplication à des polynômes de degré au plus 127. Soient \(A(X)\) et \(B(x)\) deux polynômes de degré au plus 127, on peut toujours écrire
\(A(x)=A_1(x)x^{64}+A_0(x)\) et \(B(x)=B_1(x)x^{64}+B_0(x)\)
où \(A_i(x)\) et \(B_i(x)\) sont deux polynômes de degré au plus 63. Les polynômes \(A(x)\) et \(B(x)\) sont stockés dans 2 registres de 128 bits, les parties hautes et basses représentant respectivement les polynômes \(A_1(x)\), \(A_0(x)\), \(B_1(x)\) et \(B_0(x)\). La multiplication de Karatsuba permet de réaliser le produit \(C(x)=A(x)B(x)\) en utilisant uniquement 3 multiplications:
\(\scriptsize C(x)= A_1(x)B_1(x)x^{128} + ((A_1(x)+A_0(x))(B_1(x)+B_0(x))-A_1(x)B_1(x)-A_0(x)B_0(x))x^{64} +A_0(x)B_0(x)\)
Le résultat sera stocké dans 2 registres de 128 bits (nous aurions pu aussi utiliser un registre de 256 bits). Voici le code complet commenté.
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <immintrin.h>
/******************************************************************************************/
/* */
/* Programme pour réaliser la multiplication de 2 polynômes de degré au plus 127 */
/* en utilisant la méthode de Karatsuba. Un polynôme est codé sur un registre de 128 bits */
/* le résultat est renvoyé sur 2 registres de 128 bits */
/* */
/******************************************************************************************/
typedef union
{
__m128i vec128;
uint64_t vec64[2];
uint32_t vec32[4];
uint16_t vec16[8];
uint8_t vec8[16];
} __reg128;
void print128(char *s, int mode, __reg128 A)
{
switch(mode)
{
case 8 :
printf("%s=(",s);
for(int i =0; i < 15; i++) printf("'%2x',",A.vec8[i]);
printf("'%2x')\n",A.vec8[15]);
break;
case 16 :
printf("%s=(",s);
for(int i =0; i < 7; i++) printf("'%4x',",A.vec16[i]);
printf("'%4x')\n",A.vec16[7]);
break;
case 32 :
printf("%s=(",s);
for(int i =0; i < 3; i++) printf("'%8x',",A.vec32[i]);
printf("'%8x')\n",A.vec32[3]);
break;
case 64 :
printf("%s=('0x%lx','0x%lx')\n",s,A.vec64[0],A.vec64[1]);
}
}
void karat_mult_F2(__m128i *C, __m128i A, __m128i B) {
// A(x) = A[1]x^(64) + A[0]
// B(x) = B[1]x^(64) + B[0]
// C(X) = A[1]B[1]x^(128) + ((A[1]+A[0])(B[1]+B[0])-A[1]B[1]-A[0]B[0])x^(64) + A[0]B[0]
// = C[1]x^(128)+C[0]
//
// 255 191 127 63 0
// |--------------------|--------------------|--------------------|--------------------|
// |____________________|____________________|____________________|____________________|
// 192 128 64
// |-----------------------------------------||----------------------------------------|
// A[1]B[1] A[0][B0]
// |----------------------------------------|
// (A[1]+A[0])(B[1]+B[0])-A[1]B[1]-A[0]B[0]
//
// Calcul de A[0]B[0]
//
// On utilise l'intruction clmul permettant de multiplier 2 polynômes de degré <= 63 sur F_2
// le paramètre 0 permet de spécifier que les 2 polynômes se trouvent dans les 64 premiers
// bits de A et B.
// On place ce premier résultat qui vaut A[0]B[0] dans C[0]
C[0] = _mm_clmulepi64_si128(A, B, 0);
//
// On utilise l'intruction clmul permettant de multiplier 2 polynômes de degré <= 63 sur F_2
// le paramètre 17 permet de spécifier que les 2 polynômes se trouvent dans les 64 derniers
// bits de A et B.
// On place ce résultat qui vaut A[1]B[1] dans C[1]
C[1] = _mm_clmulepi64_si128(A, B, 17);
// On calcule la partie du milieu
// Pour faire le xor entre A1 et A0 et B1 et B0, on construit à la volée à partir de A=(A1 A0)
// le registe (A0 A1) (swap des mots de 64 bits) et on fait la même chose pour B --> (B0 B1)
// on fera alors un xor entre A et (A0 A1) --> (A1^A0 A0^A1) et de même pour B et on
// récupérera les 64 derniers bits de ces deux résultats pour appliquer clmul afin de calculer
// (A[1]+A[0])(B[1]+B[0]). Ce swap se fait avec l'instruction _mm_shuffle_epi32 qui permet de
// permuter les paquets de 32 bits.
// Numérotons ces paquets de 0 (32 bits de poids faible) à 3 (32 bits de poids fort).
// L'instruction _mm_shuffle_epi32 utilise un paramètre qui code la permutation à utiliser.
// Ici la permutation à utiliser est 3 2 1 0 --> 1 0 3 2 (les 32 bits de pois faible sont
// positionnés dans le bloc 2 et les 32 bits suivants dans le bloc 3. Ainsi les 64 bits de
// poids faible constituent maintenant les 64 bits de poids fort du résultat et ainsi de suite.
// Chaque valeur de la permutation est codée sur 2 bits
// 1 -> 01, 0 -> 00, 3-> 11, 2->10, ce qui donne 0100 1110, soit 0x4e en hexadécimal.
// A = (A1 A0) = (A11 A10 A01 A00)
// Aij 32 bits ---> _mm_shuffle_epi32(A,0x4e) = (A01 A00 A11 A10)
// _mm_shuffle_epi32(A,0x4e) --> (A0 A1)
// _mm_xor_si128(A,_mm_shuffle_epi32(A,0x4e)) --> (A1^A0, A0^A1)
// Calcul de (A[1]+A[0])(B[1]+B[0])
__m128i middle = _mm_clmulepi64_si128(_mm_xor_si128(A,_mm_shuffle_epi32(A,0x4e)),
_mm_xor_si128(B,_mm_shuffle_epi32(B,0x4e)),0);
// Calcul de (A[1]+A[0])(B[1]+B[0])-A[1]B[1]
middle = _mm_xor_si128(middle,C[1]);
// Calcul de (A[1]+A[0])(B[1]+B[0])-A[1]B[1]-A[0]B[0]
middle = _mm_xor_si128(middle,C[0]);
// On récupère les 64 bits de poids faible du milieu et on xor sur les 64 bits de poids fort
// de C[0]. La fonction unpacklo créé un entier de 128 bits en plaçant les 64 bits de poids
// faible de l'argument 1 dans les 64 bits de poids faible du résultat et les 64 bits de poids
// faible de l'argument 2 dans les 64 bits de poids fort du résultat
C[0] = _mm_xor_si128(C[0], _mm_unpacklo_epi64(_mm_setzero_si128(), middle));
// On récupère les 64 bits de poids fort du milieu et on xor sur les 64 bits de poids faible
// de C[1]. La fonction unpackhi créé un entier de 128 bits en plaçant les 64 bits de poids
// fort de l'argument 1 dans les 64 bits de poids faible du résultat et les 64 bits de poids
// forts de l'argument 2 dans les 64 bits de poids fort du résultat
C[1] = _mm_xor_si128(C[1], _mm_unpackhi_epi64(middle, _mm_setzero_si128()));
}
int main(void)
{
__reg128 tempo;
__m128i *C = calloc(2,sizeof(__reg128));
// on initialise les registres 128 bits
// {64 bits poids faible, 64 bits poids fort}
__m128i A = {0xffffaa1256ee1234,0xfffabfffeeffffff};
__m128i B ={0xea0d362010800099,0xbfeefffdffffffff};
tempo = (__reg128)A;
printf("Données : \n");
printf("A0=0x%lx\n",tempo.vec64[0]);
printf("A1=0x%lx\n",tempo.vec64[1]);
tempo = (__reg128)B;
printf("B0=0x%lx\n",tempo.vec64[0]);
printf("B1=0x%lx\n\n",tempo.vec64[1]);
karat_mult_F2(C, A, B);
tempo = (__reg128)C[0];
printf("Resultat : \n");
print128("C0",64,tempo);
tempo = (__reg128)C[1];
print128("C1",64,tempo);
}
Voici un code Sage permettant de vérifier l’exactitude des résultats.
def tolistebin(h):
L=[]
for i in range(64):
L += [h & 1]
h >>= 1
return L
# coller ici les données du programme C
A0=
A1=
B0=
B1=
# coller ici les résulats du programme C
C0=
C1=
# on construit l'anneau des polynômes GF2[x]
R.<x> = PolynomialRing(GF(2))
AL = tolistebin(A0)
AH = tolistebin(A1)
BL = tolistebin(B0)
BH = tolistebin(B1)
A = R(AL)+(x^(64))*R(AH)
B = R(BL)+(x^(64))*R(BH)
C = A*B
print(C)
C=list(C)
# on complete par des 0 pour avoir 256 bits
C = C+[0]*(256-len(C))
# la liste est ordonnée des bits de poids faible
# vers les bits de poids fort
# int transforme une chaine de caractères en hexa
# en considérant que le bit de poids faible est le
# dernier caractère de la chaine, il faut donc
# d'abord renverser la liste
result_sage_C0 = (hex(int(''.join(list(map(str,C[0:64][::-1]))),2)),hex(int(''.join(list(map(str,C[64:128][::-1]))),2)))
result_sage_C1 = (hex(int(''.join(list(map(str,C[128:192][::-1]))),2)),hex(int(''.join(list(map(str,C[192:256][::-1]))),2)))
print('C0 : ', result_sage_C0)
print('C1 : ', result_sage_C1)
# on vérifie que les résultats renvoyés par le programme C
# sont identiques aux résultats renvoyés par Sage
print(result_sage_C0 == C0)
print(result_sage_C1 == C1)
La multiplication de Karatsuba (version 256 bits)¶
La fonction précédente peut nous servir de base pour obtenir une multiplication polynomiale pour des polynômes de degré au plus 255. Dans ce cas, les polynômes \(A(X)\) et \(B(x)\) seront stockés dans des registres 256 bits. On utilisera la structure suivante pour faciliter la gestion de ces registres :
typedef union
{
__m256i vec256;
__m128i vec128[2];
uint64_t vec64[4];
uint32_t vec32[8];
uint16_t vec16[16];
uint8_t vec8[32];
} __reg256;
Le code ci-dessous réalise la multiplication polynomiale :
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <immintrin.h>
/******************************************************************************************/
/* */
/* Programme pour réaliser la multiplication de 2 polynômes de degré au plus 255 */
/* en utilisant la méthode de Karatsuba. Un polynôme est codé sur un registre de 256 bits */
/* le résultat est renvoyé sur 2 registres de 256 bits */
/* */
/******************************************************************************************/
typedef union
{
__m256i vec256;
__m128i vec128[2];
uint64_t vec64[4];
uint32_t vec32[8];
uint16_t vec16[16];
uint8_t vec8[32];
} __reg256;
void print256(char *s, int mode, __reg256 A)
{
switch(mode)
{
case 8 :
printf("%s=(",s);
for(int i =0; i < 31; i++) printf("'%2x',",A.vec8[i]);
printf("'%2x')\n",A.vec8[31]);
break;
case 16 :
printf("%s=(",s);
for(int i =0; i < 15; i++) printf("'%4x',",A.vec16[i]);
printf("'%4x')\n",A.vec16[15]);
break;
case 32 :
printf("%s=(",s);
for(int i =0; i < 7; i++) printf("'%8x',",A.vec32[i]);
printf("'%8x')\n",A.vec32[7]);
break;
case 64 :
printf("%s=('0x%lx','0x%lx', '0x%lx','0x%lx')\n",s,A.vec64[0],A.vec64[1], A.vec64[2], A.vec64[3]);
break;
case 128 :
printf("%s=('0x%lx%lx', '0x%lx%lx')\n",s,A.vec64[1],A.vec64[0], A.vec64[3], A.vec64[2]);
}
}
/* Multiplication pour les polynomes de degre 127 */
void karat_mult_F2(__m128i *C, __m128i A, __m128i B) {
C[0] = _mm_clmulepi64_si128(A, B, 0);
C[1] = _mm_clmulepi64_si128(A, B, 17);
__m128i middle = _mm_clmulepi64_si128(_mm_xor_si128(A,_mm_shuffle_epi32(A,0x4e)),
_mm_xor_si128(B,_mm_shuffle_epi32(B,0x4e)),0);
middle = _mm_xor_si128(middle,C[1]);
middle = _mm_xor_si128(middle,C[0]);
C[0] = _mm_xor_si128(C[0], _mm_unpacklo_epi64(_mm_setzero_si128(), middle));
C[1] = _mm_xor_si128(C[1], _mm_unpackhi_epi64(middle, _mm_setzero_si128()));
}
void karat_mult_F2_256(__m256i *C, __reg256 A, __reg256 B) {
__reg256 middle;
// A(x) = A[1]x^(128) + A[0]
// B(x) = B[1]x^(128) + B[0]
// C(X) = A[1]B[1]x^(256) + ((A[1]+A[0])(B[1]+B[0])-A[1]B[1]-A[0]B[0])x^(128) + A[0]B[0]
// = C[1]x^(256)+C[0]
//
//
// Calcul de A[0]B[0] en appelant la fonction
// karat_mult_F2, le résultat est placé dans C[0]
karat_mult_F2((__m128i *)C,A.vec128[0],B.vec128[0]);
// Calcul de A[1]B[1]
karat_mult_F2((__m128i *)(C+1),A.vec128[1],B.vec128[1]);
// Calcul de (A[1]+A[0])(B[1]+B[0])
karat_mult_F2(middle.vec128, _mm_xor_si128(A.vec128[1],A.vec128[0]), _mm_xor_si128(B.vec128[1],B.vec128[0]));
// Calcul de (A[1]+A[0])(B[1]+B[0])-A[1]B[1]
middle.vec256 = _mm256_xor_si256(middle.vec256,C[1]);
// Calcul de (A[1]+A[0])(B[1]+B[0])-A[1]B[1]-A[0]B[0]
middle.vec256 = _mm256_xor_si256(middle.vec256,C[0]);
// On récupère les 128 bits de poids faible du milieu et on xor sur les 128 bits de poids fort
// de C[0]. La fonction _mm256_set_m128i créé un entier de 256 bits à partir de 2 paramètres
// de 128 bits.
C[0] = _mm256_xor_si256(C[0], _mm256_set_m128i(middle.vec128[0],_mm_setzero_si128 ()));
// On récupère les 128 bits de poids fort du milieu et on xor sur les 128 bits de poids faible
// de C[1].
C[1] = _mm256_xor_si256(C[1], _mm256_set_m128i(_mm_setzero_si128 (), middle.vec128[1]));
}
int main(void)
{
__reg256 tempo, A, B;
__m256i *C = calloc(2,sizeof(__m256i));
// on initialise les registres 256 bits
// des poids faibles vers les poids forts
A.vec256 = (__m256i){0xffffaa1256ee1234,0xfffabfffeeffffff,0xffffaa1256ee0000,0xfffabfffee111111};
B.vec256 = (__m256i){0xea0d362010800099,0xbfeefffdffffffff,0xea0d362010811199,0xbfee00000000};
printf("Données : \n");
print256("A",128,A);
print256("B",128,B);
karat_mult_F2_256(C, A, B);
tempo = (__reg256)C[0];
printf("Resultat : \n");
print256("C0",128,tempo);
tempo = (__reg256)C[1];
print256("C1",128,tempo);
}
Le code Sage de vérification :
def tolistebin(h):
L=[]
for i in range(128):
L += [h & 1]
h >>= 1
return L
# coller ici les données du programme C
A=
B=
# coller ici les résulats du programme C
C0=
C1=
# on construit l'anneau des polynômes GF2[x]
R.<x> = PolynomialRing(GF(2))
AL = tolistebin(int(A[0],16))
AH = tolistebin(int(A[1],16))
BL = tolistebin(int(B[0],16))
BH = tolistebin(int(B[1],16))
A = R(AL)+(x^(128))*R(AH)
B = R(BL)+(x^(128))*R(BH)
C = A*B
print(C)
C=list(C)
# on complete par des 0 pour avoir 512 bits
C = C+[0]*(512-len(C))
# la liste est ordonnée des bits de poids faible
# vers les bits de poids fort
# int transforme une chaine de caractères en hexa
# en considérant que le bit de poids faible est le
# dernier caractère de la chaine, il faut donc
# d'abord renverser la liste
result_sage_C0 = (hex(int(''.join(list(map(str,C[0:128][::-1]))),2)),hex(int(''.join(list(map(str,C[128:256][::-1]))),2)))
result_sage_C1 = (hex(int(''.join(list(map(str,C[256:384][::-1]))),2)),hex(int(''.join(list(map(str,C[384:len(C)][::-1]))),2)))
print('C0 : ', result_sage_C0)
print('C1 : ', result_sage_C1)
# on vérifie que les résultats renvoyés par le programme C
# sont identiques aux résultats renvoyés par Sage
print(result_sage_C0 == C0)
print(result_sage_C1 == C1)