/*
   fichier crc16-64.c
   version du 16 janvier 2007 : remise en forme
   14 mars : version temporaire pour diffusion
*/

#include <string.h>

#include "def64.h" /* pour utiliser les intrinsèques de GCC */

/* Directives de vorbosage */
/*
#define PRINTCRC
#define PRINTCRC2
#define TESTS_TABLE
*/

#ifdef PRINTCRC
  #include <stdio.h>
  #define PRINT_ULL(x) printf("%016llX ", x )
#else
  #define PRINT_ULL(x) {}
#endif

/* Les registres parallèles de CRC sont
   situés en mémoire globale pour simplifier l'adressage.
   Attention, ce n'est donc ni réentrant, ni thread-safe */
U64 CRC16_64_table[16];

volatile U64 moins_un= (U64)-1LL; /* utile pour les masques, rapide pour x86 */

/* Constantes arbitraire pour l'initialisation des registres de CRC
  + 7 rotations de la table à 16 entrées en 64 bits, soit 1024 octets
  seules les premières entrées sont initialisées : */
U64 CRC16_64_const[16*8]={
#ifdef TESTS_TABLE /* c'est plus facile à lire pour vérifier l'algorithme */
  0x0102030405060708ULL, 0x1112131415161718ULL,
  0x2122232425262728ULL, 0x3132333435363738ULL,
  0x4142434445464748ULL, 0x5152535455565758ULL,
  0x6162636465666768ULL, 0x7172737475767778ULL,
  0x8182838485868788ULL, 0x9192939495969798ULL,
  0xA1A2A3A4A5A6A7A8ULL, 0xB1B2B3B4B5B6B7B8ULL,
  0xC1C2C3C4C5C6C7C8ULL, 0xD1D2D3D4D5D6D7D8ULL,
  0xE1E2E3E4E5E6E7E8ULL, 0xF1F2F3F4F5F6F7F8ULL
#else
  0x62f08f0a115bac5bULL, 0xdb93dbb089865d38ULL,
  0x81b68d1814441139ULL, 0x912a63004ad95b66ULL,
  0x2622a84506da1598ULL, 0x5d1bd29770306f74ULL,
  0xf4c7f1d32f046ca0ULL, 0x8b4f965bdf10ab50ULL,
  0x29ab650f5ed5446dULL, 0xef8cf77dc5632827ULL,
  0x23e71d7d2afa86dcULL, 0xd6fcf6213866c513ULL,
  0xfb59fc246bd3b6f3ULL, 0x4171947756a1e9afULL,
  0x36465e621bf4308dULL, 0x01451507db3dc05dULL
#endif
/* la suite est complétée par rotation_table() */
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
};

/* ------------- version alignée de base ------------- */

/*
 aligned_CRC16_64 : fonction de référence du calcul de CRC16-64.
       *pointeur est limité à des mots alignés sur 64 bits.
 paramètres : pointeur du début du bloc et taille du bloc en octets.
 */
U16 aligned_CRC16_64(U32 taille_du_bloc, U64 *pointeur) {
  int taille = taille_du_bloc >> 3;
  int index=0; /* modulo 16 */
  U64 t, u, mask; /* pour l'épilogue */

  memcpy(CRC16_64_table,CRC16_64_const,sizeof(CRC16_64_table));

  while (taille--) {
    CRC16_64_table[index] = U64_XOR(
      U64_XOR(CRC16_64_table[index], *(pointeur++)),
      U64_XOR(CRC16_64_table[(index+1)&15],
              CRC16_64_table[(index+14)&15] ));
    PRINT_ULL(CRC16_64_table[index]);
    index=(index+1) & 15; /* rotation des registres */
  }

/* épilogue */
  taille = taille_du_bloc & 7;
  if (taille) {

    #ifdef PRINTCRC
      printf(" (%d octets restant) ", taille);
    #endif

    t = CRC16_64_table[index];
    u = U64_XOR(
      U64_XOR(t , *pointeur),
      U64_XOR(CRC16_64_table[(index+1)&15],
              CRC16_64_table[(index+14)&15] ));

    /* construction du masque */
    mask = moins_un;
    taille <<= 3;
    mask = LOCAL_SHLr(mask, taille);  /* attention : LittleEndian*/

    /* sélectionne les bons octets */
    t    = U64_AND(t,mask);
    mask = U64_ANDN(mask, u);
    t    = U64_OR(mask, t) ;

    CRC16_64_table[index] = t;
    PRINT_ULL(t);
  }

  return 0; /* pas de compaction ici */
}


