/*

Progam: Autocorrelations.c computes the set of valid autocorrelation
vectors for strings of length <q>

Author: Eric Rivals
Date:   14/07/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

**********************************************************************
*                    DESCRIPTION OF THE ALGORITHM
**********************************************************************
 Algo: ComputeAutocorrelations

Variables:

paralell arrays  from 0 to q/2+1 = NB_SIZES

int NbVectors[NB_SIZES+1]
(Bitarr**) Vectors[NB_SIZES+1]

// same for size q
NbVectorsQ
(Bitarr**) VectorsQ


int
i     // index for VectorsQ
ivnq  // index for Vectors vector_index_new_q }
*Ones // dynamic array 
iones // index for Ones
p1    // first non-zero period 
j
half_q
beta  // [(greatest multiple of p1) < q] minus p1
new_q // new length of strings for which Autocorrelations are recursively computed

initialise for length 1

initialise for p1 = 1
half_q = q/2;

// new_q <= q/2 + 1;
create array Ones with q/2 + 1 places and initialise to 0

for p1 ( 2..half_q ) { p1 <= half_q }

  if ( q mod p1 == 0 )
    beta = q - p1;
  else
    beta = q - p1 - (q mod p1);
  endif
  new_q = q - beta;
  if (NbVectors[new_q] == NB_INIT_VALUE) { Vectors for length new_q not already computed }
    ComputeAutocorrelations( new_q );
  endif

  i = NbVectorsQ; 
  // create and prepare VectorsQ[i]
  create and clear VectorsQ[i];
  set bits 0, p1, 2p1, ..., beta of VectorsQ[i] to 1;
  if ( new_q > p1 )
    set bit beta+p1 of VectorsQ[i] to 1;
  endif
  
  for ( ivnq = 0; ivnq < NbVectors[new_q]; ivnq++ )

    iones = 0;
    if ( new_q > p1 ) && ( bit p1 of of Vectors[new_q][ivnq] to 0 )
      continue;
    endif

    // search for bit to 1 ( and store them) and if no period divides p1
    for ( j=1; j < new_q; j++ ) 
      if ( bit j of Vectors[new_q][ivnq] to 1 )
        if ( p1 mod j == 0 )
	  iones = 0;
	  continue;
	endif
        Ones[iones++] = j;
      endif
    endfor j

    // set to 1 the bits of VectorsQ  corresponding to beta+bits of Vectors[new_q][ivnq]
    for ( j=0; j < iones; j++ )
      set bit beta+Ones[j] of VectorsQ[i] to 1;
    endfor j

    i++; // increment the nb of VectorsQ
    // create and prepare VectorsQ[i]
    create and clear VectorsQ[i];
    set bits 0, p1, 2p1, ..., beta of VectorsQ[i] to 1;
    if ( new_q > p1 )
      set bit beta+p1 of VectorsQ[i] to 1;
    endif

  endfor ivnq
  NbVectorsQ = i;
  destroy or reset VectorsQ[i];

endfor p1


for p1 ( half_q..q-1 ) { p1 > half_q }

  new_q = q-p1;
  if (NbVectors[new_q] == NB_INIT_VALUE) { Vectors for length new_q not already computed }
    ComputeAutocorrelations( new_q );
  endif
  
  lim = NbVectorsQ + NbVectors[new_q];
  resize VectorsQ with NbVectors[new_q] new Bitarr*

  for( i = NbVectorsQ; i < lim; i++ ) // i index for VectorsQ

    ivnq =  i - NbVectorsQ   //{ ivnq index for Vectors vector_index_new_q }
     create and clear VectorsQ[i];
     set bits 0, p1 of VectorsQ[i] to 1

     for ( j=1; j < new_q; j++ ) 

       if ( bit j of Vectors[new_q][ivnq] to 1 )
         set bit p1+j of VectorsQ[i] to 1
       endif
     endfor j
  endfor i
    
  NbVectorsQ +=  NbVectors[new_q];  // update Nb of vectors of length q

endfor p1

**********************************************************************
*                    END OF ALGORITHM
**********************************************************************
 */

