Solution

Il s’agit de la multiplication polynomiale sur \(\mathbb{F}_2[x]\) pour des polynômes de degré au plus 127 en utilisant la méthode classique schoolbook qui nécessite 4 multiplications intermédiaires que l’on peut réaliser en une seule instruction en AVX 512. Le code 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 schoolbook. Un polynôme est codé sur un registre de 128 bits   */
/* le résultat est renvoyé dans 1 registre de 256 bits.                                   */  
/* Le code tire partie de l'extension VPCLMULQDQ  disponible sur certains processeurs     */
/* AVX512                                                                                 */
/*                                                                                        */
/******************************************************************************************/
 
typedef union
{
  __m128i vec128;
  uint64_t vec64[2];
  uint32_t vec32[4];
  uint16_t vec16[8];
  uint8_t vec8[16];
} __reg128;

typedef union
{
  __m256i vec256;
  __m128i vec128[2];
  uint64_t vec64[4];
  uint32_t vec32[8];
  uint16_t vec16[16];
  uint8_t vec8[32];
} __reg256;

typedef union
{
  __m512i vec512;
  __m256i vec256[2];
  __m128i vec128[4];
  uint64_t vec64[8];
  uint32_t vec32[16];
  uint16_t vec16[32];
  uint8_t vec8[64];
} __reg512;

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(__reg256 *C, __reg128 A, __reg128 B) {

__reg512 aux1,aux2;
__reg128 middle;

 // 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]B[0]+A[0]B[1])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]B[0]+A[0]B[1]
 //
 // On initialise les registres aux1 et aux2 de la façon suivante afin de pouvoir effectuer
 // les 4 multiplications en 1 seule instruction
 // 
 //  511       447       383       319       255       191       127       63        0  
 //  |---------|---------|---------|---------|---------|---------|---------|---------|
 //  |____0____|___A0____|____0____|___A1____|___0_____|___A1____|____0____|___A0____|
 //          448       384       320       256       192       128        64 
 //
 //  511       447       383       319       255       191       127       63        0  
 //  |---------|---------|---------|---------|---------|---------|---------|---------|
 //  |____0____|___B1____|____0____|___B0____|___0_____|___B1____|____0____|___B0____|
 //          448       384       320       256       192       128        64 
 // 
aux1.vec512 = _mm512_set_epi64(0,A.vec64[0],0,A.vec64[1], 0, A.vec64[1],0, A.vec64[0]);
aux2.vec512 = _mm512_set_epi64(0,B.vec64[1],0,B.vec64[0], 0, B.vec64[1],0, B.vec64[0]); 
//
 // On utilise l'instruction clmulepi64_epi128 avec la constante Imm8 égale à 0
 // pour spécifier que les polynômes à multiplier sont dans les 64 bits de poids
 // faible de chaque paquet de 128 bits.
aux1.vec512 = _mm512_clmulepi64_epi128(aux1.vec512, aux2.vec512, 0);
 // 
 // aux1 est de la forme suivante
 //  511                 383                 255                 127                 0  
 //  |-------------------|-------------------|-------------------|-------------------|
 //  |_______A0B1________|________A1B0_______|________A1B1_______|________A0B0_______|
 //                    384                 256                 128         
 
 // Calcul de A[0]B[1]+A[1]B[0]
 middle.vec128 = _mm_xor_si128(aux1.vec128[3],aux1.vec128[2]);
// On récupère facilement les 64 bits de poids faible de middle en utilisant la structure 
 // __reg128. Idem pour les 64 bits de poids fort de A0B0 en utilisant la structure __reg512.
 aux1.vec64[1]^=middle.vec64[0];
 C->vec128[0] = aux1.vec128[0];
 // idem pour additionner les 64 bits de poids fort de middle aux 64 bits de poids faible
 // de A1B1.
 aux1.vec64[2]^=middle.vec64[1];
 C->vec128[1] = aux1.vec128[1];        
 }
 
 int main(void)
 {
  __reg128 A, B;
  __reg256 *C = calloc(1,sizeof(__reg256));
  
   // on initialise les registres 128 bits
   // {64 bits poids fort, 64 bits de poids faible}
   A.vec128 = _mm_set_epi64x(0xfffabfffeeffffff,0xffffaa1256ee1234);
   B.vec128 = _mm_set_epi64x(0xbfeefffdffffffff,0xea0d362010800099);

   printf("Données : \n");
   printf("A0=0x%lx\n", A.vec64[0]);
   printf("A1=0x%lx\n", A.vec64[1]);
   printf("B0=0x%lx\n", B.vec64[0]);
   printf("B1=0x%lx\n\n", B.vec64[1]);
   karat_mult_F2(C, A, B);
   printf("Resultat : \n");
   print128("C0",64,(__reg128)(C->vec128[0]));
   print128("C1",64,(__reg128)(C->vec128[1]));
}