/*
 * Decouverte des episodes frequents
 *
 * A compiler par :  gcc -O4 e5.s frequent.c -lpcre
 *
 */

#include <string.h>
#include <stdio.h>
#include <pcre.h>

#define pattern  "(\\w+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)|(\\d+)"
#define myfile   "/data3/christophe/DOC/MICHEL/data/BigTable.dat"
/* On donne ici le nombre de colonnes du fichier texte ci dessus */
#define COL 11

#define SUPPORT 68567		/* Ca represente 10% du nombre de lignes de la table */
#define VEC_SIZE 72		/* should be a multiple of 3 */

/* Bit-twiddling lexicographic combination generator
   Doug Moore (dougm@rice.edu) (C) 1993 */

/* The procedures implemented here provide a method for enumerating
   combinations, expressed as bitmasks, of n items taken k at a time in
   lexicographic order.  The types, procedures and their purposes are: */

typedef unsigned long long bitmask_t;
/* represents a subset of the integers ({0, .., 63} (on most machines)
   E.g., bitmask = 00110 (in binary) represents set {1,2} */

bitmask_t lexFirstComb (unsigned k);
/* Returns the lexicographically first combination of k items (i.e. {0,1,..,k-1}) */

bitmask_t leastItem (bitmask_t comb);
/* Returns the least item in a combination (i. e. leastItem({2,4,5}) == {2} */

bitmask_t lexLastComb (unsigned n, unsigned k);
/* Returns the lexicographically last combination of k items from n
   (i.e. {n-k,..,n-1}) */


#ifdef COMPUTE_PARITY
#define PARITY_ARG , int* parity
#define PARITY_PARAM , &parity
#else
#define PARITY_ARG
#define PARITY_PARAM
#endif

bitmask_t lexNextComb (bitmask_t comb PARITY_ARG);
/* Returns the lexicographically next combination after input comb; updates
   parity to value of parity of returned combination */


#include <assert.h>

static const bitmask_t evenBits = ~0ULL / 3;	/* == ...0101010101 in binary */

#define Lshift1(i) ((bitmask_t)1 << (i)/2 << ((i)-(i)/2))

bitmask_t
lexFirstComb (unsigned k)
{
  assert (k <= 8 * sizeof (bitmask_t));
  return Lshift1 (k) - 1;
}

bitmask_t
leastItem (bitmask_t comb)
{
  return comb & -comb;
}

bitmask_t
lexLastComb (unsigned n, unsigned k)
{
  assert (k <= n);
  return Lshift1 (n) - Lshift1 (n - k);
}

bitmask_t
lexNextComb (bitmask_t comb PARITY_ARG)
{
  bitmask_t hibit;
  bitmask_t lobit = leastItem (comb);
  comb += lobit;
  hibit = leastItem (comb);
#ifdef COMPUTE_PARITY
  *parity ^= ((hibit | lobit) & evenBits) != 0;
#endif
  hibit -= lobit;
  while (!(hibit & 1))
    hibit >>= 1;
  comb |= hibit >> 1;
  return comb;
}

bitmask_t
lexPrevComb (bitmask_t comb)
{
  bitmask_t lobits = leastItem (~comb) - 1;
  comb &= ~lobits;
  bitmask_t hibit = leastItem (comb);
  while (lobits)
    lobits >>= 1, hibit >>= 1;
  comb -= hibit >> 1;
  return comb;
}

