/*!
 **************************************************************************
 \Firm:    QUALCOMM    
 \File:     b-m.cpp
 \Project   Berlekamp-Massey algorithm
 \Function  The following code tests an implementation of
 \Function  the Berlekamp-Massey algorithm for finding the
 \Function  linear complexity of a finite binary sequence.
 \Function  This is useful in analyzing linear feedback
 \Function  shift registers (LFSRs).see Menezes et al, P.201 
 \Version   $Id: b-m.cpp, v 1.4 09/02/02 $
 **************************************************************************
 **************************************************************************
 \History
  02.09.2002  V0.0- V1.3 by GGR at QUALCOMM; Algorithm itself, test enviroment and printPolynom function
  07.01.2002  V1.4 Ported to C++ by Heiko Mahr, Changed s from unsigned char to int.
              Renamed variables according german literature,Added most of the comments
              Made BMA dynamical, and prgrammed the BMA to a parametrable function
  
  For variablenames i use hungarian notation.
  So the leading smallcharacters say something of the type of the variable
  p : Pointer
  a:  Array
  i:  integer
  f:  float
  pai is a Pointer on a Array of Integer. aso.
 **************************************************************************
*/

/**************************************************************************
  includes
 **************************************************************************/
#include <stdio.h>
#include <string.h>
#include <memory.h>
#define MAX 100000

//0000 1001 0111 0001 0011  // testpolynom for argc[2] 

/**************************************************************************
  global variables
 **************************************************************************/
unsigned char	inbuf[MAX / 8]; // byte-binary input buffer 

/*!
 ***************************************************************************
 \FctName          printPolynom
 \FctDescription   Prints Polynom to Screen
 \Parameter         char *pcString
 \Description       Pointer to C-String
 \Parameter         int  *piPolynom
 \Description       Pointer on Polynom to print 
 \Parameter         int iL
 \Description       Length of Polynom
 \Parameter         bool bOptions
 \Description       true: print as verbose, else Print as Polynom  
 ***************************************************************************
*/
void printPolynom(char *pcString,int *piPolynom, int iL, bool bOptions)
/*!
 ***************************************************************************
 \FctName          BMA
 \FctDescription   Berlekamp-Massey Algorithm
 \Parameter         int *paiBitsequenceS
 \Description       paiBitsequenceS [0..iLenghtBitSeq] is the Input Sequence 
 \Parameter         int iLenghtBitSeq
 \Description       Maximum Length of Bitseq 
 \Parameter         bool bOption
 \Description       true: print as verbose, else Print as Polynom  
 \Parameter         int *paiFeedbackpolyC
 \Description       Pointer on Array with Feedbackpolynom of the LFSR
 ***************************************************************************
*/
void  BMA(int *paiBitsequenceS,
          int iLenghtBitSeq,
          bool bOption,
          int *paiFeedbackpolyC);

void main ( int argc, char * argv[])
{
  bool bBMCOptions=false;// for output of the sequence as verbose or polynom
  int  aiS[MAX];
  int  iN;		// current / total number of input bits 
//Read In first Main Parameter
  while (argc >= 2 && argv[1][0] == '-')
  {
	  if (strcmp(argv[1], "-v") == 0)//Display verbose 
    {
      bBMCOptions=true;++argv; --argc;
    }
	  else if (strcmp(argv[1], "-p") == 0)//Display as polynom 
    {
	    bBMCOptions=false;  ++argv; --argc;
    }
  }

  if (argc > 2) // if more as two arguments have been given over
  {
	  fprintf(stderr, "Usage: %s [-v|-p] [bits]\n  (or binary on stdin)\n");	  
  }
  else if (argc == 2) //only two arguments ?
  {
	  iN = strlen(argv[1]);
	  for (int i = 0; i < iN; ++i)// write S(i) 
	    argv[1][i]== '0' ? aiS[i]=0 : aiS[i]=1;
  }
  else// less than two arguments ?
  {
    //This maybe wontīt work anymore
    // because s was first of all implemented as array of unsigned char
    // and itīs now implemented as array of int
	  iN = fread(inbuf, 1, MAX/8, stdin);
	  if (iN <= 0)
	    fprintf(stderr, "%s: No Input to BMA");
    // ReadIn from st
	  for (int i = 0; i < iN; ++i)
      for (int j = 0; j < 8; ++j)
        aiS[8*i + j] = (inbuf[i] >> (7 - j)) & 0x1;
    iN *= 8;
  }
/*Programmpart of HM*/
  int *paiC= new int [iN+1];//Max Length of C is N+1
  memset(C,0,(iN+1)*sizeof(int));
  BMA(aiS,iN,bBMCOptions,paiC);
  delete [] paiC; 
}


