/*

Progam: Autocorr.C computes the set of valid autocorrelation
vectors for strings of length <q>,  
C++ version with BitString implementation

Author: Eric Rivals
Date:   12/08/1999

% Author: Eric Rivals
% Address: LIRMM, 161, rue Ada, F-34392 Montpellier Cedex 5
% Email: rivals@lirmm.fr


Parameters: 
- <q>: length of strings of which the autocorrelations are wanted
- <print_kappas>: flag to print out the kappas, the number of autocorrelations
                  for a given size
 */


#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>
#include <BitString.h>

//#define _DEBUG 
//#define _DEBUG_TERM 
//#define _DEBUG_BS 
#define    INIT_NB_VECTORS -1


typedef BitString *pBS;

int    *NbVectors;
int    *K2;
int    NbVQ = 0;     /* nb of vectors of size q */
int    K2Q = 0;      /* K2 for size q */
pBS    *Vectors;

int    FirstNotYetComputedQ = 4;
const char Kappas_fname1[] = "Kappas.txt";
const char Kappas_fname2[] = "Kappas_";

ofstream OsKappasGen;

//////////////////////////////////////////////////////////////////////
//                    InitAutocorrelations
//////////////////////////////////////////////////////////////////////
void InitAutocorrelations (int q){

    int j;

#ifdef _DEBUG_INIT
    fprintf(stderr,"InitAutocorrelations Begin\n");
#endif

    NbVectors = new int[q+1];
    K2 = new int[q+1];
    Vectors = new pBS[q+1];

    /* initialise K2 and NbVectors */
    K2[0] = 0; K2[1] = 1; K2[2] = 1; K2[3] = 2;

    for (j=FirstNotYetComputedQ; j <= q; ){
      K2[j] = 0;
      NbVectors[j++] = INIT_NB_VECTORS;
    }

    /* initialise for word sizes 0, 1, 2, 3 */
    NbVectors[0] = 1;

    NbVectors[1] = 1;
    Vectors[1] = new BitString[1];
    Vectors[1][0].set(0);            /* 1 */

#ifdef _DEBUG_INIT
    fprintf(stderr,"size 1 fini\n");
#endif

    NbVectors[2] = 2;
    Vectors[2] = new BitString[2];
    Vectors[2][0].set(0,1);          /* 11 */
    Vectors[2][1].clear(1);          /* 10 */
    Vectors[2][1][0] = 1;            // ATTENTION such [] access cannot used as
                                     // the first access to a bit
#ifdef _DEBUG_INIT
    fprintf(stderr,"size 2 fini\n");
#endif

    NbVectors[3] = 3;
    Vectors[3] = new BitString[3];
    Vectors[3][0].set(0,2);          /* 111 */
    Vectors[3][1].set(0,2);          /* 101 */
    Vectors[3][1][1] = 0;
    Vectors[3][2].set(0);            /* 100 */
    Vectors[3][2].clear(1,2);    

#ifdef _DEBUG_INIT
    fprintf(stderr,"InitAutocorrelations End\n");
#endif

    return;
}

