/*
  Author:  Pate Williams (c) 1997

  Exercise
  "5.11 It can be shown that the matrix H shown
  below is a parity-check matrix for a [15, 7, 5]
  code called a BCH code.

      1 0 0 0 1 0 0 1 1 0 1 0 1 1 1
      0 1 0 0 1 1 0 1 0 1 1 1 1 0 0
      0 0 1 0 0 1 1 0 1 0 1 1 1 1 0
  H = 0 0 0 1 0 0 1 1 0 1 0 1 1 1 1
      1 0 0 0 1 1 0 0 0 1 1 0 0 0 1
      0 0 0 1 1 0 0 0 1 1 0 0 0 1 1
      0 0 1 0 1 0 0 1 0 1 0 0 1 0 1
      0 1 1 1 1 0 1 1 1 1 0 1 1 1 1

  Decode, if possible, each of the following
  received vectors r using the syndrome decoding
  method.
  (a) r = (1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0).
  (b) r = (1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0).
  (c) r = (1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0)."
  -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by Douglas
  R. Stinson page 202.
*/

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

#define MAX_ATTEMPTS 11025

void generate(long n, long w, long *e)
/* generates an error vector of weight w */
{
  long i;

  for (i = 0; i < n; i++) e[i] = 0;
  for (i = 0; i < w; i++)
    e[rand() % n] = 1;
}

void mat_vec_mul(long k, long n,
                 long *r, long *s,
                 long **h)
{
  long i, j, sum;

  for (i = 0; i < k; i++) {
    sum = 0;
    for (j = 0; j < n; j++)
      sum += h[i][j] * r[j];
    s[i] = sum % 2;
  }
}

int syndrome(long d, long k, long n,
             long *e, long *r, long *s,
             long *x, long *y,
             long **h)
{
  int equal, nonzero = 0;
  long i, j, l;

  mat_vec_mul(k, n, r, s, h);
  printf("the syndrome vector is:\n");
  for (i = 0; i < k; i++)
    printf("%ld ", s[i]);
  printf("\n");
  for (i = 0; !nonzero && i < k; i++)
    nonzero = s[i] != 0;
  if (!nonzero) {
    for (i = 0; i < n; i++) x[i] = r[i];
    return 1;
  }
  /* generate and test all error vectors of weight 1 */
  for (i = 0; i < n; i++) {
    for (j = 0; j < i; j++) e[j] = 0;
    e[i] = 1;
    for (j = i + 1; j < k; j++) e[j] = 0;
    mat_vec_mul(k, n, e, y, h);
    equal = 1;
    for (j = 0; equal && j < k; j++)
      equal = s[j] == y[j];
    if (equal) {
      for (j = 0; j < n; j++) {
        x[j] = r[j] - e[j];
        if (x[j] < 0) x[j] += 2;
      }
      return 1;
    }
  }
  for (i = 2; i <= d; i++) {
    for (j = 0; j < MAX_ATTEMPTS; j++) {
      generate(n, i, e);
      mat_vec_mul(k, n, e, y, h);
      equal = 1;
      for (l = 0; equal && l < k; l++)
        equal = s[l] == y[l];
      if (equal) {
        for (l = 0; l < n; l++) {
          x[l] = r[l] - e[l];
          if (x[l] < 0) x[l] += 2;
        }
        printf("the error vector is:\n");
        for (l = 0; l < n; l++)
          printf("%ld ", e[l]);
        printf("\n");
        return 1;
      }
    }
  }
  return 0;
}

int main(void)
{
  long H[8][15] = {{1,0,0,0,1,0,0,1,1,0,1,0,1,1,1},
                   {0,1,0,0,1,1,0,1,0,1,1,1,1,0,0},
                   {0,0,1,0,0,1,1,0,1,0,1,1,1,1,0},
                   {0,0,0,1,0,0,1,1,0,1,0,1,1,1,1},
                   {1,0,0,0,1,1,0,0,0,1,1,0,0,0,1},
                   {0,0,0,1,1,0,0,0,1,1,0,0,0,1,1},
                   {0,0,1,0,1,0,0,1,0,1,0,0,1,0,1},
                   {0,1,1,1,1,0,1,1,1,1,0,1,1,1,1}};
  long r[3][15] = {{1,1,0,0,0,0,0,0,0,0,0,0,0,0,0},
                   {1,1,0,1,1,1,1,0,1,0,1,1,0,0,0},
                   {1,0,1,0,1,0,0,1,0,1,1,0,0,0,0}};
  long d = 5, i, j, k = 7, n = 15;
  long e[15], s[15], x[15], y[15];
  long **h = calloc(k + 1, sizeof(long *));

  if (!h) {
    fprintf(stderr, "*error*\ninsufficient memory");
    exit(1);
  }
  printf("the parity check matrix is:\n");
  for (i = 0; i < k + 1; i++) {
    h[i] = calloc(n, sizeof(long *));
    if (!h[i]) {
      fprintf(stderr, "*error*\ninsufficient memory");
      exit(1);
    }
    for (j = 0; j < n; j++) {
      printf("%ld ", H[i][j]);
      h[i][j] = H[i][j];
    }
    printf("\n");
  }
  for (i = 0; i < 3; i++) {
    printf("the vector to be decoded is:\n");
    for (j = 0; j < n; j++)
      printf("%ld ", r[i][j]);
    printf("\n");
    if (syndrome(d, k + 1, n, e, r[i], s, x, y, h)) {
      printf("the decoded vector is:\n");
      for (j = 0; j < n; j++)
        printf("%ld ", x[j]);
      printf("\n");
    }
    else
      printf("vector was not decoded\n");
  }
  for (i = 0; i < k + 1; i++) free(h[i]);
  free(h);
  return 0;
}