/* ------------- version alignée plus rapide ------------- */

/*
 aligned_CRC16_64_2 : fonction de CRC16-64 accélérée, déroulée.
   Elle travaille sur des blocs de 16 mots de 64 bits
   (128 octets à la fois par itération), pointeur aligné sur 8 octets.
 paramètres : pointeur du début du bloc et taille du bloc en octets.
*/
U16 aligned_CRC16_64_2(U32 taille_du_bloc, U64 *p) {
  int taille = taille_du_bloc >> 7;
  U64 r0, r1, r2, r3;

  /* initialisation */
  memcpy(CRC16_64_table,CRC16_64_const,sizeof(CRC16_64_table));

  r0 = CRC16_64_table[0];
/*r1 = CRC16_64_table[1]; effectué lors de la première itération */
  r2 = CRC16_64_table[14];
  r3 = CRC16_64_table[15];

  #define CRC_ROUND(x, a, b, c) \
      r##b = CRC16_64_table[(x+1) & 15]; \
      r##a = U64_XOR( U64_XOR( r##a ,p[x]), U64_XOR( r##b , r##c )); \
      CRC16_64_table[(x) & 15] = r##a ; \
      PRINT_ULL(r##a);
  
  while (taille--) { /*
 l'instanciation : a = (x) % 4
                   |  b = (x+1) % 4
                   |  |  c = (x+14) % 4  */
    CRC_ROUND( 0, 0, 1, 2)
    CRC_ROUND( 1, 1, 2, 3)
    CRC_ROUND( 2, 2, 3, 0)
    CRC_ROUND( 3, 3, 0, 1)
    CRC_ROUND( 4, 0, 1, 2)
    CRC_ROUND( 5, 1, 2, 3)
    CRC_ROUND( 6, 2, 3, 0)
    CRC_ROUND( 7, 3, 0, 1)
    CRC_ROUND( 8, 0, 1, 2)
    CRC_ROUND( 9, 1, 2, 3)
    CRC_ROUND(10, 2, 3, 0)
    CRC_ROUND(11, 3, 0, 1)
    CRC_ROUND(12, 0, 1, 2)
    CRC_ROUND(13, 1, 2, 3)
    CRC_ROUND(14, 2, 3, 0)
    CRC_ROUND(15, 3, 0, 1)
    p+=16;
  }

  return 0; /* pas encore de compaction */
}

#undef CRC_ROUND

/* ------------- version non alignée ------------- */

/* Précalcul des rotations de la table */
void rotation_table(void) {
  #if defined (__MMX__) /* version 64 bits forcée */
      register U64 t asm("mm0");
      register U64 u asm("mm1");
      register U64 v asm("mm2");
  #else
    U64 t, u, v;
  #endif
  int i=8, i2=56, /* à la fois compteurs et opérandes de décalage */
      j, /* compteur et index dans la table source */
      k=16; /* index dans la table de destination */

  /* recopie chaque entrée et effectue la rotation dans le bon sens */
  do {
    j=0; do {
      t=CRC16_64_const[j];         /* entrée originale */
      u=CRC16_64_const[(j-1)&15];  /* entrée décalée d'un cran */
      j++;
     /* attention : pas de correction du sens du shift ici ! */
      v=U64_OR( U64_SHLr(t,i), U64_SHRr(u,i2));
      CRC16_64_const[k]= CONVERT_LE64(v);
      k++;
    } while (j<16);
    i+=8, i2-=8;
  } while (i2!=0);

  /* Attention : entrée initiale = pas de rotation
      à faire sauf si Big Endian mais en dernier */
  #ifdef WORDS_BIGENDIAN
  for (k=0; k<16; k++)
    CRC16_64_const[k] = CONVERT_LE64(CRC16_64_const[k]);
  #endif

  #ifdef PRINTCRC2
    for (k=0; k<8*16; k+=16)
      print_table(&(CRC16_64_const[k]));
  #endif
}


void unaligned_CRC16_64(U32 taille_du_bloc, U8 *q) {
/* attention : changement de la taille des accès, de U64 à U8,
  car la granularité d'accès est maintenant d'un octet
  (sinon ce con de gcc vire les 3 LSB pour aligner) */

/* déclarations : utilise 6 entiers 64 bits (+2 temporaires) */
  U64 r0, r1, r2, r3; /* variables pour le calcul du CRC */
  U64 mask, t; /* pour construire les masques */

/* 4 variables + 2 paramètres d'appel
 = 6 registres 32 bits (quasi-limite du x86)  */
  unsigned offset = PTR_CAST q & 7; /* calcule l'alignement */
    signed rem;  /* nombre d'octets restant dans le mot */
  U64 *ptr_table = CRC16_64_const + (offset<<4);
  int index; /* peut utiliser le même registre que ptr_table */

  /* Recopie la bonne partie de la table des constantes
  memcpy(CRC16_64_table,CRC16_64_table_rotations+(offset<<7),128);
    problème : memcpy est "blocant" donc on "déroule"
    et "entrelace" (saupoudre) les instructions.
  Lit en premier les données qui nous intéressent d'abord : */
    r0 = ptr_table[0];  CRC16_64_table[0] = r0;
    r1 = ptr_table[1];  CRC16_64_table[1] = r1;

  /* puis remplit le pipeline avec plein d'instructions,
  pour "masquer" la latence d'exécution des autres instructions.
  On compte sur le renommage des registres par le CPU
  pour compenser le codage sub-optimal ci-dessous : */

  /* reste de l'initialisation */
  mask=moins_un;
    r2 = ptr_table[14];  CRC16_64_table[14] = r2;

  rem = taille_du_bloc - 8;  /* décompte un mot entier */
    r3 = ptr_table[15];  CRC16_64_table[15] = r3;

/* Problème lors de la copie : lorsque les données ne
   sont pas touchées ou réutilisées dans un calcul,
   gcc3.3.4 passe stupidement par une paire de registres
   32 bits, au lieu d'utiliser un registre MMX pour déplacer
   les valeurs 64 bits. Une petite macro fait l'affaire.
   Et pour booster, une version SSE (128bits) est fournie : */

#if defined (__SSE__)  /* version 128 bits */
  #define MOVTABLE(index) {                                      \
     __builtin_ia32_storeups ((float*)(CRC16_64_table+index) ,   \
      __builtin_ia32_loadups ((float*)     (ptr_table+index)));  }
#else
  #if defined (__MMX__) /* version 64 bits forcée */
    /* aucune intrinsèque pour movq disponible */
    #define MOVTABLE(index) {      \
      register U64 y asm("mm0");   \
      register U64 z asm("mm1");   \
      y = ptr_table[index];        \
      z = ptr_table[index+1];      \
      CRC16_64_table[index]   = y; \
      CRC16_64_table[index+1] = z; }
  #else /* version normale */
    #define MOVTABLE(index) {                       \
      CRC16_64_table[index]   = ptr_table[index];   \
      CRC16_64_table[index+1] = ptr_table[index+1]; }
  #endif