/* Les colonnes qui nous interessent sont (de 1 a 11) */
/* 3 : CPU usage        */
/* 5 : memory Available */
/* 7 : Nb Processus */
/* 10: IP datagrams */
/* On a donc 4 classes ; on divise ces 4 classes en 10 catégories */
/* Au total, on se retrouve avec 40 items. Il y a potentiellement 2^40 */
/* = 1.099.511.627.776 frequent episodes ! */
/* A chacun de ces 40 items est associé un bitset */
/* Comme la table /home/christophe/DOC/MICHEL/data/BigTable.dat a */
/* 685673 lignes, on prend des tailles de bitsets de 524288*2 bits, soit */
/* 65536*2 bytes. Comme on en a 40, on obtient une occupation memoire de */
/* 20971520*2 bytes... ce qui correspond a peu pret a la taille de la table */
/* On peut ne jamais travailler sur disque : apres la */
/* premiere phase de l'algorithme de calcul des supports, il va y avoir */
/* elimination de bitsets : on decide alors de reutiliser l'emplacement */
/* memoire des items definitivement eliminés */

/*
 * L'objet strucBitset contient un champ Nom, un champs NbElements et
 * un bitset de 131072 bytes
 */

struct structBitset
{
  int NbElements;		/* Nombre d'éléments dans le bitset : le support */
  unsigned char bitset[131072];
  int i, j, k, l, m, n, o, p;	/* pour retrouver qui est qui */
};

/*
 * On fabrique maintenant un tableau de 40 structures structBitset
 * qui sont nos 40 items
 */

struct structBitset TabBitset[40];
unsigned char MyRes[131072];
struct structBitset TabBitsetCopy[20];
struct structBitset TabBitsetOne[20];