void  BMA(int *paiBitsequenceS, int iLengthBitSeq,bool bOption,int *paiFeedbackpolyC)
{
  // paiBitsequenceS [0..iLenghtBitSeq] is the Input Sequence 
  // paiFeedbackpolyC [0...iLenghtBitSeq+1] is the Feedbackpolynom of the LFSR
  int	k=0;        // Iterationcounter ( and Number of Bits which are read in)
  int	L= 0;		    // actual linear Complexity L of the Input Sequence
  bool warned=false;   // For detection of jumps in the linear Complexity
  int	m=-1;       // Number of the last Iteration where C(x) has been changed
  
  int	d=paiBitsequenceS[k];//Diskrepanz            
  int	*B= new  int[iLengthBitSeq+1];//in Literature sometimes named C*(X)
  int *T= new  int[iLengthBitSeq+1];//in  Literature sometimes T(x)
  // INIT
  memset(T,0,(iLengthBitSeq+1)*sizeof(int));
  memset(B,0,(iLengthBitSeq+1)*sizeof(int));
  paiFeedbackpolyC[0]=B[0] = 1;

/*Loop from 0 to Length of Inputsequnce
  In every Iteration :
  -Get 1 Bit more from the Inputsequence
  -Prove if the sequence could be created by a LFSR with the feedbackpolynom C(x)
  -If condition not matched: Change C(x)
*/
  while (k < iLengthBitSeq) 
  {
	  // Calculate the discrepancy d
	  d = paiBitsequenceS[k];
	  for (int i = 1; i <= L; ++i)
      d ^= *(paiFeedbackpolyC+i) & *(paiBitsequenceS+(k - i));
	  
    if (d==1)// if d is 1 then C(x) has to be corrected
    {
      // "save" old C(x) to T(x)
      memcpy(T,paiFeedbackpolyC,(k+1)*sizeof(int));
      // build new Cj(X) with old B(X)
	    for (int i = k-m, j = 0; i <= k+1; ++i, ++j)
        paiFeedbackpolyC[i] ^= B[j];
		  if (L <= k/2) // if 2*L <= k then linear complexity raises
      {
		    //Jump in linear komplexity, Test if is unusual 
		    if ((k/2 - L) > 6 && L > 2 && !warned)
        {
		      printf("Large jump (%d) in linear complexity detected\n",k - 2*L);
		      warned = true;
		    }
		    L = k + 1 - L;
		    m = k;
        memcpy(B,T,(k+1)*sizeof(int));		    
	    }
      // Print T,C,B L,k m to Screen
		  printf("k=%d s[%d]=%d d=%d L=%d m=%d\n", k, k-1, paiBitsequenceS[k-1], d, L, m);
		  printPolynom("T=", T,L, bOption);
		  printPolynom("C=", paiFeedbackpolyC,L, bOption);
		  printPolynom("B=", B,L,bOption);
	  }
	++k;
  }// END WHILE 

  // count weight of polynomial -- don't include last bit 
  int j=0;
  for (int i=0; i < L; ++i)
	  j += paiFeedbackpolyC[i];
  printf("L=%d W=%d\n", L, j);
  printPolynom("", paiFeedbackpolyC,L,bOption);
  
  if (L > ( (iLengthBitSeq/2) -4))
  {
	  printf("Sequence is probably nonlinear ");
	  if (iLengthBitSeq == MAX)
	    printf(" up to MAX (%d) tested \n", MAX);
	  printf("\n");
  }                  
  delete [] B;
  delete [] T;
}

//Print Routine
void printPolynom(char *pcString,  int *piPolynom, int iL, bool bOptions)
{
  register int	i;
  int	  first = 1;

  if (bOptions==true)
  {
	  printf("%s", pcString);
	  for (i = 0; i <= iL; ++i)
	    printf("%d", piPolynom[i]);	
  }
  else 
  {
	  for (i  = 0; i <= iL; ++i) 
    {
	    if (piPolynom[i]) 
      {
	      if (first)
          first = 0;
	      else
		      printf(" + ");
      }
	    if (i == 0 && piPolynom[i])
		    printf("1");
	    else if (i == 1 && piPolynom[i])
		    printf("x");
	    else if (piPolynom[i])
		    printf("x^%d", i);
    }
  }
  printf("\n");
}