/* CSC 310 ASSIGNMENT #4 - TEST DECODING METHODS FOR PRODUCT OF HAMMING CODES. 
 
   Radford M. Neal, 2002 */

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


/* ADJUSTABLE CONSTANTS. */

#define N_sims 10000	/* Number of transmissions done in test */
#define E_prob 0.08	/* Channel error probablity */


/* DECODING METHODS.  The N_METHODS constant is the number of methods. */

enum { ROWCOL, ROWCOLROWCOL, ENUMERATE, N_METHODS };


/* FIND CHECK BITS FOR [7,4] HAMMING CODE. */

void encode_hamming 
( char x1, char x2, char x3, char x4, char *c1, char *c2, char *c3
)
{
  *c1 = x2^x3^x4;
  *c2 = x1^x3^x4;
  *c3 = x1^x2^x4;
}


/* DECODE DATA FROM [7,4] HAMMING CODE. */

void decode_hamming 
( char *x1, char *x2, char *x3, char *x4, char *c1, char *c2, char *c3
)
{ 
  char s1, s2, s3;

  /* Compute the syndrome bits. */

  s1 = *x2 ^ *x3 ^ *x4 ^ *c1;
  s2 = *x1 ^ *x3 ^ *x4 ^ *c2;
  s3 = *x1 ^ *x2 ^ *x4 ^ *c3;

  /* Use syndrome to determine which bit to flip (if any). */

  switch (4*s1 + 2*s2 + s3)
  {
    case 0: break;
    case 1: *c3 ^= 1; break;
    case 2: *c2 ^= 1; break;
    case 3: *x1 ^= 1; break;
    case 4: *c1 ^= 1; break;
    case 5: *x2 ^= 1; break;
    case 6: *x3 ^= 1; break;
    case 7: *x4 ^= 1; break;
  }
}


/* FIND CHECK BITS FOR PRODUCT OF TWO [7,4] HAMMING CODES. */

void encode_product 
( char d[49]
)
{   
  /* Find check bits for rows. */

  encode_hamming ( d[0], d[1], d[2], d[3], &d[4], &d[5], &d[6]);
  encode_hamming ( d[7], d[8], d[9],d[10],&d[11],&d[12],&d[13]);
  encode_hamming (d[14],d[15],d[16],d[17],&d[18],&d[19],&d[20]);
  encode_hamming (d[21],d[22],d[23],d[24],&d[25],&d[26],&d[27]);

  /* Find check bits for columns. */

  encode_hamming ( d[0], d[7], d[14], d[21], &d[28], &d[35], &d[42]);
  encode_hamming ( d[1], d[8], d[15], d[22], &d[29], &d[36], &d[43]);
  encode_hamming ( d[2], d[9], d[16], d[23], &d[30], &d[37], &d[44]);
  encode_hamming ( d[3],d[10], d[17], d[24], &d[31], &d[38], &d[45]);
  encode_hamming ( d[4],d[11], d[18], d[25], &d[32], &d[39], &d[46]);
  encode_hamming ( d[5],d[12], d[19], d[26], &d[33], &d[40], &d[47]);
  encode_hamming ( d[6],d[13], d[20], d[27], &d[34], &d[41], &d[48]);
}


/* APPLY HAMMING DECODING TO ROWS OF PRODUCT CODE. */

void decode_rows
( char d[49]
) 
{
  decode_hamming ( &d[0], &d[1], &d[2], &d[3], &d[4], &d[5], &d[6]);
  decode_hamming ( &d[7], &d[8], &d[9],&d[10],&d[11],&d[12],&d[13]);
  decode_hamming (&d[14],&d[15],&d[16],&d[17],&d[18],&d[19],&d[20]);
  decode_hamming (&d[21],&d[22],&d[23],&d[24],&d[25],&d[26],&d[27]);
  decode_hamming (&d[28],&d[29],&d[30],&d[31],&d[32],&d[33],&d[34]);
  decode_hamming (&d[35],&d[36],&d[37],&d[38],&d[39],&d[40],&d[41]);
  decode_hamming (&d[42],&d[43],&d[44],&d[45],&d[46],&d[47],&d[48]);
}


/* APPLY HAMMING DECODING TO COLUMNS OF PRODUCT CODE. */

void decode_cols
( char d[49]
) 
{
  decode_hamming ( &d[0], &d[7], &d[14], &d[21], &d[28], &d[35], &d[42]);
  decode_hamming ( &d[1], &d[8], &d[15], &d[22], &d[29], &d[36], &d[43]);
  decode_hamming ( &d[2], &d[9], &d[16], &d[23], &d[30], &d[37], &d[44]);
  decode_hamming ( &d[3],&d[10], &d[17], &d[24], &d[31], &d[38], &d[45]);
  decode_hamming ( &d[4],&d[11], &d[18], &d[25], &d[32], &d[39], &d[46]);
  decode_hamming ( &d[5],&d[12], &d[19], &d[26], &d[33], &d[40], &d[47]);
  decode_hamming ( &d[6],&d[13], &d[20], &d[27], &d[34], &d[41], &d[48]);
}


/* DECODE RECEIVED DATA FROM A PRODUCT CODE.  The parameters are the received
   data bits, the place to store the decoded data bits, and the decoding 
   method to use. */