int
main (int argc, char *argv[])
{
  pcre *re;
  int rc;
  const char *err_msg;
  int err_loc;
  int ovector[VEC_SIZE];
  int options;
  int offset;
  int i, j, l, p, r, colonne, myint, myligne, max, max1, max2;
  int verif, ind;
  long long TOTAL;
  /**** Pour lageneraton des combinaisons ************/
  unsigned n, k;
  bitmask_t lastComb, prev;
  bitmask_t comb;
  int parity = 0;
  /***************************************************/
  int done;
  FILE *f;
  int Candidate[40];		/* les k-itemset sont indicés par ce tableau */
  char mybuf[1024];
  char mycopy[1024];

  int m, MyBitset[2];

  re = pcre_compile (pattern, 0, &err_msg, &err_loc, NULL);
  if (re == NULL)
    {
      //    printf("regex failed: %s xat %d\n", err_msg, err_loc);
      printf ("PCRE compilation failed at offset %d: %s\n", err_loc, err_msg);
      exit (1);
    }

  done = 0;
  myligne = 0;

  for (i = 0; i < 40; i++)
    {
      TabBitset[i].NbElements = 0;
    }
  for (i = 0; i < 20; i++)
    {
      TabBitsetCopy[i].NbElements = 0;
      TabBitsetOne[i].NbElements = 0;
    }

  f = fopen (myfile, "r");
  while (!feof (f))
    {
      /* On lit une ligne du fichier texte */
      fgets (mybuf, 1024, f);
      rc =
	pcre_exec (re, NULL, mybuf, (int) strlen (mybuf), 0, 0, ovector,
		   VEC_SIZE);
      /* Si on trouve une ligne blanche alors on suppose que c'est fini */
      if (mybuf[0] == '\n')
	break;

      colonne = 0;
      done = 0;
      do
	{
	  /* on examine les colonnes 3, 5, 7, 10 */
	  if (!done
	      && ((colonne == 2 * 2) || (colonne == 4 * 2)
		  || (colonne == 6 * 2) || (colonne == 9 * 2)))
	    {
	      strncpy (mycopy, mybuf + ovector[0], ovector[1] - ovector[0]);
	      strcpy (mycopy + ovector[1] - ovector[0], "\0");
	      myint = atoi (mycopy);
	      // printf("%d: %d\n",myligne,myint);
	      switch (colonne)
		{
		case 4:	/* CPU */
		  if (myint < 10)
		    {
		      TabBitset[0].NbElements++;
		      SetBitTo1 (TabBitset[0].bitset, myligne);
		    }
		  else if (myint < 20)
		    {
		      TabBitset[1].NbElements++;
		      //    printf("%d: %d\n",myligne,TabBitset[1].NbElements);

		      SetBitTo1 (TabBitset[1].bitset, myligne);
		    }
		  else if (myint < 30)
		    {
		      TabBitset[2].NbElements++;
		      SetBitTo1 (TabBitset[2].bitset, myligne);
		    }
		  else if (myint < 40)
		    {
		      TabBitset[3].NbElements++;
		      SetBitTo1 (TabBitset[3].bitset, myligne);
		    }
		  else if (myint < 50)
		    {
		      TabBitset[4].NbElements++;
		      SetBitTo1 (TabBitset[4].bitset, myligne);
		    }
		  else if (myint < 60)
		    {
		      TabBitset[5].NbElements++;
		      SetBitTo1 (TabBitset[5].bitset, myligne);
		    }
		  else if (myint < 70)
		    {
		      TabBitset[6].NbElements++;
		      SetBitTo1 (TabBitset[6].bitset, myligne);
		    }
		  else if (myint < 80)
		    {
		      TabBitset[7].NbElements++;
		      SetBitTo1 (TabBitset[7].bitset, myligne);
		    }
		  else if (myint < 90)
		    {
		      TabBitset[8].NbElements++;
		      SetBitTo1 (TabBitset[8].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[9].NbElements++;
		      SetBitTo1 (TabBitset[9].bitset, myligne);
		    }
		  break;
		case 8:	/* MEM available */
		  if (myint < 1000)
		    {
		      TabBitset[10].NbElements++;
		      SetBitTo1 (TabBitset[10].bitset, myligne);
		    }
		  else if (myint < 8000)
		    {
		      TabBitset[11].NbElements++;
		      SetBitTo1 (TabBitset[11].bitset, myligne);
		    }
		  else if (myint < 10000)
		    {
		      TabBitset[12].NbElements++;
		      SetBitTo1 (TabBitset[12].bitset, myligne);
		    }
		  else if (myint < 20000)
		    {
		      TabBitset[13].NbElements++;
		      SetBitTo1 (TabBitset[13].bitset, myligne);
		    }
		  else if (myint < 40000)
		    {
		      TabBitset[14].NbElements++;
		      SetBitTo1 (TabBitset[14].bitset, myligne);
		    }
		  else if (myint < 60000)
		    {
		      TabBitset[15].NbElements++;
		      SetBitTo1 (TabBitset[15].bitset, myligne);
		    }
		  else if (myint < 80000)
		    {
		      TabBitset[16].NbElements++;
		      SetBitTo1 (TabBitset[16].bitset, myligne);
		    }
		  else if (myint < 100000)
		    {
		      TabBitset[17].NbElements++;
		      SetBitTo1 (TabBitset[17].bitset, myligne);
		    }
		  else if (myint < 300000)
		    {
		      TabBitset[18].NbElements++;
		      SetBitTo1 (TabBitset[18].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[19].NbElements++;
		      SetBitTo1 (TabBitset[19].bitset, myligne);
		    }
		  break;
		case 12:	/* N processes */
		  if (myint < 10)
		    {
		      TabBitset[20].NbElements++;
		      SetBitTo1 (TabBitset[20].bitset, myligne);
		    }
		  else if (myint < 20)
		    {
		      TabBitset[21].NbElements++;
		      SetBitTo1 (TabBitset[21].bitset, myligne);
		    }
		  else if (myint < 30)
		    {
		      TabBitset[22].NbElements++;
		      SetBitTo1 (TabBitset[22].bitset, myligne);
		    }
		  else if (myint < 40)
		    {
		      TabBitset[23].NbElements++;
		      SetBitTo1 (TabBitset[23].bitset, myligne);
		    }
		  else if (myint < 50)
		    {
		      TabBitset[24].NbElements++;
		      SetBitTo1 (TabBitset[24].bitset, myligne);
		    }
		  else if (myint < 60)
		    {
		      TabBitset[25].NbElements++;
		      SetBitTo1 (TabBitset[25].bitset, myligne);
		    }
		  else if (myint < 70)
		    {
		      TabBitset[26].NbElements++;
		      SetBitTo1 (TabBitset[26].bitset, myligne);
		    }
		  else if (myint < 80)
		    {
		      TabBitset[27].NbElements++;
		      SetBitTo1 (TabBitset[27].bitset, myligne);
		    }
		  else if (myint < 90)
		    {
		      TabBitset[28].NbElements++;
		      SetBitTo1 (TabBitset[28].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[29].NbElements++;
		      SetBitTo1 (TabBitset[29].bitset, myligne);
		    }
		  break;
		case 18:	/* IP datagram */
		  if (myint < 1000)
		    {
		      TabBitset[30].NbElements++;
		      SetBitTo1 (TabBitset[30].bitset, myligne);
		    }
		  else if (myint < 3000)
		    {
		      TabBitset[31].NbElements++;
		      SetBitTo1 (TabBitset[31].bitset, myligne);
		    }
		  else if (myint < 5000)
		    {
		      TabBitset[32].NbElements++;
		      SetBitTo1 (TabBitset[32].bitset, myligne);
		    }
		  else if (myint < 7000)
		    {
		      TabBitset[33].NbElements++;
		      SetBitTo1 (TabBitset[33].bitset, myligne);
		    }
		  else if (myint < 9000)
		    {
		      TabBitset[34].NbElements++;
		      SetBitTo1 (TabBitset[34].bitset, myligne);
		    }
		  else if (myint < 11000)
		    {
		      TabBitset[35].NbElements++;
		      SetBitTo1 (TabBitset[35].bitset, myligne);
		    }
		  else if (myint < 14000)
		    {
		      TabBitset[36].NbElements++;
		      SetBitTo1 (TabBitset[36].bitset, myligne);
		    }
		  else if (myint < 16000)
		    {
		      TabBitset[37].NbElements++;
		      SetBitTo1 (TabBitset[37].bitset, myligne);
		    }
		  else if (myint < 18000)
		    {
		      TabBitset[38].NbElements++;
		      SetBitTo1 (TabBitset[38].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[39].NbElements++;
		      SetBitTo1 (TabBitset[39].bitset, myligne);
		    }
		  break;
		}
	    }
	  colonne++;
	  options = 0;
	  offset = ovector[1];
	  if (ovector[0] == ovector[1])
	    {
	      if (ovector[0] == strlen (mybuf))
		done = 1;
	      else
		{
		  options = PCRE_NOTEMPTY | PCRE_ANCHORED;
		}
	    }
	  rc =
	    pcre_exec (re, NULL, mybuf, strlen (mybuf), offset, options,
		       ovector, VEC_SIZE);
	  if (!done && rc == PCRE_ERROR_NOMATCH)
	    {
	      if (options == 0)
		done = 1;
	      else
		ovector[1] = offset + 1;
	    }
	  /* on examine les colonnes 3, 5, 7, 10 */
	  if (!done
	      && ((colonne == 2 * 2) || (colonne == 4 * 2)
		  || (colonne == 6 * 2) || (colonne == 9 * 2)))
	    {
	      strncpy (mycopy, mybuf + ovector[0], ovector[1] - ovector[0]);
	      strcpy (mycopy + ovector[1] - ovector[0], "\0");
	      myint = atoi (mycopy);
	      //printf("%d: %d\n",myligne,myint);
	      switch (colonne)
		{
		case 4:	/* CPU */
		  if (myint < 10)
		    {
		      TabBitset[0].NbElements++;
		      SetBitTo1 (TabBitset[0].bitset, myligne);
		    }
		  else if (myint < 20)
		    {
		      TabBitset[1].NbElements++;
		      SetBitTo1 (TabBitset[1].bitset, myligne);
		    }
		  else if (myint < 30)
		    {
		      TabBitset[2].NbElements++;
		      SetBitTo1 (TabBitset[2].bitset, myligne);
		    }
		  else if (myint < 40)
		    {
		      TabBitset[3].NbElements++;
		      SetBitTo1 (TabBitset[3].bitset, myligne);
		    }
		  else if (myint < 50)
		    {
		      TabBitset[4].NbElements++;
		      SetBitTo1 (TabBitset[4].bitset, myligne);
		    }
		  else if (myint < 60)
		    {
		      TabBitset[5].NbElements++;
		      SetBitTo1 (TabBitset[5].bitset, myligne);
		    }
		  else if (myint < 70)
		    {
		      TabBitset[6].NbElements++;
		      SetBitTo1 (TabBitset[6].bitset, myligne);
		    }
		  else if (myint < 80)
		    {
		      TabBitset[7].NbElements++;
		      SetBitTo1 (TabBitset[7].bitset, myligne);
		    }
		  else if (myint < 90)
		    {
		      TabBitset[8].NbElements++;
		      SetBitTo1 (TabBitset[8].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[9].NbElements++;
		      SetBitTo1 (TabBitset[9].bitset, myligne);
		    }
		  break;
		case 8:	/* MEM available */
		  if (myint < 1000)
		    {
		      TabBitset[10].NbElements++;
		      SetBitTo1 (TabBitset[10].bitset, myligne);
		    }
		  else if (myint < 8000)
		    {
		      TabBitset[11].NbElements++;
		      SetBitTo1 (TabBitset[11].bitset, myligne);
		    }
		  else if (myint < 10000)
		    {
		      TabBitset[12].NbElements++;
		      SetBitTo1 (TabBitset[12].bitset, myligne);
		    }
		  else if (myint < 20000)
		    {
		      TabBitset[13].NbElements++;
		      SetBitTo1 (TabBitset[13].bitset, myligne);
		    }
		  else if (myint < 40000)
		    {
		      TabBitset[14].NbElements++;
		      SetBitTo1 (TabBitset[14].bitset, myligne);
		    }
		  else if (myint < 60000)
		    {
		      TabBitset[15].NbElements++;
		      SetBitTo1 (TabBitset[15].bitset, myligne);
		    }
		  else if (myint < 80000)
		    {
		      TabBitset[16].NbElements++;
		      SetBitTo1 (TabBitset[16].bitset, myligne);
		    }
		  else if (myint < 100000)
		    {
		      TabBitset[17].NbElements++;
		      SetBitTo1 (TabBitset[17].bitset, myligne);
		    }
		  else if (myint < 300000)
		    {
		      TabBitset[18].NbElements++;
		      SetBitTo1 (TabBitset[18].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[19].NbElements++;
		      SetBitTo1 (TabBitset[19].bitset, myligne);
		    }
		  break;
		case 12:	/* N processes */
		  if (myint < 10)
		    {
		      TabBitset[20].NbElements++;
		      SetBitTo1 (TabBitset[20].bitset, myligne);
		    }
		  else if (myint < 20)
		    {
		      TabBitset[21].NbElements++;
		      SetBitTo1 (TabBitset[21].bitset, myligne);
		    }
		  else if (myint < 30)
		    {
		      TabBitset[22].NbElements++;
		      SetBitTo1 (TabBitset[22].bitset, myligne);
		    }
		  else if (myint < 40)
		    {
		      TabBitset[23].NbElements++;
		      SetBitTo1 (TabBitset[23].bitset, myligne);
		    }
		  else if (myint < 50)
		    {
		      TabBitset[24].NbElements++;
		      SetBitTo1 (TabBitset[24].bitset, myligne);
		    }
		  else if (myint < 60)
		    {
		      TabBitset[25].NbElements++;
		      SetBitTo1 (TabBitset[25].bitset, myligne);
		    }
		  else if (myint < 70)
		    {
		      TabBitset[26].NbElements++;
		      SetBitTo1 (TabBitset[26].bitset, myligne);
		    }
		  else if (myint < 80)
		    {
		      TabBitset[27].NbElements++;
		      SetBitTo1 (TabBitset[27].bitset, myligne);
		    }
		  else if (myint < 90)
		    {
		      TabBitset[28].NbElements++;
		      SetBitTo1 (TabBitset[28].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[29].NbElements++;
		      SetBitTo1 (TabBitset[29].bitset, myligne);
		    }
		  break;
		case 18:	/* IP datagram */
		  if (myint < 1000)
		    {
		      TabBitset[30].NbElements++;
		      SetBitTo1 (TabBitset[30].bitset, myligne);
		    }
		  else if (myint < 3000)
		    {
		      TabBitset[31].NbElements++;
		      SetBitTo1 (TabBitset[31].bitset, myligne);
		    }
		  else if (myint < 5000)
		    {
		      TabBitset[32].NbElements++;
		      SetBitTo1 (TabBitset[32].bitset, myligne);
		    }
		  else if (myint < 7000)
		    {
		      TabBitset[33].NbElements++;
		      SetBitTo1 (TabBitset[33].bitset, myligne);
		    }
		  else if (myint < 9000)
		    {
		      TabBitset[34].NbElements++;
		      SetBitTo1 (TabBitset[34].bitset, myligne);
		    }
		  else if (myint < 11000)
		    {
		      TabBitset[35].NbElements++;
		      SetBitTo1 (TabBitset[35].bitset, myligne);
		    }
		  else if (myint < 14000)
		    {
		      TabBitset[36].NbElements++;
		      SetBitTo1 (TabBitset[36].bitset, myligne);
		    }
		  else if (myint < 16000)
		    {
		      TabBitset[37].NbElements++;
		      SetBitTo1 (TabBitset[37].bitset, myligne);
		    }
		  else if (myint < 18000)
		    {
		      TabBitset[38].NbElements++;
		      SetBitTo1 (TabBitset[38].bitset, myligne);
		    }
		  else
		    {
		      TabBitset[39].NbElements++;
		      SetBitTo1 (TabBitset[39].bitset, myligne);
		    }
		  break;
		}
	    }
	  colonne++;
	}
      while (!done && !feof (f));
      myligne++;		/* on va lire une nouvelle ligne dans le fichier texte */
    }
  fclose (f);

  /* Affichage des résulats */
  max = -1;
  j = 0;
  for (i = 0; i < 40; i++)
    {
      if (TabBitset[i].NbElements > SUPPORT)
	{
	  printf ("%d: %d\n", i, TabBitset[i].NbElements);
	  max = j;
	  Candidate[j++] = i;	/* On range les k-itemset */
	}
    }
  printf ("Nb 1-itemset=%d, rc=%d\n", max, rc);

  /* On passe maintenant a la generation des 2-itemset */

  TOTAL = 0;			/* nombre total déléments implqués dans les intersections */
  n = max;
  k = 2;
  max1 = 0;
  prev = 0;
  lastComb = lexLastComb (n, k);
  j = 0;

  for (comb = lexFirstComb (k);; comb = lexNextComb (comb PARITY_PARAM))
    {
      unsigned i;
      bitmask_t p = lexPrevComb (comb);
      if (prev != p)
	printf ("Error: p = %llx, prev = %llx\n", p, prev);
      prev = comb;
      m = 0;
      for (i = 0; i < n; ++i)
	if (comb & Lshift1 (i))
	  {
	    /* on sélectionne les 2 bitsets */
	    //printf("%2d ", i);
	    MyBitset[m++] = i;
	  }
      /* On calcule le And sur les 2 bitsets */
      TOTAL += TabBitset[Candidate[MyBitset[0]]].NbElements +
	TabBitset[Candidate[MyBitset[1]]].NbElements;
      pandSIMDAMDOptCopy (TabBitset[Candidate[MyBitset[0]]].bitset,
		 TabBitset[Candidate[MyBitset[1]]].bitset,
		 TabBitsetCopy[j].bitset, 131072 * 8);

#ifdef COMPUTE_PARITY
      printf ("parity = %d", parity);
#endif
      //printf("\n");
      //printf("Nb Bit à 1: %d\n",NumberBit1(MyRes,131072*8));
      m = NumberBit1 (TabBitsetCopy[j].bitset, 131072 * 8);
      if (m > SUPPORT)
	{
	  j++;
	  printf ("%d-itemset: (%d, %d) of size: %d\n", k,
		  Candidate[MyBitset[0]], Candidate[MyBitset[1]], m);

	  TabBitsetCopy[max1].i = Candidate[MyBitset[0]];
	  TabBitsetCopy[max1].j = Candidate[MyBitset[1]];
	  TabBitsetCopy[max1].NbElements = m;

	  max1++;		/* il y a un itemset de plus */

	};

      if (comb == lastComb)
	break;
    }

  printf ("Nb 2-itemset=%d, rc=%d\n", max1, rc);

  /* Generation des 3-itemsets */
  /* Algo est le suivant :     */
  /* for tous les k de 1 à max1 faire */
  /*     - considerer le 2 itemset a la position k dans TabBitsetCopy */
  /*     en particulier les champs i et j qui donnent les deux 1-itemsets */
  /*     generateur. */
  /*     - verifier que tous les Candidate[p] (0< p< max) */
  /*     ne comportent pas les champs i et j ci dessus */
  /*     - Verifier que les 3 items choisis n'apparaissent pas deja (dans */
  /*     un autre ordre) : le triplet est un nouveau triplet ? */
  /*     Si c'est le cas, on procede a l'intersection et on garde si */
  /*     le support est suffisant */
  max2 = 0;
  for (k = 0; k < max1; k++)
    {
      i = TabBitsetCopy[k].i;
      j = TabBitsetCopy[k].j;
      for (l = 0; l < max; l++)
	{
	  if ((Candidate[l] != i) && (Candidate[l] != j))
	    {
              /* Verification que le triplet n'a pas deja ete genere */
              verif = 1; ind = 0;
              while(verif && (ind < max2)) { 
		if((TabBitsetOne[ind].i == i) || (TabBitsetOne[ind].i == j) || (TabBitsetOne[ind].i == Candidate[l]))
                  if((TabBitsetOne[ind].j == i) || (TabBitsetOne[ind].j == j) || (TabBitsetOne[ind].j == Candidate[l]))
		    if((TabBitsetOne[ind].l == i) || (TabBitsetOne[ind].l == j) || (TabBitsetOne[ind].l == Candidate[l])) verif = 0;
                ind++;
              }
              /* fin de la verification */
              if(verif) {
	      /* on peut proceder a l'intersection */
	      TOTAL += TabBitset[Candidate[l]].NbElements +
		TabBitsetCopy[k].NbElements;
	      pandSIMDAMDOptCopy (TabBitset[Candidate[l]].bitset,
			 TabBitsetCopy[k].bitset,
			 TabBitsetOne[max2].bitset, 131072 * 8);
	      

#ifdef COMPUTE_PARITY
	      printf ("parity = %d", parity);
#endif

	      m = NumberBit1 (TabBitsetOne[max2].bitset, 131072 * 8);
	      //              printf ("3-itemset: (%d, %d, %d) of size: %d\n",
	      //                  i, j, Candidate[l], m);
	      if (m > SUPPORT)
		{
		  printf ("3-itemset: (%d, %d, %d) of size: %d\n",
			  i, j, Candidate[l], m);
		  TabBitsetOne[max2].i = i;
		  TabBitsetOne[max2].j = j;
		  TabBitsetOne[max2].l = Candidate[l];
		  TabBitsetOne[max2].NbElements = m;
		  max2++;	/* il y a un itemset de plus */

		};
	      }
	    }
	}
    }

  printf ("Nb 3-itemsets=%d, rc=%d\n", max2, rc);

  /* Generation des 4-itemsets */
  /* Algo est le suivant :     */
  /* for tous les k de 1 à max1 faire */
  /*     - considerer le 3 itemset a la position k dans TabBitsetCopy */
  /*     en particulier les champs i, j, l qui donnent les 3 1-itemsets */
  /*     generateurs. */
  /*     - verifier que tous les Candidate[p] (0< p< max) */
  /*     ne comportent pas les champs i, j, l ci dessus */
  /*     - Verifier que les 4 items choisis n'apparaissent pas deja (dans */
  /*     un autre ordre) : le 4-itemset est un nouveau 4-itemset ? */
  /*     Si c'est le cas, on procede a l'intersection et on garde si */
  /*     le support est suffisant */
  max1 = 0;
  for (k = 0; k < max2; k++)
    {
      i = TabBitsetOne[k].i;
      j = TabBitsetOne[k].j;
      p = TabBitsetOne[k].l;
      for (l = 0; l < max; l++)
	{
	  //      printf("******* k=%d l=%d**********\n",k,l);
	  //printf("**** i=%d j=%d p=%d candidate=%d\n",i,j,p,Candidate[l]);
	  //if((i==0)&&(j==0)&&(p==0)) {printf("k=%d\n");exit(0);} 

	  if ((Candidate[l] != i) && (Candidate[l] != j)
	      && (Candidate[l] != p))
	    {

              /* Verification que le triplet n'a pas deja ete genere */
              verif = 1; ind = 0;
              while(verif && (ind < max1)) { 
		if((TabBitsetCopy[ind].i == i) || (TabBitsetCopy[ind].i == j) || (TabBitsetCopy[ind].i == Candidate[l]) ||  (TabBitsetCopy[ind].i == p))
                  if((TabBitsetCopy[ind].j == i) || (TabBitsetCopy[ind].j == j) || (TabBitsetCopy[ind].j == Candidate[l])  ||  (TabBitsetCopy[ind].j == p))
		    if((TabBitsetCopy[ind].l == i) || (TabBitsetCopy[ind].l == j) || (TabBitsetCopy[ind].l == Candidate[l])  ||  (TabBitsetCopy[ind].l == p)) verif = 0;
                ind++;
              }
              /* fin de la verification */
              if(verif) {
	      /* on peut proceder a l'intersection */
	      TOTAL += TabBitset[Candidate[l]].NbElements +
		TabBitsetOne[k].NbElements;
	      pandSIMDAMDOptCopy (TabBitset[Candidate[l]].bitset,
			 TabBitsetOne[k].bitset,
			 TabBitsetCopy[max1].bitset, 131072 * 8);

	      m = NumberBit1 (TabBitsetCopy[max1].bitset, 131072 * 8);
	      //printf ("4-itemset: (%d, %d, %d, %d) of size: %d\n",
	      //i, j, p, Candidate[l], m);
	      if (m > SUPPORT)
		{
		  printf ("4-itemset: (%d, %d, %d, %d) of size: %d\n",
			  i, j, p, Candidate[l], m);
		  TabBitsetCopy[max1].i = i;
		  TabBitsetCopy[max1].j = j;
		  TabBitsetCopy[max1].l = p;
		  TabBitsetCopy[max1].m = Candidate[l];
		  TabBitsetCopy[max1].NbElements = m;
		  max1++;	/* il y a un itemset de plus */

		};
	      }
	    }
	}
    }

  printf ("Nb 4-itemsets=%d, rc=%d\n", max1, rc);
  printf
    ("Total number of lines (elements) involved in intersections : %ld\n",
     TOTAL);
  exit (0);
}