//////////////////////////////////////////////////////////////////////
// Normal            ComputeAutocorrelations( int q, int PrintKappas )
//////////////////////////////////////////////////////////////////////
int ComputeAutocorrelations( int q, int PrintKappas ){

    int 
      i,            /* index for VectorsQ */
      ivnq,         /* index for Vectors vector_index_new_q } */
      p = 1,        /* basic period : first non-zero period  */
      r = 0,        // rest of the division of q by p
      j, k, l,
      half_q = q/2, /* integer half of q */
      beta = 0,     /* [(greatest multiple of p) < q] minus p */
      new_q = 0;    // size of nested vectors

    int K1 = 1;          // to compute K1
    ofstream os;

    BitString RefAC; /* AC of reference, prefix common to each ac with same basic period*/
    int
      RefAC_intermediate_size = 0, //real size for RefAC prefix vector in inner loop of case 1
      RefAC_size = 0;  // real size used for the current recursion
    

#ifdef _DEBUG
    fprintf(stderr,"ComputeAutocorrelations %d\n",q);
#endif	

    // Compute the set of Autocorrelations \Gamma_i for i in [2, half_q]
    // and compute the upper bound of the number of needed vectors
    for ( new_q = FirstNotYetComputedQ; new_q <= half_q; new_q++ )
      if (NbVectors[new_q] == INIT_NB_VECTORS) {// Gamma[new_q] to be computed
	FirstNotYetComputedQ = ComputeAutocorrelations( new_q, PrintKappas );
      }

    if (PrintKappas){
      char fname[50];sprintf(fname,"%s%d",Kappas_fname2,q);
      os.open(fname,ios::out); 
      if (!os){
	cerr << "Error: cannot open out file " << fname << endl;
	exit(1);
      }
    }
    
    // compute K2[q]
    if (FirstNotYetComputedQ % 2 == 0){
      K2[FirstNotYetComputedQ] = K2[FirstNotYetComputedQ-1];
      j = FirstNotYetComputedQ+1;
    }
    else 
      j = FirstNotYetComputedQ;
    k = j/2;
    while ( j <= q ){
      K2[j] = NbVectors[k] + K2[j-1];
      j++; k++;
      K2[j] = K2[j-1];
      j++;
    }

#ifdef _DEBUG
    fprintf(stderr,"K2[%d] = %d\n", q, K2[q]);
    fprintf(stderr,"END Compute sub q\n");
#endif	

    // allocate as much vectors for size q
    Vectors[q] = new BitString[2*K2[q]];

    //////////////////////////////////////////////////
    //                          Case a :  p <= q/2
    //////////////////////////////////////////////////

    /* 11...1 : initialise for p = 1; */
    Vectors[q][0].set(0,q-1);
    i = NbVectors[q] = 1;


    for ( p=2; p <= half_q; p++ ) { 

	r = q%p;
        beta = q - p - r;	/* INV beta is a multiple of p */


#ifdef _DEBUG_BS
	fprintf(stderr,"1 RefAC of size beta %d\n", beta);
#endif
	// create reference vector RefAC;
	RefAC_size = beta;
	RefAC.clear(0,RefAC_size-1);
	RefAC[0] = 1;
	for ( k = p; k < beta; k += p ) RefAC[k] = 1;  
#ifdef _DEBUG_BS
	fprintf(stderr,"end RefAC\n");
	fprintf(stderr,"p %d beta %d RefAC ", p,beta);
	for(k=0;k < RefAC.length(); k++)fprintf(stderr,"%d ", RefAC.test(k));
	fprintf(stderr,"\n");
#endif

	if (r == 0){                  // new_q == p

	  for ( ivnq = 1; ivnq < NbVectors[p]-1; ivnq++ ){ 
	    // 1 because 1...1 is not a valid nested vector
	    // NbVectors[p]-1  because last nested vector is 10...0 and is always used

#ifdef _DEBUG_BS
	    fprintf(stderr,"Vectors[%d][%d]: ", p,  ivnq);
	    for(k=0;k < p; k++)fprintf(stderr,"%d ", Vectors[p][ivnq].test(k));
	    fprintf(stderr,"\n");
#endif
	    /*search for the first bit to one : from 2 because vector 11..1 was skipped 
	      no check if j < p because there is at least  a bit to one */
	    j = Vectors[p][ivnq].index(1,2); // search the first non-zero position
	    // 2 because the first nested vector 1...1 is not considered in the loop
 
	    // INVARIANT  ( j < p ) 
	    if ( (p % j) == 0) continue;
	    // INV current vector of length p is valid 

	    /* create and prepare Vectors[q][i] */
	    cat(RefAC.before(RefAC_size), Vectors[p][ivnq], Vectors[q][i]);
#ifdef _DEBUG_BS
	    fprintf(stderr,"p %d beta %d Vectors ", p,beta);
	    for(k=0;k < Vectors[q][i].length(); k++)
	      fprintf(stderr,"%d ", Vectors[q][i].test(k));
	    fprintf(stderr,"\n");
#endif

	    i++;     /* increment the index of Vectors[q] */
	  }/*   endfor ivnq */

	  /* create vector for the nested vector 100...0 */
	  cat(RefAC.before(RefAC_size), Vectors[p][ivnq], Vectors[q][i]);
	  i++;     /* increment the index of Vectors[q] */

	} // endif r == 0
	else{                  // r > 0 

	  // For all indices j<p which can be a basic period of suffix of lg (p+r)
	  for ( j = r+2; j < p; j++){
	 
	    if (p%j == 0) continue;    // skip j if it divides p

	    // update RefAC according to j
	    RefAC_intermediate_size = RefAC_size + j;
	    RefAC.set(beta);
	    RefAC.clear(beta+1,RefAC_intermediate_size-1);
	    new_q = p + r - j;

#ifdef _DEBUG
	    fprintf(stderr,"RefAC is %d\n", RefAC_intermediate_size);
	    fprintf(stderr,"p %d r %d j %d\n", p, r, j);
#endif
	    // visit all vectors of Gamma[new_q]; 
	    //  except 10...0 => NbVectors[new_q]- 1 because  p-j>0 and bit (p-j) should be 1
	    for ( ivnq = 0; ivnq < NbVectors[new_q]-1; ivnq++ ){ 

	      if ( Vectors[new_q][ivnq][p-j] == 0 ) continue; /* vector should have period p*/

	      if ( j <= ((p+r)/2) ){

		// search for the first bit to one : from index 1
		// no check if l < new_q because there is at least  a bit to one
		l = Vectors[new_q][ivnq].index(1,1); 
#ifdef _DEBUG
		fprintf(stderr,"l %d\n", l);
#endif
		// if (l == j) OK
		if (l > j) break;
		/* INVARIANT  ( l < new_q <p ) */
		if ((l < j) && (((j % l) == 0) || 
				((j < (p+r-j)) && (Vectors[new_q][ivnq][j] == 0)))) continue;
	      }

#ifdef _DEBUG_BS
	      fprintf(stderr,"new_q %d ivnq %d p %d Vectors ", new_q, ivnq, p);
	      for(k=0;k < Vectors[new_q][ivnq].length(); k++)
		fprintf(stderr,"%d ", Vectors[new_q][ivnq].test(k));
	      fprintf(stderr,"\n");
#endif

	      /* INVARIANT current nested vector of length new_q is valid */ 
	      /* create Vectors[q][i] */
	      cat(RefAC.before(RefAC_intermediate_size), Vectors[new_q][ivnq], Vectors[q][i]);

#ifdef _DEBUG_BS
	      fprintf(stderr,"new_q %d p %d beta %d lg %d Vectors ", 
		      new_q, p,beta, Vectors[q][i].length());
	      for(k=0;k < Vectors[q][i].length(); k++)
		fprintf(stderr,"%d ", Vectors[q][i].test(k));
	      fprintf(stderr,"\n");
#endif

	      i++;     /* increment the index of Vectors[q] */
	    
	    }/*   endfor ivnq */
	  }// endfor j

	  // case j == p, and new_q == r
	  // update RefAC
	  RefAC_intermediate_size = RefAC_size + p;
	  RefAC.set(beta);
	  RefAC.clear(beta+1, RefAC_intermediate_size-1);

	  // visit all vectors of Gamma[r]
	  for ( ivnq = 0; ivnq < NbVectors[r]; ivnq++ ){ 

#ifdef _DEBUG_BS
	    fprintf(stderr,"r %d ivnq %d p %d Vectors ", r, ivnq, p);
	    for(k=0;k < Vectors[r][ivnq].length(); k++)
	      fprintf(stderr,"%d ", Vectors[r][ivnq].test(k));
	    fprintf(stderr,"\n");
#endif

	    /* INVARIANT current nested vector of length r is valid */		
	    /* create Vectors[q][i] */
	    cat(RefAC.before(RefAC_intermediate_size), Vectors[r][ivnq], Vectors[q][i]);
#ifdef _DEBUG_BS
	    fprintf(stderr,"r %d p %d beta %d lg %d Vectors ", 
		    r, p,beta, Vectors[q][i].length());
	    for(k=0;k < Vectors[q][i].length(); k++)
	      fprintf(stderr,"%d ", Vectors[q][i].test(k));
	    fprintf(stderr,"\nEND\n");
#endif

	    i++;     /* increment the index of Vectors[q] */
	    
	  }/*   endfor ivnq */
	} // end else r != 0

#ifdef _DEBUG
	fprintf(stderr,"p %d i %d\n", p, i);
#endif
	if (PrintKappas){
	  os << q << " " << p << " " << (q%p) << " " << i-NbVectors[q] << " " << NbVectors[p] << endl;
	  K1 += i-NbVectors[q];
	  // 	fprintf(stdout,"NbVectors[ %5d %5d ] %6d\n",q,p,i-NbVectors[q]);
	}
	NbVectors[q] = i;  /* update the nb of vectors of length q */

    } /* endfor p */


    //////////////////////////////////////////////////
    //                          Case b :  p > q/2
    //////////////////////////////////////////////////

    // INVARIANT: 	NbVectors[q] == i && p == half_q+1
    for ( ; p < q; p++ ) { // { p > half_q }

	new_q = q-p;

	if (PrintKappas){
	  os << q << " " << p << " " << NbVectors[new_q] << endl;
	  // 	fprintf(stdout,"NbVectors[ %5d %5d ] %6d\n",q,p,NbVectors[new_q]);
	}
#ifdef _DEBUG_BS
	fprintf(stderr,"RefAC of size p %d\n", p);
#endif
	// create reference vector RefAC
	RefAC_size = p;
	RefAC.clear(1,RefAC_size-1);
	RefAC[0] = 1;

#ifdef _DEBUG_BS
	fprintf(stderr,"end RefAC\n");
	fprintf(stderr,"new_q %d p %d  RefAC ", new_q, p);
	for(k=0;k < RefAC.length(); k++)fprintf(stderr,"%d ", RefAC.test(k));
	fprintf(stderr,"\n");
#endif
	    
	for ( ivnq = 0; ivnq < NbVectors[new_q]; ivnq++, i++ ){ // build a vector for every nested vector
	  // create and prepare Vectors[q][i]
	  cat(RefAC.before(RefAC_size), Vectors[new_q][ivnq], Vectors[q][i]);

#ifdef _DEBUG_BS
	  fprintf(stderr,"new_q %d p %d Vectors ", new_q, p);
	  for(k=0;k < Vectors[q][i].length(); k++)
	    fprintf(stderr,"%d ", Vectors[q][i].test(k));
	  fprintf(stderr,"\n");
#endif

	}// endfor ivnq
    
	NbVectors[q] +=  NbVectors[new_q];  //  update Nb of vectors of length q 
    } // endfor p 

    // 10...0 : case p == q
    Vectors[q][i].set(0);    Vectors[q][i].clear(1,q-1);
#ifdef _DEBUG_BS
       fprintf(stderr,"Vector p == q ");
       for(k=0;k < Vectors[q][i].length(); k++)fprintf(stderr,"%d ", Vectors[q][i].test(k));
       fprintf(stderr,"\n");
#endif

    i++;
    NbVectors[q]++;

#ifdef _DEBUG
    fprintf(stderr,"End ComputeAutocorrelations %d %d\n", q, NbVectors[q]);
#endif
    if (PrintKappas){	
      os << q << " " << NbVectors[q] << " " << K1 << " " << K2[q] << endl;
      //     fprintf(stdout,"\nNbVectors[ %5d ]      %6d\n",q,NbVectors[q]);
      os.close();
    }

    return q;
   
}