#endif

    MOVTABLE(2);

  /* efface les 3 LSB, donc aligne le pointeur sur 8 octets : */
  q = (U8*) ((PTR_CAST q) & -8L);
    MOVTABLE(4);

  /* création du masque pour traiter le premier mot */
  mask = LOCAL_SHLr(mask, (offset<<3));
    MOVTABLE(6);

  /* lecture des 8 premiers octets (alignés) */
  t = *(U64*)q;
  q+=8;
    MOVTABLE(8);

  /* calcule un pas de CRC avec t comme donnée */
  t  = U64_XOR( U64_XOR(r0,t), U64_XOR(r1,r2));
    MOVTABLE(10);

  /* Gestion du cas où trop peu d'octets sont à traiter */
  rem += offset; /* compense avec l'offset */
    MOVTABLE(12);

  if (rem <= 0) {  /* version court-circuit */

    /* Il faut encore tronquer le masque. Deux décalages
    en sens opposés sont plus simples qu'un masque supplémentaire. */
    rem = ((-rem) << 3);
    #if defined (__MMX__) /* version MMX optimisée */
    {
      U64 tmp;
      U64_MOVD(tmp, rem)
      mask = U64_SHR(U64_SHL(mask, tmp), tmp);
    }
    #else
      mask = LOCAL_SHR(LOCAL_SHL(mask, rem), rem);
    #endif

    /* réécrit le résultat, mais uniquement les octets choisis */
    t    = U64_AND(t,mask);
    mask = U64_ANDN(mask, r0);
    r0   = U64_OR(mask, t) ;
    CRC16_64_table[0] = r0;
    PRINT_ULL(r0);

    /* à ce moment, on pourrait tenter une "compaction rapide"
    avec des valeurs précalculées plus petites... plus tard. */
  }

  /* cas normal où plus d'un mot doit être traité */
  else {
    /* réécrit le résultat, mais uniquement les octets choisis */
    /* note : ce code est une copie du précédent bloc,
          pour réduire le nombre de branchements */
    t    = U64_AND(t,mask);
    mask = U64_ANDN(mask, r0);
    r0   = U64_OR(mask, t) ;
    CRC16_64_table[0] = r0;
    PRINT_ULL(r0);

    #define CRC_ROUND_U8(x, a, b, c) \
       r##b = CRC16_64_table[(x+1) & 15]; \
       r##a = U64_XOR( U64_XOR( r##a , ((U64*)q)[x-1]), \
                       U64_XOR( r##b , r##c )); \
       CRC16_64_table[(x) & 15] = r##a ; \
       PRINT_ULL(r##a);

    /* la version ultra-rapide */
    while (rem >= 128) {
   /* l'instanciation : a = (x) % 4
                        |  b = (x+1) % 4
                        |  |  c = (x+14) % 4  */
      CRC_ROUND_U8( 1, 1, 2, 3)
      CRC_ROUND_U8( 2, 2, 3, 0)
      CRC_ROUND_U8( 3, 3, 0, 1)
      CRC_ROUND_U8( 4, 0, 1, 2)
      CRC_ROUND_U8( 5, 1, 2, 3)
      CRC_ROUND_U8( 6, 2, 3, 0)
      CRC_ROUND_U8( 7, 3, 0, 1)
      CRC_ROUND_U8( 8, 0, 1, 2)
      CRC_ROUND_U8( 9, 1, 2, 3)
      CRC_ROUND_U8(10, 2, 3, 0)
      CRC_ROUND_U8(11, 3, 0, 1)
      CRC_ROUND_U8(12, 0, 1, 2)
      CRC_ROUND_U8(13, 1, 2, 3)
      CRC_ROUND_U8(14, 2, 3, 0)
      CRC_ROUND_U8(15, 3, 0, 1)
      CRC_ROUND_U8(16, 0, 1, 2)
 /* note : c'est décalé par rapport à la version
 initiale car la première ronde est déjà calculée
 mais ça décale l'accès à la mémoire, d'où le (x-1)
 dans la macro */
      rem -= 128;
      q += 128;
    }

    index=1;
    /* la version... juste rapide */
    while (rem >= 8) {
      CRC16_64_table[index] = U64_XOR(
        U64_XOR(CRC16_64_table[index], *(U64*)q),
        U64_XOR(CRC16_64_table[(index+1)&15],
                CRC16_64_table[(index+14)&15]));
      PRINT_ULL(CRC16_64_table[index]);
      index=(index+1) & 15; /* rotation des registres */
      rem -= 8;
      q += 8;
    }

/* finalement, on traite le petit bout qui dépasse à la fin : */
    if (rem != 0) {
#ifdef PRINTCRC
printf(" (%d octets restant) ", rem);
#endif
      r0 = CRC16_64_table[index];

      /* dernier calcul */
      t = U64_XOR(
          U64_XOR(r0 , *(U64*)q),
          U64_XOR(CRC16_64_table[(index+1)&15],
                  CRC16_64_table[(index+14)&15]));

      /* construction du nouveau masque */
      mask = moins_un;
      rem = (rem << 3);
      mask = LOCAL_SHLr(mask, rem); /* Attention à l'endian */

      /* sélectionne les bons octets */
      r0   = U64_AND(r0,mask);
      mask = U64_ANDN(mask, t);
      r0   = U64_OR(mask, r0) ;
      CRC16_64_table[index] = r0;
      PRINT_ULL(r0);
    }
  }
}

#undef MOVTABLE
#undef CRC_ROUND_U8

#undef PRINT_ULL
