genmat_ebch2.c

error correction code

/*
 * File:        genmat_ebch2.c
 * Description: Compute the generator polynomial of a cyclic BCH code and the
 *              generator matrix of its extended (2^m,k) code.
 * Input:       m and t, where m <= M_MAX (defined below).
 * Output:      Generator polynomial (stdout) and file (matrix.out) containing
 *              the generator matrix.
 * Revised:     February 5, 1997.
 * Author:      Robert Morelos-Zaragoza. Copyright 1997. All rights reserved.
 *
 */
// ------------------------------------------------------------------------
// This program is complementary material for the book:
//
// R.H. Morelos-Zaragoza, The Art of Error Correcting Coding, Wiley, 2002.
//
// ISBN 0471 49581 6
//
// This and other programs are available at http://the-art-of-ecc.com
//
// You may use this program for academic and personal purposes only. 
// If this program is used to perform simulations whose results are 
// published in a journal or book, please refer to the book above.
//
// The use of this program in a commercial product requires explicit
// written permission from the author. The author is not responsible or 
// liable for damage or loss that may be caused by the use of this program. 
//
// Copyright (c) 2002. Robert H. Morelos-Zaragoza. All rights reserved.
// ------------------------------------------------------------------------

#define M_MAX 9
#define K_MAX 511
#define N_MAX 512

#include <math.h>
#include <stdio.h>

int m, n_ext, n, length, k, t, d;
int p[10], alpha_to[1024], index_of[1024], g[1024], G[K_MAX][N_MAX];
FILE *fp;

void read_p()
/*
 *  Get precomputed coefficients p[] of p(x) the primitive polynomial
 *  used to compute the Galois field GF(2**m). Compute the code length.
 */
{
  int i;

  for (i=1; i<m; i++)
    p[i] = 0;
  p[0] = p[m] = 1;

  if (m == 2)			p[1] = 1;
  else if (m == 3)	p[1] = 1;
  else if (m == 4)	p[1] = 1;
  else if (m == 5)	p[2] = 1;
  else if (m == 6)	p[1] = 1;
  else if (m == 7)	p[1] = 1;
  else if (m == 8)	p[4] = p[5] = p[6] = 1;
  else if (m == 9)	p[4] = 1;
  else if (m == 10)	p[3] = 1;
  else if (m == 11)	p[2] = 1;
  else if (m == 12)	p[3] = p[4] = p[7] = 1;
  else if (m == 13)	p[1] = p[3] = p[4] = 1;
  else if (m == 14)	p[1] = p[11] = p[12] = 1;
  else if (m == 15)	p[1] = 1;
  else if (m == 16)	p[2] = p[3] = p[5] = 1;
  else if (m == 17)	p[3] = 1;
  else if (m == 18)	p[7] = 1;
  else if (m == 19)	p[1] = p[5] = p[6] = 1;
  else if (m == 20)	p[3] = 1;
  n = 1;
  for (i = 0; i <= m; i++)
    n *= 2;
  n = n / 2 - 1;
  length = n;
  n_ext = n+1;
}


void generate_gf()
/*
 * Generate GF(2**m) using a primitive polynomial p(X) with
 * coefficients in p[0]..p[m].
 *
 * Lookup tables:
 *   index->polynomial form: alpha_to[] contains j=alpha^i;
 *   polynomial form -> index form:	index_of[j=alpha^i] = i
 *
 * alpha=2 is the primitive element of GF(2**m) 
 */
{
  register int i, mask;

  mask = 1;
  alpha_to[m] = 0;
  for (i = 0; i < m; i++)
    {
      alpha_to[i] = mask;
      index_of[alpha_to[i]] = i;
      if (p[i] != 0)
	alpha_to[m] ^= mask;
      mask <<= 1;
    }
  index_of[alpha_to[m]] = m;
  mask >>= 1;
  for (i = m + 1; i < n; i++) 
    {
      if (alpha_to[i - 1] >= mask)
	alpha_to[i] = alpha_to[m] ^ ((alpha_to[i - 1] ^ mask) << 1);
      else
	alpha_to[i] = alpha_to[i - 1] << 1;
      index_of[alpha_to[i]] = i;
    }
  index_of[0] = -1;
}