#include <stdio.h>
#include <stdlib.h>
#include <memoire.h>

/*#define _DEBUG */
#define    INIT_NB_VECTORS -1
#define    MAX_SIZES 101
#define    MAX_VECTORS  35201

typedef  int *pint;

int 	NbVectors[MAX_SIZES];

pint    Vectors[MAX_SIZES][MAX_VECTORS];

void InitAutocorrelations (){

    int i, j;

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

    for (j=0; j < MAX_SIZES; )	{
	for (i = 0; i < MAX_VECTORS; )
	    Vectors[j][i++] = NULL;
	NbVectors[j++] = INIT_NB_VECTORS;
    }
    /* initialise for word sizes 0, 1 and 2 */
    NbVectors[0] = 0;    
    NbVectors[1] = 1;
    Vectors[1][0] = callocI(1); Vectors[1][0][0] = 1;
    NbVectors[2] = 2;
    Vectors[2][0] = callocI(2); Vectors[2][0][0] = 1;    Vectors[2][0][1] = 0;
    Vectors[2][1] = callocI(2); Vectors[2][1][0] = 1;    Vectors[2][1][1] = 1;

    return;
}

int ComputeAutocorrelations( int q ){


    int 
	i,     /* index for VectorsQ */
	ivnq,  /* index for Vectors vector_index_new_q } */
	p1 = 1,    /* first non-zero period  */
	j, k, 
	half_q,
	beta,  /* [(greatest multiple of p1) < q] minus p1 */
	new_q; /* new length of strings for which Autocorrelations are recursively computed 
		  new_q <= q/2 + 1;     */

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

    /*     initialise for p1 = 0 and 1; */
    Vectors[q][0] = callocI(q); Vectors[q][0][0] = 1;
    Vectors[q][1] = callocI(q); for (j=0; j < q; ) Vectors[q][1][j++] = 1;
    i = NbVectors[q] = 2;
    half_q = q/2;


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

	if ( q % p1 == 0 )
	    beta = q - p1;
	else
	    beta = q - p1 - (q % p1);

	new_q = q - beta;
	if (NbVectors[new_q] == INIT_NB_VECTORS) /* { Vectors for length new_q not already computed } */
	    ComputeAutocorrelations( new_q );

	/*  Separate in two subcases (new_q == p1) and (new_q > p1) */ 
	if ( new_q == p1 ) {

	    /*  create a new vector for the nested vector 100...0 */
	    Vectors[q][i] = callocI(q); Vectors[q][i][0] = 1;
	    for ( k = p1; k <= beta; k += p1 ) Vectors[q][i][k] = 1;  
	    i++;

	    for ( ivnq = 2; ivnq < NbVectors[new_q]; ivnq++ ){

//#ifdef _DEBUG
       fprintf(stderr,"%d %d: ", new_q, p1);
       for(k=0;k < new_q; k++)fprintf(stderr,"%d ", Vectors[new_q][ivnq][k]);
       fprintf(stderr,"\n");
//#endif
		/*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 */
		for ( j=2; (Vectors[new_q][ivnq][j] == 0); j++);

		/* INVARIANT  ( j < new_q=p1 ) */
		if ( (p1 % j) == 0) continue;
fprintf(stderr,"here %d %d\n", j, p1);
		/* INV current vector of length new_q is valid */
		
		/* create and prepare Vectors[q][i] */
		Vectors[q][i] = callocI(q); Vectors[q][i][0] = 1;
		for ( k = p1; k <= beta; k += p1 ) Vectors[q][i][k] = 1;  

		Vectors[q][i][beta+j] = 1;  /* update for the current period j */
		
		/*  search for bit to 1  and update the new Vectors[q][i]  */
		for ( j++; j < new_q; j++ ) 
		    if ( Vectors[new_q][ivnq][j] ) Vectors[q][i][beta+j] = 1;  

		i++;     /* increment the index of Vectors[q] */
	    }/*   endfor ivnq */
	}
	else { /* new_q > p1 */

	    for ( ivnq = 2; ivnq < NbVectors[new_q]; ivnq++ ){

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

		/*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 */
		for ( j=2; (Vectors[new_q][ivnq][j] == 0); j++);

		if (j > p1) break;
		else 
		    if ((j<p1) && ((p1 % j) == 0)) continue;

		/* INVARIANT  ( j < new_q ) i.e. current vector of length new_q is valid */
		
		/* create and prepare Vectors[q][i] */
		Vectors[q][i] = callocI(q); Vectors[q][i][0] = 1;
		for ( k = p1; k <= beta; k += p1 ) Vectors[q][i][k] = 1;  
		/* AFTER: k == beta+p1 */
		Vectors[q][i][k] = 1;

		Vectors[q][i][beta+j] = 1;  /* update for the current period j */
		
		/*  search for bit to 1  and update the new Vectors[q][i]  */
		for ( j++; j < new_q; j++ ) 
		    if ( Vectors[new_q][ivnq][j] ) Vectors[q][i][beta+j] = 1;  

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


#ifdef _DEBUG
	fprintf(stderr,"p1 %d i %d\n", p1, i);
#endif
	fprintf(stdout,"NbVectors[ %5d %5d ] %6d\n",q,p1,i-NbVectors[q]);
	NbVectors[q] = i;  /* update the nb of vectors of length q */

    } /* endfor p1 */


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

	new_q = q-p1;
	if (NbVectors[new_q] == INIT_NB_VECTORS) /* { Vectors for length new_q not already computed } */
	    ComputeAutocorrelations( new_q );

	fprintf(stdout,"NbVectors[ %5d %5d ] %6d\n",q,p1,NbVectors[new_q]);

	for ( ivnq = 0; ivnq < NbVectors[new_q]; ivnq++, i++ ){ /* "map" a vector for q to every vector of new_q */

	    /* 	create and prepare Vectors[q][i] : bit 0 and p1 to 1*/
	    Vectors[q][i] = callocI(q); Vectors[q][i][0] = 1; Vectors[q][i][p1] = 1;

	    for ( j=1; j < new_q; j++ ) 
		if ( Vectors[new_q][ivnq][j] ) Vectors[q][i][p1+j] = 1;

	    /* endfor j */
	}/* endfor i */
    
	NbVectors[q] +=  NbVectors[new_q];  /*  update Nb of vectors of length q */
    } /* endfor p1 */

#ifdef _DEBUG
    fprintf(stderr,"End ComputeAutocorrelations %d %d\n", q, NbVectors[q]);
#endif
    fprintf(stdout,"\nNbVectors[ %5d ]      %6d\n",q,NbVectors[q]);

    return NbVectors[q];
   
}

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

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

    if (( argc < 2 ) || (atoi(argv[1]) < 1) || (atoi(argv[1]) > MAX_SIZES))  {
	fprintf(stderr,"usage: %s word_size(<%d) [-all] [-save_q | -save_all]\n", argv[0], MAX_SIZES);
	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 (argc >= 3) { 
	if (!strcmp(argv[2],StringOptionAll)) {
	    OptionAll = 1;

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


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

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

    if (OptionAll) {
	printf("Compute ALL sets of autocorrelation vector for word size q GOING FROM 1 TO %d\n",q);
	for ( i = 3; i <= q; i++ ) /* for all sizes */
	    nb_vectors_q = ComputeAutocorrelations(i);
    }
    else{
	printf("Compute the set of autocorrelation vector for word size q == %d\n",q);
	nb_vectors_q = ComputeAutocorrelations(q);
    }

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



    strcpy(fname,"AutocorrelationVectors.");
    if ( OptionSaveAll ) i = 3;
    else /*  OptionSaveQ  save only q to file */
	i = q;

    /* output results to file */
    for ( ; 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][k++]);
		fprintf(fout,"\n");
	    }
	    fclose(fout);
	}
    }
    return 0;
}