/*
   fichier gfsr32.c
  créé le jeudi 5 janvier 2006 par Yann GUIDON (whygee@f-cpu.org)
  version 29 janvier 2006

Ce fichier crée un "Registre à Décalage à Rétroaction Généralisée"
plus ou moins paramétrique fournissant des nombres pseudo-aléatoires
avec une très longue période et un encombrement modeste en mémoire.

Le LFSR principal contient 32 registres de N bits chacun,
avec comme trinôme générateur x^N + x + 1. Les 32 registres
sont brassés par un LFSR utilisant le polynôme de CRC32.
La période théorique de l'ensemble est de 2^(N*32 -5).

Au besoin, N (défini par TAILLE_GFSR) peut être fixé à 4, 6, 9, 15, 22,
28, 30, 46, 60, 63, 127, 153, 172, 303, 471, 532, 865, 900...
La période pour N=3 a été mesurée à 4294967296 (2^32) (pas bien)
La période pour N=4 est inconnue mais supérieure à 2,5*10^12 (intéressant)
La période pour N=15 est inconnue mais supérieure à 3*10^12.

Autotest :
  gcc -DTEST_GFSR -Os -march=pentium3 -fomit-frame-pointer -o gfsr32 gfsr32-15.c
(définir d'abord le test à exécuter)

avec N=15:
Iteration 3262647581826 : Match partiel de 1 nombres.
iteration 3266591064067 ...Command terminated by signal 2
31483.10user 3.36system 8:45:27elapsed 99%CPU
===> La vitesse de recherche est d'environ 10^8 mots par seconde @700MHz
*/

#define chrono

#define CRC32_standard  0x04C11DB7UL
#define TAILLE_GFSR     4

typedef unsigned int U32; /* DOIT être sur 32 bits */
typedef   signed int S32; /* DOIT être sur 32 bits */

U32 GFSR_table[TAILLE_GFSR], GFSR_temp, GFSR_index;

void init_GFSR(void) {
  int i=0;
  U32 j, t=0x1A2B3C4D;  /* seed */

  do {
    /* applique une petite couche d'uniformisation : */
    t ^= (t>>5) ^ (t<<1);

    /* calcule un pas de LFSR : */
    j=t>>31;
    t<<=1;
    if (j) /* on teste l'ancien MSB de t */
      t^=CRC32_standard;

    GFSR_table[i++]=t;
  } while (i<TAILLE_GFSR);

  GFSR_temp=GFSR_table[TAILLE_GFSR-1];
  GFSR_index=0;
}

/* gcc-3.3.4 -Os -march=pentium3 -fomit-frame-pointer : 95,22s / 2^32 iterations */
inline U32 random_GSFR (void) {
  S32 t;

  GFSR_temp^=GFSR_table[GFSR_index];     /* lecture et XOR */

  /* étape de LFSR */
  t=((signed int)GFSR_temp) >> 31;
  GFSR_temp<<=1; /* ou GFSR_temp+=GFSR_temp; si shifters limités */
  GFSR_temp^=t & CRC32_standard;

  GFSR_table[GFSR_index]=GFSR_temp;  /* écriture */

  if (++GFSR_index >= TAILLE_GFSR)
    GFSR_index = 0;

  return GFSR_temp;
}

/* gcc-3.3.4 -Os -march=pentium3 -fomit-frame-pointer : 50s / 2^32 iterations */
#define GFSR_MACRO(index)     {        \
  signed int t;                        \
  GFSR_temp^=GFSR_table[index];        \
  t=((signed int)GFSR_temp) >> 31;     \
  GFSR_temp<<=1;                       \
  GFSR_temp^=t & CRC32_standard;       \
  GFSR_table[index]=GFSR_temp;  }


/* Remplit un tampon donné avec de la bouillie sur 32 bits.
   Ici on prend le parti que tout doit être fait "en interne"
   pour éviter que le compilo ne fasse trop d'accès en mémoire,
   donc on maximise les inférences de valeurs. En particulier,
   on élimine la dernière assignation GFSR_table[index]=GFSR_temp
   et on effectue un renommage cyclique : */

