bch_erasures.c

error correction code

  printf("\n");

  if (syn_error) 
    {
    //
    // Compute the error location polynomial with the Euclidean algorithm
    // 

    for (i=0; i<=d; i++) {
      r[0][i] = 0;
      r[1][i] = 0;
      b[0][i] = 0;
      b[1][i] = 0;
      qt[i] = 0;
     }

  b[1][0] = 1; degb[0] = 0; degb[1] = 0;

  r[0][d] = 1; // x^{2t+1}
  degr[0] = d;

  for (i=0; i<=t2; i++)
    {
      if (s[i]!=-1) {
        r[1][i] = alpha_to[s[i]];
        degr[1] = i;
      }
      else
        r[1][i] = 0;
    }

  j = 1;

  if( (degr[0]-degr[1]) < t ) {

  do {

    j++;

printf("\n************************  j=%3d\n", j);
    // ----------------------------------------
    // Apply long division to r[j-2] and r[j-1]
    // ----------------------------------------

    // Clean r[j]
    for (i=0; i<=d; i++) r[j][i] = 0;

    for (i=0;i<=degr[j-2];i++) 
      r[j][i] = r[j-2][i]; 
    degr[j] = degr[j-2];

    temp = degr[j-2]-degr[j-1];
    for (i=temp; i>=0; i--) {
      u = degr[j-1]+i;
      if (degr[j] == u)
        {
        if ( r[j][degr[j]] && r[j-1][degr[j-1]] )
          qt[i] = alpha_to[(index_of[r[j][degr[j]]]
                                       +n-index_of[r[j-1][degr[j-1]]])%n];

//printf("r[j][degr[j]]] = %d,  r[j-1][degr[j-1]] = %d\n",
//index_of[r[j][degr[j]]], index_of[r[j-1][degr[j-1]]]);
printf("\nqt[%d]=%d\n", i, index_of[qt[i]]);

        for (u=0; u<=d; u++) aux[u] = 0;

        temp = degr[j-1];
        for (u=0; u<=temp; u++)
          if ( qt[i] && r[j-1][u] )
            aux[u+i] = alpha_to[(index_of[qt[i]]+index_of[r[j-1][u]])%n];
          else
            aux[u+i] = 0;

printf("r    = ");
for (u=0; u<=degr[j]; u++) printf("%4d ", index_of[r[j][u]]);
printf("\n");
printf("aux  = ");
for (u=0; u<=degr[j-1]+i; u++) printf("%4d ", index_of[aux[u]]);
printf("\n");

        for (u=0; u<=degr[j]; u++)
          r[j][u] ^= aux[u];
        u = d;
        while ( !r[j][u] && (u>0)) u--;
        degr[j] = u;
        }
      else
        qt[i] = 0;

printf("r    = ");
for (u=0; u<=degr[j]; u++) printf("%4d ", index_of[r[j][u]]);
printf("\n");

      }

printf("\nqt = ",j);
temp = degr[j-2]-degr[j-1];
for (i=0; i<=temp; i++) printf("%4d ", index_of[qt[i]]);
printf("\nr = ");
for (i=0; i<=degr[j]; i++) printf("%4d ", index_of[r[j][i]]);
printf("\nb = ");

    // Compute b(x)

    for (i=0; i<=d; i++) 
      aux[i] = 0; 

    temp = degr[j-2]-degr[j-1];
    for (i=0; i<=temp; i++)
      for (u=0; u<=degb[j-1]; u++)
        if ( qt[i] && b[j-1][u] )
          aux[i+u] ^= alpha_to[(index_of[qt[i]]+index_of[b[j-1][u]])%n];

    for (i=0; i<=d; i++) 
      b[j][i] = b[j-2][i] ^ aux[i];

    u = d;
    while ( !b[j][u] && (u>0) ) u--;
    degb[j] = u;

for (i=0; i<=degb[j]; i++) printf("%4d ", index_of[b[j][i]]);
printf("\n");

  } while (degr[j]>t); 

}

  u=1;
  temp = degb[j];
  // Normalize error locator polynomial
  for (i=0;i<=temp;i++) {
    elp[u][i] = alpha_to[(index_of[b[j][i]]-index_of[b[j][0]]+n)%n];
    }
  l[u] = temp;

  if (l[u] <= t) {
    /* put elp into index form */
    for (i = 0; i <= l[u]; i++)
      elp[u][i] = index_of[elp[u][i]];

    printf("sigma(x) = ");
    for (i = 0; i <= l[u]; i++)
      printf("%3d ", elp[u][i]);
    printf("\n");
    printf("Roots: ");

    /* Chien search: find roots of the error location polynomial */
    for (i = 1; i <= l[u]; i++)
      reg[i] = elp[u][i];
    count = 0;
    for (i = 1; i <= n; i++) {
      q = 1;
      for (j = 1; j <= l[u]; j++)
        if (reg[j] != -1) {
          reg[j] = (reg[j] + j) % n;
          q ^= alpha_to[reg[j]];
        }
      if (!q) {
        root[count] = i;
        loc[count] = n - i;
        count++;
        printf("%3d ", n - i);
      }
    }
  printf("\n");

  if (count == l[u])	
    /* no. roots = degree of elp hence <= t errors */
    for (i = 0; i < l[u]; i++)
      recd[loc[i]] ^= 1;
  else
    printf("Incomplete decoding: errors detected\n");
  }
 }
}