//////////////////////////////////////////////////////////////////////
//        ComputeTerminalAC(int q, int PrintKappas, int OptionSaveQ )
//////////////////////////////////////////////////////////////////////

int ComputeTerminalAC( int q, int PrintKappas, int OptionSaveQ ){

    int 
      i,            // index for VectorsQ 
      new_q = 0,    // size of nested vectors: new_q <= q/2 + 1
      ivnq,         // index for Vectors of size new_q vector_index_new_q 
      p = 1,        // basic period of current vector (first non-zero period)
      r = 0,        // rest of the division of q by p
      j, k, l,
      half_q = q/2, // integer half of q
      //lowest_stored_sub_q = 2*(half_q/3)+half_q%3 - 1,
      beta = 0;     /* [(greatest multiple of p) < q] minus p */


    int K1 = 1;          // to compute K1
    ofstream os;
    BitString RefAC;     // AC of reference, prefix common to each ac with same basic period
    int
      RefAC_intermediate_size = 0, //real size for RefAC prefix vector in inner loop of case 1
      RefAC_size = 0;  // real size used for the current recursion
    char 
      fnameK[50],
      fname[50];
    FILE
      *fout;

    if ( OptionSaveQ ) { //  save only q to file
      strcpy(fname,"AutocorrelationVectors.");
      sprintf(fname+23,"%d",q);
      fout = fopen(fname,"w");       // open output file 

      /* print periods */
      for ( k=0; k < q; ) fprintf(fout,"%d ", (k++)%10 );
      fprintf(fout,"\n\n");

    }


    // Compute the set of Autocorrelations \Gamma_i for i in [2, half_q]
    // and compute the upper bound of the number of needed vectors
    for ( new_q = FirstNotYetComputedQ; new_q <= half_q; new_q++ )
      if (NbVectors[new_q] == INIT_NB_VECTORS) {/* { Vectors for length new_q not already computed } */
	FirstNotYetComputedQ = ComputeAutocorrelations( new_q, PrintKappas );
      }

#ifdef _DEBUG_TERM
    fprintf(stderr,"End ComputeAutocorrelations for nested AC until %d\n",half_q);
#endif	

    if (PrintKappas){
      sprintf(fnameK,"%s%d",Kappas_fname2,q);
      os.open(fnameK,ios::out); 
      if (!os){
	cerr << "Error: cannot open out file " << fnameK << endl;
	exit(1);
      }
    }

#ifdef _DEBUG_TERM
    fprintf(stderr,"ComputeAutocorrelations %d\n",q);
#endif	

    //////////////////////////////////////////////////
    //                          Case a :  p <= q/2
    //////////////////////////////////////////////////


    // 1...1 : case p == 1
    RefAC.set(0,q-1);           // set all bits to one in a single operation

    if (OptionSaveQ){
      // print current vector 
      for ( k=0; k < q; ) fprintf(fout,"%d ", RefAC.test(k++));
      fprintf(fout,"\n");
    }

    i = NbVQ = 1;

    for ( p=2; p <= half_q; p++ ) { 

      r = q % p;                      
      beta = q - p - r;	              // INV beta is a multiple of p

#ifdef _DEBUG_TERM
    fprintf(stderr,"Loop p == %d\n",p);
#endif	
 

#ifdef _DEBUG_BS
      fprintf(stderr,"RefAC of size beta %d\n", beta);
#endif
      // create reference vector RefAC; set bit beta to one not needed because already in the nested vector
      RefAC_size = beta;
      RefAC.clear(1,RefAC_size-1);
      // CHECK
      for ( k = p; k < beta; k += p ) RefAC[k] = 1;  
#ifdef _DEBUG_BS
      fprintf(stderr,"end RefAC\n");
      fprintf(stderr,"new_q %d p %d beta %d RefAC ", new_q, p,beta);
      for(k=0;k < RefAC.length(); k++)fprintf(stderr,"%d ", RefAC.test(k));
      fprintf(stderr,"\n");
#endif

      /*  Separate in two subcases (new_q == p) and (new_q > p) */ 
      if ( r == 0 ) {                  //new_q == p

	for ( ivnq = 1; ivnq < NbVectors[p]-1; ivnq++ ){
	  // 1 because 1...1 already treated:  it is not a valid nested vector
	  // NbVectors[p]-1  because last nested vector is 10...0 and is always used


#ifdef _DEBUG_BS
	  fprintf(stderr,"Vectors[%d newq][%d ivnq]: ", p, p);
	  for(k=0;k < p; k++)fprintf(stderr,"%d ", Vectors[p][ivnq].test(k));
	  fprintf(stderr,"\n");
#endif
	  /*search for the first bit to one : from 2 because vector 11..1 was skipped 
	    no check if j < p because there is at least  a bit to one */
	  j = Vectors[p][ivnq].index(1,2); // search the first non-zero position
	  /* INVARIANT  ( j < p ) */
	  if ( (p % j) == 0) continue;
	  /* INV current vector of length p is valid */
		
	  /* create and prepare Vectors[q][i] */
	  cat(RefAC.before(RefAC_size), Vectors[p][ivnq], RefAC);

	  if (OptionSaveQ){	    // print current vector 
	    for ( k=0; k < q; ) fprintf(fout,"%d ", RefAC.test(k++));
	    fprintf(fout,"\n");
	  }

	  i++;     /* increment the index of Vectors[q] */
	}/*   endfor ivnq */

	/* nested vector  10...0 */
	cat(RefAC.before(RefAC_size), Vectors[p][ivnq], RefAC);

	if (OptionSaveQ){	  // print current vector 
	  for ( k=0; k < q; ) fprintf(fout,"%d ", RefAC.test(k++));
	  fprintf(fout,"\n");
	}


	i++;     /* increment the index of Vectors[q] */

      }// end if r == 0

      else{                  // r > 0 

	// For all indices j<p which can be a basic period of suffix of lg (p+r)
	for ( j = r+2; j < p; j++){
	 
	  if (p%j == 0) continue;    // skip j if it divides p

	  // update RefAC according to j
	  RefAC_intermediate_size = RefAC_size + j;
	  RefAC.set(beta);
	  RefAC.clear(beta+1,RefAC_intermediate_size-1);
	  new_q = p + r - j;

#ifdef _DEBUG
	  fprintf(stderr,"RefAC is %d\n", RefAC_intermediate_size);
	  fprintf(stderr,"p %d r %d j %d\n", p, r, j);
#endif
	  // visit all vectors of Gamma[new_q]; 
	  //  except 10...0 => NbVectors[new_q]- 1 because  p-j>0 and bit (p-j) should be 1
	  for ( ivnq = 0; ivnq < NbVectors[new_q]-1; ivnq++ ){ 

	    if ( Vectors[new_q][ivnq][p-j] == 0 ) continue; /* vector should have period p*/

	    if ( j <= ((p+r)/2) ){

	      // search for the first bit to one : from index 1
	      // no check if l < new_q because there is at least  a bit to one
	      l = Vectors[new_q][ivnq].index(1,1); 
#ifdef _DEBUG
	      fprintf(stderr,"l %d\n", l);
#endif
	      // if (l == j) OK
	      if (l > j) break;
	      /* INVARIANT  ( l < new_q <p ) */
	      if ((l < j) && (((j % l) == 0) || 
			      ((j < (p+r-j)) && (Vectors[new_q][ivnq][j] == 0)))) continue;
	    }

#ifdef _DEBUG_BS
	    fprintf(stderr,"new_q %d ivnq %d p %d Vectors ", new_q, ivnq, p);
	    for(k=0;k < Vectors[new_q][ivnq].length(); k++)
	      fprintf(stderr,"%d ", Vectors[new_q][ivnq].test(k));
	    fprintf(stderr,"\n");
#endif

	    /* INVARIANT current nested vector of length new_q is valid */ 
	    // create current vectors of size q
	    cat(RefAC.before(RefAC_intermediate_size), Vectors[new_q][ivnq], RefAC);

	    i++;     /* increment the index of Vectors[q] */
	    
	  }/*   endfor ivnq */
	}// endfor j

	// case j == p, and new_q == r
	// update RefAC
	RefAC_intermediate_size = RefAC_size + p;
	RefAC.set(beta);
	RefAC.clear(beta+1, RefAC_intermediate_size-1);

	// visit all vectors of Gamma[r]
	for ( ivnq = 0; ivnq < NbVectors[r]; ivnq++ ){ 

#ifdef _DEBUG_BS
	  fprintf(stderr,"r %d ivnq %d p %d Vectors ", r, ivnq, p);
	  for(k=0;k < Vectors[r][ivnq].length(); k++)
	    fprintf(stderr,"%d ", Vectors[r][ivnq].test(k));
	  fprintf(stderr,"\n");
#endif

	  /* INVARIANT current nested vector of length r is valid */		
	  /* create Vectors[q][i] */
	  cat(RefAC.before(RefAC_intermediate_size), Vectors[r][ivnq], RefAC);

	  i++;     /* increment the index of Vectors[q] */
	    
	}/*   endfor ivnq */
      } // end else r != 0

// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//        else { /* new_q > p */

//  	for ( ivnq = 1; ivnq < NbVectors[new_q]-1; ivnq++ ){
//  	  // 1 because 1...1 already treated:  it is not a valid nested vector
//  	  // NbVectors[new_q]-1  because last nested vector is 10...0 and not valid here

//  	  if ( Vectors[new_q][ivnq][p] == 0 ) continue; // vector should have period p

//  	  /*search for the first bit to one : from 2 because vector 11..1 was skipped
//  	    no check if j < new_q because there is at least  a bit to one */
//  	  j = Vectors[new_q][ivnq].index(1,2); // search the first non-zero position

//  	  if (j > p) break; // because no vectors with j > p can contain p which is required
//  	  else 
//  	    if ((j<p) && ((p % j) == 0)) continue;
	  
//  	  /* INVARIANT current nested vector of length new_q is valid */		
//  	  /* create Vectors[q][i] */
//  	  cat(RefAC.before(RefAC_size), Vectors[new_q][ivnq], RefAC);

//  	  if (OptionSaveQ){
//  	    // print current vector 
//  	    for ( k=0; k < q; ) fprintf(fout,"%d ", RefAC.test(k++));
//  	    fprintf(fout,"\n");
//  	  }

//  	  i++;     /* increment the index of Vectors[q] */

//  	}/*   endfor ivnq */
//        }/*    end for else */
//======================================================================
      if (PrintKappas){
	os << q << " " << p << " " << (q%p) << " " << i-NbVQ << " " << NbVectors[p] << endl;
	K1 += i-NbVQ;
      }
      NbVQ = i;  /* update the nb of vectors of length q */

    } /* endfor p */


#ifdef _DEBUG_TERM
    fprintf(stderr,"End case 1  p <= q %d\n",q);
#endif	


    //////////////////////////////////////////////////
    //                          Case b :  p > q/2
    //////////////////////////////////////////////////

    /* INVARIANT: 	NbVQ == i && p == half_q+1 */
    for ( ; p < q; p++ ) { /* { p > half_q } */

      new_q = q-p;
      
      if (PrintKappas){
	os << q << " " << p << " " << NbVectors[new_q] << endl;
      }

#ifdef _DEBUG_BS
      fprintf(stderr,"RefAC of size p %d\n", p);
#endif
      // create reference vector RefAC
      RefAC_size = p;
      RefAC.clear(RefAC_size-1);
      RefAC.clear(); 
      RefAC[0] = 1;
#ifdef _DEBUG_BS
      fprintf(stderr,"end RefAC\n");
      fprintf(stderr,"new_q %d p %d  RefAC ", new_q, p);
      for(k=0;k < RefAC.length(); k++)fprintf(stderr,"%d ", RefAC.test(k));
      fprintf(stderr,"\n");
#endif
	    
      // build a vector for each nested vector
      for ( ivnq = 0; ivnq < NbVectors[new_q]; ivnq++, i++){
	 
	// create and prepare Vectors[q][i]
	cat(RefAC.before(RefAC_size), Vectors[new_q][ivnq], RefAC);

	if (OptionSaveQ){
	  // print current vector 
	  for ( k=0; k < q; ) fprintf(fout,"%d ", RefAC.test(k++));
	  fprintf(fout,"\n");
	}

      }/* endfor ivnq */

      K2Q +=   NbVectors[new_q];
      NbVQ +=  NbVectors[new_q];  /*  update Nb of vectors of length q */

//        if (new_q > lowest_stored_sub_q){ // deallocate Vectors
//  	delete [] (Vectors[new_q]);
//        }

    } /* endfor p */

    // 10...0 : case p == q
    RefAC.clear(1,q-1); RefAC[0] = 1;
    NbVQ++;

    if (OptionSaveQ){
      // print current vector 
      for ( k=0; k < q; ) fprintf(fout,"%d ", RefAC.test(k++));
      fprintf(fout,"\n");

      fclose(fout);
    }

#ifdef _DEBUG
    fprintf(stderr,"End ComputeAutocorrelations %d %d\n", q, NbVQ);
#endif
    if (PrintKappas){
      os << q << " " << NbVQ << " " << K1 << " " << K2Q << endl;
      os.close();
    }

    return q;
   
}