#define GFSR_MACRO2(current, last) { \
  S32 t;                             \
  current^=last;                     \
  t=((signed int)current) >> 31;     \
  current<<=1;                       \
  current^=t & CRC32_standard;  }

/* size et ptr sont alignés sur 4 octets par simplicité
     33.86s / 2^32 mots @700MHz, ou  */
void random_buffer(U32 *ptr, int size) {
  U32 GFSR0, GFSR1, GFSR2, GFSR3;
  int index2, size2;

/* Chargement de la table, modulo l'index de départ : */
  index2 = (GFSR_index+1) & 3;
                                    /* +  00 01 10 11 */
  GFSR0 =GFSR_table[GFSR_index];    /* 00 00 01 10 11 */
  GFSR1 =GFSR_table[index2];        /* 01 01 10 11 00 */
  GFSR2 =GFSR_table[GFSR_index^2];  /* 10 10 11 00 01 */
  GFSR3 =GFSR_table[index2^2];      /* 11 11 00 01 10 */
/* nb: par définition, GFSR_temp=GFSR_table[GFSR_index+3],
   donc pas besoin d'y toucher */

  size2=size>>2;
  while (size2--) {
    GFSR_MACRO2(GFSR0, GFSR3)   ptr[0]=GFSR0;
    GFSR_MACRO2(GFSR1, GFSR0)   ptr[1]=GFSR1;
    GFSR_MACRO2(GFSR2, GFSR1)   ptr[2]=GFSR2;
    GFSR_MACRO2(GFSR3, GFSR2)   ptr[3]=GFSR3;
    ptr+=4;
  }

  /* Sauve la table */
            GFSR_table[GFSR_index]   =GFSR0;
            GFSR_table[index2]       =GFSR1;
            GFSR_table[GFSR_index^2] =GFSR2;
  GFSR_temp=GFSR_table[index2^2]     =GFSR3;

  size &= 3;
  while (size--) {
    *(ptr++) = random_GSFR();
  }
}

/***************Auto-test******************/
#ifdef TEST_GFSR

#include <stdio.h> /* pour printf() */

#define TEST_REMPLISSAGE

unsigned long long int iteration; /* 64 bits */

#ifdef test_pattern2
#define TAILLE_SEQUENCE (10)
  U32 sequence[TAILLE_SEQUENCE], debut, j;
#endif

#ifdef TEST_REMPLISSAGE
U32 BUFFER[4096];
#endif

char binaire[]="%02d                                 \n";

int affiche_binaire(U32 n) {
  int i, j=0;

  for (i=31; i>=0; i--) {
    if (n & 1) {
      j++;
      binaire[i+5] = '*';
    }
    else
      binaire[i+5] = ' ';
    n>>=1;
  }

  printf(binaire, j);
  return j;
}

int affiche_decimal(U32 n) {
  int i, j=0, k=n;

  for (i=31; i>=0; i--) {
    if (n & 1)
      j++;
    n>>=1;
  }

  printf("%02d %8X\n", j, k);
  return j;
}

/* Affichage binaire de l'état du LFSR */
void affiche_GFSR(void) {
  int j=0, k;

  for (k=0; k<TAILLE_GFSR; k++)
    j+=affiche_binaire(GFSR_table[k]);

  printf("%d bits à 1 sur %d\n",j,32*TAILLE_GFSR);
}