void gen_poly()
/*
 * Compute the generator polynomial of a binary BCH code. Fist generate the
 * cycle sets modulo 2**m - 1, cycle[][] =  (i, 2*i, 4*i, ..., 2^l*i). Then
 * determine those cycle sets that contain integers in the set of (d-1)
 * consecutive integers {1..(d-1)}. The generator polynomial is calculated
 * as the product of linear factors of the form (x+alpha^i), for every i in
 * the above cycle sets. (Note: this is not very efficient.)
 */
{
  register int ii, jj, ll, kaux;
  register int test, aux, nocycles, root, noterms, rdncy;
  int cycle[1024][21], size[1024], min[1024], zeros[1024];
  
  /* Generate cycle sets modulo n, n = 2**m - 1 */
  cycle[0][0] = 0;
  size[0] = 1;
  cycle[1][0] = 1;
  size[1] = 1;
  jj = 1; /* cycle set index */
  if (m > 9)
    {
      printf("Computing cycle sets modulo %d\n", n);
      printf("(This may take some time)...\n");
    }
  do {
    /* Generate the jj-th cycle set */
    ii = 0;
    do {
      ii++;
      cycle[jj][ii] = (cycle[jj][ii - 1] * 2) % n;
      size[jj]++;
      aux = (cycle[jj][ii] * 2) % n;
    } while (aux != cycle[jj][0]);
    /* Next cycle set representative */
    ll = 0;
    do {
      ll++;
      test = 0;
      for (ii = 1; ((ii <= jj) && (!test)); ii++)	
	/* Examine previous cycle sets */
	for (kaux = 0; ((kaux < size[ii]) && (!test)); kaux++)
	  if (ll == cycle[ii][kaux])
	    test = 1;
    } while ((test) && (ll < (n - 1)));
    if (!(test)) {
      jj++;	/* next cycle set index */
      cycle[jj][0] = ll;
      size[jj] = 1;
    }
  } while (ll < (n - 1));
  nocycles = jj;		/* number of cycle sets modulo n */
  
  if (!t)
    { /* Parity check code */
      d = 1;
      rdncy = 1;
      k = length;
      g[0] = g[1] = 1;
    }
  else
    { /* Non trivial BCH code */
      d = 2 * t + 1;
      
      /* Search for roots 1, 2, ..., d-1 in cycle sets */
      kaux = 0;
      rdncy = 0;
      for (ii = 1; ii <= nocycles; ii++)
	{
	  min[kaux] = 0;
	  test = 0;
	  for (jj = 0; ((jj < size[ii]) && (!test)); jj++)
	    for (root = 1; ((root < d) && (!test)); root++)
	      if (root == cycle[ii][jj])
		{
		  test = 1;
		  min[kaux] = ii;
		}
	  if (min[kaux]) 
	    {
	      rdncy += size[min[kaux]];
	      kaux++;
	    }
	}
      
      noterms = kaux;
      kaux = 1;
      for (ii = 0; ii < noterms; ii++)
	for (jj = 0; jj < size[min[ii]]; jj++)
	  {
	    zeros[kaux] = cycle[min[ii]][jj];
	    kaux++;
	  }
      
      k = length - rdncy;
      
      /* Compute the generator polynomial */
      g[0] = alpha_to[zeros[1]];
      g[1] = 1;		/* g(x) = (X + zeros[1]) initially */
      for (ii = 2; ii <= rdncy; ii++)
	{
	  g[ii] = 1;
	  for (jj = ii - 1; jj > 0; jj--)
	    if (g[jj] != 0)
	      g[jj] = g[jj - 1] ^ alpha_to[(index_of[g[jj]] + zeros[ii]) % n];
	    else
	      g[jj] = g[jj - 1];
	  g[0] = alpha_to[(index_of[g[0]] + zeros[ii]) % n];
	}
    }
}

void gen_matrix()
{
	
/* 
 * Given the generator polynomial of a binary cyclic code, calculate
 * the generator matrix of its extended code, in systematic form.
 */

  register int i,j,row,parity; 
  
  int w = 0; /* Hamming weight of generator polynomial */
  
  for (i=0; i<n_ext; i++)
    G[0][i] = 0;
  
  /* First row */
  row=0;
  
  while (row<k)
    {
      
      if (row)
	{
	  /* shift to right */
	  for (i=0;i<n;i++)
	    G[row][i+1] = G[row-1][i];
	}
      else
	{
	  for (i=0;i<n_ext-k;i++)
	    G[row][i] = g[i];
	  for (i=n_ext-k+1; i<n; i++)
	    G[row][i] = 0;
	}
      
      /* Parity */
      parity = 0;
      for (i=0;i<n; i++)
        parity ^= G[row][i];
      G[row][n_ext-1] = parity;
      
      row++;
      
    }
  
  /* Put in systematic form */
  
  row = 0;
  
  while (row<k)
    {
      
      for (i=row+1; i<k; i++)
        if (G[row][i])
          for (j=0; j<n_ext; j++)
            G[row][j] ^= G[i][j];
      
      /* Print row */
      for (i=0; i<n_ext; i++)
        fprintf(fp,"%1d ", G[row][i]);
      fprintf(fp,"\n");
      
      row++;
    }
}

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

  if (!(argc==4))
    {
      printf("USAGE:\na.out m t filename\n");
      printf("where 2<m<20 and 0<=t<2^{m-2}\n");
      exit(0);
    }

  m = atoi(argv[1]);
  t = atoi(argv[2]);
  strcpy(filename,argv[3]);

  if (m>M_MAX)
    {
      printf("ERROR:\nvalue of m too large\n");
      exit(0);
    }
  
  read_p();
  generate_gf();
  gen_poly();
  fp = fopen(filename,"w");
  gen_matrix();
  fclose(fp);
  printf("Finished eBCH(%d,%d,%d) code\n", n_ext, k, d+1);
}