/* CSC 310 - Practical Assignment #2 Solution */

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


/**** PART A - SYNDROME DECODING ***/


/* READ A PARITY CHECK MATRIX. */

void read_pchk (int n_pchks, int cw_len, int pchk[])
{
  int i, j, b;

  if (n_pchks > sizeof(int)*8-1)
  { fprintf(stderr,"Too many rows in parity check matrix\n");
    exit(1);
  }

  if (cw_len > 256)
  { fprintf(stderr,"Too many columns in parity check matrix\n");
    exit(1);
  }

  for (j = 0; j<cw_len; j++) pchk[j] = 0;

  for (i = 0; i<n_pchks; i++)
  { 
    for (j = 0; j<cw_len; j++)
    { 
      if (scanf("%d",&b)!=1 || b!=0 && b!=1)
      { fprintf(stderr,
           "Bad input in row %d col %d of parity check matrix\n",i+1,j+1);
        exit(1);
      }

      pchk[j] |= b<<(n_pchks-i-1);
    }
  }
}

/* Test program - change name to "main" and remove other "main" to enable. */

int test_read_pchk/*main*/(void)
{
  int n_pchks, cw_len, j;

  scanf("%d",&n_pchks);
  scanf("%d",&cw_len);

  int pchk[cw_len];

  read_pchk(n_pchks,cw_len,pchk);
  for (j = 0; j<cw_len; j++) printf("%08x\n",pchk[j]);
  return 0;
}


/* ERROR PATTERN GENERATION. */

void first_err_pat (int n_errors, int err_pos[])
{
  int i;
  for (i = 0; i<n_errors; i++) err_pos[i] = i;
}

int next_err_pat (int cw_len, int n_errors, int err_pos[])
{
  int k;
  for (k = 1; err_pos[n_errors-k] == cw_len-k; k++)
  { if (k==n_errors) return 0;
  }

  err_pos[n_errors-k] += 1;

  for (k = k-1; k > 0; k--)
  { err_pos[n_errors-k] = err_pos[n_errors-k-1] + 1;
  }

  return 1;
}

/* Test program - change name to "main" and remove other "main" to enable. */

int test_err_pat/*main*/(void)
{
  int err_pos[3];
  first_err_pat(3,err_pos);
  do
  { int i;
    for (i = 0; i<3; i++) printf("%d ",err_pos[i]);
    printf("\n");
  } while (next_err_pat(6,3,err_pos));
  return 0;
}


/* FILL IN A SYNDROME DECODING TABLE */

void make_syndrome_table (int n_pchks, int cw_len, int pchk[], 
  int syn_err_count[], unsigned char syn_err_pos[][cw_len])
{ 
  unsigned int table_size = 1U<<n_pchks; /* Size of syndrome decoding table */

  unsigned int z;
  int err_pos[cw_len];
  unsigned filled;
  int n_errors;
  int j;

  /* Initialize all syndrom table entries to -1, indicating not filled yet. */

  for (z = 0; z<table_size; z++) syn_err_count[z] = -1;

  /* Set up entry for syndrome of zero; initialize for number of errors = 1. */

  syn_err_count[0] = 0;
  filled = 1;
  n_errors = 1;
  first_err_pat(n_errors,err_pos);

  /* Look at error patterns by non-decreasing weight until all looked at
     or all syndrome table entries filled. */
 
  while (filled<table_size)
  { 
    z = find_syndrome(pchk,n_errors,err_pos);

    if (syn_err_count[z]==-1)
    { filled += 1;
      syn_err_count[z] = n_errors;
      if (syn_err_pos!=0)
      { for (j = 0; j<n_errors; j++) syn_err_pos[z][j] = err_pos[j];
      }
    }

    if (!next_err_pat(cw_len,n_errors,err_pos))
    { n_errors += 1;
      if (n_errors>cw_len) break;
      first_err_pat(n_errors,err_pos);
    }
  } 
}

/* Test program - change name to "main" and remove other "main" to enable. */

int test_syndrom_table/*main*/(void)
{
  int n_pchks, cw_len;
  unsigned z;
  int j;

  scanf("%d",&n_pchks);
  scanf("%d",&cw_len);

  int pchk[cw_len];

  read_pchk(n_pchks,cw_len,pchk);

  int syn_err_count[1U<<n_pchks];
  unsigned char syn_err_pos[1U<<n_pchks][cw_len];

  make_syndrome_table (n_pchks, cw_len, pchk, syn_err_count, syn_err_pos);

  for (z = 0; z<(1U<<n_pchks); z++)
  { printf("%8x %d :",z,syn_err_count[z]);
    for (j = 0; j<syn_err_count[z]; j++) printf(" %d",syn_err_pos[z][j]);
    printf("\n");
  }
  
  return 0;
}


/*** PART B - FINDING THE PROBABILITY OF DECODING ERROR ***/


/* COMPUTE SYNDROME. */

int find_syndrome (int pchk[], int n_errors, int err_pos[])
{
  int syn, i;

  syn = 0;
  for (i = 0; i<n_errors; i++) syn ^= pchk[err_pos[i]];
  
  return syn;
}


/* MAIN PROGRAM TO COMPUTE AND DISPLAY DECODING ERROR PROBABILITIES. */

int main(int argc, char **argv)
{
  int n_pchks, cw_len;
  unsigned z;

  /* Read parity check matrix. */

  if (scanf("%d",&n_pchks)!=1 || scanf("%d",&cw_len)!=1)
  { fprintf(stderr,"Bad header in parity check matrix\n");
    exit(1);
  }

  int pchk[cw_len];

  read_pchk(n_pchks,cw_len,pchk);

  /* Compute syndrome decoding table. */

  int syn_err_count[1U<<n_pchks];

  make_syndrome_table (n_pchks, cw_len, pchk, syn_err_count, 0);

  /* Find and print the histogram of weights for corrected error patterns. */

  int hist[cw_len+1], max_count, imp_count;
  int j;
  
  for (j = 0; j<=cw_len; j++) hist[j] = 0;
  max_count = 0;
  imp_count = 0;

  for (z = 0; z<(1U<<n_pchks); z++) 
  { if (syn_err_count[z]==-1)
    { imp_count += 1;
    }
    else 
    { hist[syn_err_count[z]] += 1;
      if (syn_err_count[z]>max_count) max_count = syn_err_count[z];
    }
  }

  printf("\nHistogram of error counts in syndrome decoding table:\n\n");
  int tot;
  tot = 1;
  for (j = 0; j<=max_count; j++) 
  { printf("%3d  %d / %d\n",j,hist[j],tot);
    tot = tot * (cw_len-j) / (j+1);
  }

  if (imp_count>0)
  { printf("\nNumber of impossible syndromes: %d\n",imp_count);
  }

  /* Find and print the probabilities of incorrect decoding. */

  double f, g;
  int k;

  printf(
"\n\nDecoding error probabilities for given channel error probabilities:\n\n");

  for (k = 1; k<argc; k++)
  {
    if (sscanf(argv[k],"%lf",&f)!=1) 
    { fprintf(stderr,"Argument %s ignored\n",argv[k]);
      continue;
    }

    g = 0.0;
    for (j = 0; j<=max_count; j++)
    { g += hist[j] * pow(f,j) * pow(1-f,cw_len-j);
    }
    printf("%.4f %.6f\n",f,1-g);
  }
  
  return 0;
}