int main(void) {

#ifdef TEST_REMPLISSAGE
  int i=1<<17;  /* 512MB= 2^29 - 12 = 17 = 128K */
  int j, k;
#endif

  if (sizeof(U32) != 4) {
    printf("Erreur : taille de U32 incorrecte\n");
    return -1;
  }

  init_GFSR(); /* initialise aussi GFSR_temp et GFSR_index */

#ifdef test_pattern1
  affiche_GFSR();
#endif

#ifdef test_pattern2
  /* enregistre le début de la séquence */
  debut=random_GSFR();
  for (iteration=1; iteration<TAILLE_SEQUENCE; iteration++)
    sequence[iteration]=random_GSFR();

  do {
    /* séquence d'alignement */
    while (GFSR_index!=0) {
      if (random_GSFR()==debut)
        goto match;
    }
    /* ici, GFSR_index=0; */

    while(1) { /* c'est la boucle qui consomme 99,9999% du CPU */
      GFSR_MACRO(0)
      if (debut==GFSR_temp) goto match0;
      GFSR_MACRO(1)
      if (debut==GFSR_temp) goto match1;
      GFSR_MACRO(2)
      if (debut==GFSR_temp) goto match2;
      GFSR_MACRO(3)
      if (debut==GFSR_temp) goto match3;
#ifdef loop15
      GFSR_MACRO(4)
      if (debut==GFSR_temp) goto match4;
      GFSR_MACRO(5)
      if (debut==GFSR_temp) goto match5;
      GFSR_MACRO(6)
      if (debut==GFSR_temp) goto match6;
      GFSR_MACRO(7)
      if (debut==GFSR_temp) goto match7;
      GFSR_MACRO(8)
      if (debut==GFSR_temp) goto match8;
      GFSR_MACRO(9)
      if (debut==GFSR_temp) goto match9;
      GFSR_MACRO(10)
      if (debut==GFSR_temp) goto match10;
      GFSR_MACRO(11)
      if (debut==GFSR_temp) goto match11;
      GFSR_MACRO(12)
      if (debut==GFSR_temp) goto match12;
      GFSR_MACRO(13)
      if (debut==GFSR_temp) goto match13;
      GFSR_MACRO(14)
      if (debut==GFSR_temp) goto match14;
#endif

      iteration+=TAILLE_GFSR;

#ifndef chrono
      if ((iteration & 0xFFFFFF0)==0) {
        printf("%citeration %Ld ...", 0xd,iteration);
        fflush(NULL);
      }
#else
      if (iteration >> 32) {
        printf("\nfini !\n");
        return 0;
      }
#endif
    }

#ifdef loop15
    match14: iteration  ++, GFSR_index ++;
    match13: iteration  ++, GFSR_index ++;
    match12: iteration  ++, GFSR_index ++;
    match11: iteration  ++, GFSR_index ++;
    match10: iteration  ++, GFSR_index ++;
    match9:  iteration  ++, GFSR_index ++;
    match8:  iteration  ++, GFSR_index ++;
    match7:  iteration  ++, GFSR_index ++;
    match6:  iteration  ++, GFSR_index ++;
    match5:  iteration  ++, GFSR_index ++;
    match4:  iteration  ++, GFSR_index ++;
#endif
    match3:  iteration  ++, GFSR_index ++;
    match2:  iteration  ++, GFSR_index ++;
    match1:  iteration  ++, GFSR_index ++;
    match0:  iteration  ++, GFSR_index ++;

    if (GFSR_index >= TAILLE_GFSR)
      GFSR_index -=TAILLE_GFSR;
    
match:
    printf ("%cIteration %Ld : ",0xd,iteration);
    j=0;
    do {
      if (++j==TAILLE_SEQUENCE) {
        printf ("%cBoucle trouvée !\n%c",7,7); /* réveiller l'opérateur*/
        return 0;
      }
      iteration++;
    } while (random_GSFR()==sequence[j]);

    printf ("Match partiel de %d nombres.\n", j);
  } while (1);
#endif /* test_pattern2 */

#ifdef TEST_REMPLISSAGE
#ifdef chrono
  do {
    random_buffer(BUFFER,1024);
  } while (--i);
#else
/* vérifie que le résultat de la routine optimisée est le même que
   celui de la routine de référence. points sensibles : routines
   d'alignement et jointure entre deux appels. */
  k=0;
  for (j=1; j<60; j++) {
    init_GFSR(); /* réinitialise tout */
    random_buffer(BUFFER,j);
    BUFFER[j]=random_GSFR();
    random_buffer(BUFFER+j+1,j);
    init_GFSR(); /* re-réinitialise tout */
    for (i=0; i<(j*2)+1; i++) {
      k+=affiche_binaire(BUFFER[i] ^random_GSFR());
      /* le XOR de deux données identiques devrait donner
         0 partout (ça permet de coder plus vite le test) */
    }
    printf ("------------------------------------\n");
  }
  printf (" =%d bits de différence\n",k);
#endif /* chrono */
#endif /* TEST_REMPLISSAGE */

  return 0;
}

#endif