int main(int argc, char* argv[]){

    int 
      i,j,k,
      OptionAll = 0,
      OptionSaveQ = 0,
      OptionSaveAll = 0,
      PrintKappas = 0,
      highest_sub_q = 0,
      q = 0;
    char
      fname[50];
    const char
      StringOptionSaveQ[8] = "-save_q",
      StringOptionSaveAll[10] = "-save_all",
      StringOptionAll[5] = "-all";
    FILE
      *fout;

    if (( argc < 3 ) || (atoi(argv[1]) < 1))  {
	fprintf(stderr,"usage: %s word_size(>1) print_kappas(0/1) [-all] [-save_q | -save_all]\n", argv[0]);
	fprintf(stderr,"       computes the set of autocorrelation vectors for word size of length <word_size>\n");
	fprintf(stderr,"       outputs one file per word size with the all autocorrelation vectors and their number\n");
	fprintf(stderr,"       with option [-all] computes the set for all sizes until <word_size>\n");
	fprintf(stderr,"       with option [-save_q] save only the set for <word_size> q to file\n");
	fprintf(stderr,"       with option [-save_all] save all sets to files\n");

	exit(1);
    }
    q = atoi(argv[1]);
    if (atoi(argv[2]) == 0)          // argument to print kappas in files
      PrintKappas = 0;
    else 
      PrintKappas = 1;

    if (argc >= 4) { 
	if (!strcmp(argv[3],StringOptionAll)) {
	    OptionAll = 1;

	    if (argc >= 5) {
		// printf("4 %s\n", argv[4]);
		if (!strcmp(argv[4],StringOptionSaveAll))
		    OptionSaveAll = 1;
		else
		    if (!strcmp(argv[4],StringOptionSaveQ))
			OptionSaveQ = 1;
	    }
	}
	else{
	    if (!strcmp(argv[3],StringOptionSaveAll))
		OptionSaveAll = 1;
	    else
		if (!strcmp(argv[3],StringOptionSaveQ))
		    OptionSaveQ = 1;
	}
    }


    printf("PrintKappas %d OptionAll %d OptionSaveQ %d OptionSaveAll %d\n", 
	   PrintKappas, OptionAll, OptionSaveQ, OptionSaveAll);

    if (OptionAll) {
	printf("Compute ALL sets of autocorrelation vector for word size q GOING FROM 1 TO %d\n",q);
	highest_sub_q = q;
	InitAutocorrelations(highest_sub_q);
#ifdef _DEBUG
	fprintf(stderr," End InitAutocorrelations\n");
#endif

	for ( i = FirstNotYetComputedQ; i <= q; i++ ){ /* for all sizes */
	  FirstNotYetComputedQ = ComputeAutocorrelations(i, 0);
	  //FirstNotYetComputedQ = ComputeAutocorrelations(i, PrintKappas);
	  // options kappas supprimes pour l'instant
	}
    }
    else{
	printf("Compute the set of autocorrelation vector for word size q == %d\n",q);

	highest_sub_q = q/2;
	// highest_sub_q = 2*(q/3)+ q%3 -1; // previous algo's value
	InitAutocorrelations(highest_sub_q);

	FirstNotYetComputedQ = ComputeTerminalAC(q, 0, OptionSaveQ);
	//	FirstNotYetComputedQ = ComputeTerminalAC(q, PrintKappas, OptionSaveQ);
	// options kappas supprimes pour l'instant
    }

    if (PrintKappas){  // open output file for output of kappas
      OsKappasGen.open(Kappas_fname1,ios::out); 
      if (!OsKappasGen){
	cerr << "Error: cannot open out file " << Kappas_fname1 << endl;
	exit(1);
      }

      /* output general kappas  to file */
      OsKappasGen << "1 0 1 1" << endl;              // size 1
      OsKappasGen << "2 1 1 2" << endl;              // size 2
      OsKappasGen << "3 1 2 3" << endl;              // size 3

      for ( i=4 ; i <= highest_sub_q; i++ ){  /* for all sizes */
	if ( NbVectors[i] != INIT_NB_VECTORS ) {
	  OsKappasGen << i << " " << NbVectors[i] - K2[i] << " " << K2[i] << " " << NbVectors[i] << endl; // size i
	}
      }
      if (!OptionAll)
	OsKappasGen << q << " " << NbVQ - K2Q << " " << K2Q << " " << NbVQ << endl; // size q

      OsKappasGen.close();
    }



    if (!OptionSaveAll) /* no save to file */
	return 0;

    if ( !OptionAll ) q = q/2;

    strcpy(fname,"AutocorrelationVectors.");
    /* output results to file */
    for ( i=3 ; i <= q; i++ ){  /* for all sizes */

	if ( NbVectors[i] != INIT_NB_VECTORS ) {

	    /* open output file */
	    sprintf(fname+23,"%d",i); /*itoa(i, fname+23,10); */
	    fprintf(stdout,"fname %s %d %d\n", fname, q, NbVectors[i]);
	    fout = fopen(fname,"w");

	    fprintf(fout,"word size = %d  nb of autocorrelation vectors %d\n\n", i, NbVectors[i]);

	    /* print periods */
	    for ( k=0; k < i; ) fprintf(fout,"%d ", (k++)%10 );
	    fprintf(fout,"\n\n");

	    for ( j = 0; j < NbVectors[i]; j++ ){ /* for all vectors of that size */

		/* print jth vector */
		for ( k=0; k < i; ) fprintf(fout,"%d ", Vectors[i][j].test(k++));
		fprintf(fout,"\n");
	    }
	    fclose(fout);
	}
    }

    return 0;
}