main()
{
  int i, j, OK;
  int buffer[32768], recd0[32768];
  int dist0, dist1;


	read_p();               /* Read m */
	generate_gf();          /* Construct the Galois Field GF(2**m) */
	gen_poly();             /* Compute the generator polynomial of BCH code */


	/* Randomly generate DATA */
	seed = 131073;
	srandom(seed);
	for (i = 0; i < k; i++)
		data[i] = ( random() & 65536 ) >> 16;

	encode_bch();           /* encode data */

	/*
	 * recd[] are the coefficients of c(x) = x**(length-k)*data(x) + b(x)
	 */
	for (i = 0; i < length - k; i++)
		recd[i] = bb[i];
	for (i = 0; i < k; i++)
		recd[i + length - k] = data[i];
	printf("Code polynomial:\nc(x) = ");
	for (i = 0; i < length; i++) {
		printf("%1d", recd[i]);
		if (i && ((i % 50) == 0))
			printf("\n");
	}
	printf("\n");



    printf("Enter the number of earsures:\n"); scanf("%d", &numera);
    printf("Enter erasure locations (integers between");
    printf(" 0 and %d): ", length-1);
    for (i = 0; i < numera; i++)
      scanf("%d", &erapos[i]);

    printf("Enter the number of errors:\n"); scanf("%d", &numerr);	
    printf("Enter error locations (integers between");
    printf(" 0 and %d): ", length-1);
    for (i = 0; i < numerr; i++)
      scanf("%d", &errpos[i]);
    if (numerr)
      for (i = 0; i < numerr; i++)
        recd[errpos[i]] ^= 1;

	printf("r(x) = ");
	for (i = 0; i < length; i++) {
		printf("%1d", recd[i]);
		if (i && ((i % 50) == 0))
			printf("\n");
	}
	printf("\n");




    for (i = 0; i < length; i++)  buffer[i] = recd[i];
    // Substitute ZEROS in the erased positions and decode
    for (i = 0; i < numera; i++)  recd[erapos[i]] = 0;
    decode_bch(); 
    for (i = 0; i < length; i++)  recd0[i] = recd[i];

    for (i = 0; i < length; i++)  recd[i] = buffer[i];
    // Substitute ONES in the erased positions and decode
    for (i = 0; i < numera; i++)  recd[erapos[i]] = 1;
    decode_bch(); 

    // Choose closest of the two codewords to the received word

    dist0 = 0;
    for (i = 0; i < length; i++) {
      // Compare only not erased postions
      OK = 1;
      for (j=0; j<numera; j++) if (i == erapos[j]) OK = 0;
      if (OK)
        if (buffer[i] != recd0[i]) dist0++;
      }

    dist1 = 0;
    for (i = 0; i < length; i++) {
      // Compare only not erased postions
      OK = 1;
      for (j=0; j<numera; j++) if (i == erapos[j]) OK = 0;
      if (OK)
        if (buffer[i] != recd[i]) dist1++;
      }

    if (dist0 < dist1)
      for (i = 0; i < length; i++) recd[i] = recd0[i];
    



	printf("Results:\n");
	printf("original data  = ");
	for (i = 0; i < k; i++) {
		printf("%1d", data[i]);
		if (i && ((i % 50) == 0))
			printf("\n");
	}
	printf("\nrecovered data = ");
	for (i = length - k; i < length; i++) {
		printf("%1d", recd[i]);
		if ((i-length+k) && (((i-length+k) % 50) == 0))
			printf("\n");
	}
	printf("\n");


	/*
	 * DECODING ERRORS? Comparing only the data portion
	 */
	for (i = length - k; i < length; i++)
		if (data[i - length + k] != recd[i])
			decerror++;
	if (decerror)
	   printf("There were %d decoding errors in message positions\n", decerror);
	else
	   printf("Succesful decoding\n");
}