void decode_product
( char r[49],
  char d[49],
  int method
)
{ 
  int i;

  switch (method)
  { 
    case ROWCOL:  /* decode rows, then columns */
    { 
      for (i = 0; i<49; i++) d[i] = r[i];
      decode_rows(d);
      decode_cols(d);
      break;
    }

    case ROWCOLROWCOL:  /* decode rows, columns, rows again, columns again */
    { 
      for (i = 0; i<49; i++) d[i] = r[i];
      decode_rows(d);
      decode_cols(d);
      decode_rows(d);
      decode_cols(d);
      break;
    }
 
    case ENUMERATE:  /* Decode optimally by checking all 2^16 codewords */
    {
      char e[49];
      int v, w, md;

      md = 50;

      for (v = 0; v<(1<<16); v++)
      {
        w = v<<1;

         e[0] = (w>>=1)&1; e[1] = (w>>=1)&1; e[2] = (w>>=1)&1; e[3] = (w>>=1)&1;
         e[7] = (w>>=1)&1; e[8] = (w>>=1)&1; e[9] = (w>>=1)&1;e[10] = (w>>=1)&1;
        e[14] = (w>>=1)&1;e[15] = (w>>=1)&1;e[16] = (w>>=1)&1;e[17] = (w>>=1)&1;
        e[21] = (w>>=1)&1;e[22] = (w>>=1)&1;e[23] = (w>>=1)&1;e[24] = (w>>=1)&1;
        
        encode_product(e);

        if (hamming_distance(e,r)<md) 
        { md = hamming_distance(e,r);
          for (i = 0; i<49; i++) d[i] = e[i];
        }
      }
    }
  }
}


/* FIND THE HAMMING DISTANCE BETWEEN BLOCKS. */

int hamming_distance
( char d1[49],
  char d2[49]
)
{
  int s, i;
  s = 0;
  for (i = 0; i<49; i++)
  { s += d1[i]^d2[i];
  }
  return s;
}


/* GENERATE A RANDOM BIT WITH 0 AND 1 EQUALLY PROBABLE. */

int rbit()
{ 
  return drand48()<0.5;
}


/* MAIN PROGRAM. */

main(void)
{ 
  static int tries[11], all_ok[11][N_METHODS], mess_ok[11][N_METHODS],
             total_all_ok[N_METHODS], total_mess_ok[N_METHODS];
  char t[49], r[49], d[49];
  int errs, i, j, m;

  /* Simulate N_sims transmissions. */

  for (i = 0; i<N_sims; i++)
  {
    /* Generate a random message. */

     t[0] = rbit();  t[1] = rbit();  t[2] = rbit();  t[3] = rbit();
     t[7] = rbit();  t[8] = rbit();  t[9] = rbit(); t[10] = rbit();
    t[14] = rbit(); t[15] = rbit(); t[16] = rbit(); t[17] = rbit();
    t[21] = rbit(); t[22] = rbit(); t[23] = rbit(); t[24] = rbit();

    /* Encode it with the product code. */

    encode_product(t);

    /* Add noise to produce received data. */

    for (j = 0; j<49; j++) 
    { r[j] = t[j] ^ (drand48()<E_prob);
    }

    errs = hamming_distance(t,r);
    if (errs>10) errs = 10;

    tries[errs] += 1;

    /* Decode by each method. */

    for (m = 0; m<N_METHODS; m++)
    { 
      int ok, mok;

      decode_product(r,d,m);

      ok = hamming_distance(t,d)==0;
      all_ok[errs][m] += ok;
      total_all_ok[m] += ok;

      mok =  t[0]==d[0]  &&  t[1]==d[1]  &&  t[2]==d[2]  &&  t[3]==d[3]  &&
             t[7]==d[7]  &&  t[8]==d[8]  &&  t[9]==d[9]  && t[10]==d[10] &&
            t[14]==d[14] && t[15]==d[15] && t[16]==d[16] && t[17]==d[17] &&
            t[21]==d[21] && t[22]==d[22] && t[23]==d[23] && t[24]==d[24];
      mess_ok[errs][m] += mok;
      total_mess_ok[m] += mok;
    }
  }

  /* Print the table of results. */

  printf(
"                 Results:  all-ok-count/fraction  message-ok-count/fraction");
  printf("\n\n");
  printf("errs  tries |");
  for (m = 1; m<=N_METHODS; m++) printf("       Method %d      |",m);
  printf("\n");

  for (errs = 0; errs<=10; errs++)
  { 
    if (errs==10) printf(">=10 ");
    else printf("%4d ",errs);

    printf(" %5d |",tries[errs]);
 
    for (m = 0; m<N_METHODS; m++)
    { if (tries[errs]!=0)
      { printf("%4d/%.3f %4d/%.3f|", 
               all_ok[errs][m], (double)all_ok[errs][m]/tries[errs],
               mess_ok[errs][m], (double)mess_ok[errs][m]/tries[errs]);
      }
      else
      { printf("                     |");
      }
    }
    
    printf("\n");
  }

  printf("total");
  printf(" %5d |",N_sims);
  
  for (m = 0; m<N_METHODS; m++)
  { printf("%4d/%.3f %4d/%.3f|", 
           total_all_ok[m], (double)total_all_ok[m]/N_sims,
           total_mess_ok[m], (double)total_mess_ok[m]/N_sims);
  }

  printf("\n");

  exit(0);
}