/*
  Author:  Pate Williams (c) 1997

  Exercise 8.1
  "Suppose the Blom Scheme with k = 1 is implemented
  with four users U, V, W, and X. Suppose that
  p = 7873, r_u = 2365, r_v = 6648, r_w = 1837, and
  r_x = 2186. The secret g polnomials are as follows:
  g_u(x) = 6018 + 6351x
  g_v(x) = 3749 + 7121x
  g_w(x) = 7061 + 7802x
  g_x(x) = 635 + 6828x
  (a) Compute the key for each pair of users,
      verifying that each pair of users obtains a
      common key (that is K_u,v = K_v,u, etc.).
  (b) Show that W and X together can compute
      K_u,v." -Douglas R. Stinson-
  See "Cryptography: Theory and Practice" by
  Douglas R. Stinson page 281.
*/

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

#define SIZE 4

long Extended_Euclidean(long b, long n)
{
  long b0 = b, n0 = n, t = 1, t0 = 0, temp, q, r;

  q = n0 / b0;
  r = n0 - q * b0;
  while (r > 0) {
    temp = t0 - q * t;
    if (temp >= 0) temp = temp % n;
    else temp = n - (- temp % n);
    t0 = t;
    t = temp;
    n0 = b0;
    b0 = r;
    q = n0 / b0;
    r = n0 - q * b0;
  }
  if (b0 != 1) return 0;
  else return t % n;
}


void Blom_Scheme(long k, long m, long p,
                 long *r, long **g, long **K)
{
  long i, j, l, s, t;

  for (i = 0; i < m; i++) {
    for (j = i + 1; j < m; j++) {
      t = r[j];
      s = g[i][k];
      for (l = k - 1; l >= 0; l--)
        s = t * s + g[i][l];
      K[i][j] = s % p;
    }
  }
  for (i = 0; i < m; i++) {
    for (j = 0; j < i; j++) {
      t = r[j];
      s = g[i][k];
      for (l = k - 1; l >= 0; l--)
        s = t * s + g[i][l];
      s %= p;
      if (s != K[j][i]) {
        fprintf(stderr, "fatal error\nin key calculation\n");
        exit(1);
      }
      K[i][j] = s;
    }
  }
}

void gaussian_elimination(long n, long p,
                          long *b, long *x,
                          long **m)
{
  int found;
  long *d = calloc(n, sizeof(long)), ck, dj;
  long i, j, k, l, sum, t;

  if (!d) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    fprintf(stderr, "from gaussian_elimination\n");
    exit(1);
  }
  for (j = 0; j < n; j++) {
    found = 0, i = j;
    while (!found && i < n) {
      found = m[i][j] != 0;
      if (!found) i++;
    }
    if (!found) {
      fprintf(stderr, "fatal error\nnon-invertible matrix\n");
      fprintf(stderr, "from gaussian_elimination\n");
      for (k = 0; k < n; k++) {
        for (l = 0; l < n; l++)
          printf("%ld ", m[k][l]);
        printf("\n");
      }
      exit(1);
    }
    if (i > j) {
      /* swap elements */
      for (l = j; l < n; l++)
        t = m[i][l], m[i][l] = m[j][l], m[j][l] = t;
      t = b[i], b[i] = b[j], b[j] = t;
    }
    dj = d[j] = Extended_Euclidean(m[j][j], p);
    if (dj == 0) {
      fprintf(stderr, "fatal error\nnon-invertlible element\n");
      fprintf(stderr, "from gaussian elimination");
      fprintf(stderr, "element %ld mod %ld\n", m[j][j], p);
      exit(1);
    }
    for (k = j + 1; k < n; k++) {
      ck = (dj * m[k][j]) % p;
      for (l = j + 1; l < n; l++) {
        m[k][l] = (m[k][l] - (ck * m[j][l]) % p) % p;
        if (m[k][l] < 0) m[k][l] += p;
      }
      b[k] = (b[k] - (ck * b[j]) % p) % p;
      if (b[k] < 0) b[k] += p;
    }
  }
  for (i = n - 1; i >= 0; i--) {
    sum = 0;
    for (j = i + 1; j < n; j++)
      sum += (m[i][j] * x[j]) % p;
    if (sum < 0) sum += p;
    x[i] = (d[i] * ((b[i] - sum) % p)) % p;
    if (x[i] < 0) x[i] += p;
  }
}

void print_matrix(long m, long n, long **matrix)
{
  long i, j;

  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++)
      printf("%4ld ", matrix[i][j]);
    printf("\n");
  }
}

int main(void)
{
  long i, j, k = 1, m = 4, n = 3, p = 7873, s;
  long b[SIZE], c[SIZE], r[SIZE], x[SIZE];
  long **a, **g, **h, **K;
  long A = 2537, B = 4128, C = 6701;

  a = calloc(m, sizeof(long *));
  g = calloc(m, sizeof(long *));
  h = calloc(m, sizeof(long *));
  K = calloc(m, sizeof(long *));
  if (!a || !g || !h || !K) {
    fprintf(stderr, "fatal error\ninsufficient memory\n");
    exit(1);
  }
  for (i = 0; i < m; i++) {
    a[i] = calloc(n, sizeof(long));
    g[i] = calloc(k, sizeof(long));
    h[i] = calloc(n, sizeof(long));
    K[i] = calloc(m, sizeof(long));
    if (!a[i] || !g[i] || !h[i] || !K[i]) {
      fprintf(stderr, "fatal error\ninsufficient memory\n");
      exit(1);
    }
  }
/*
  g[0][0] = 6018, g[0][1] = 6351;
  g[1][0] = 3749, g[1][1] = 7121;
  g[2][0] = 7601, g[2][1] = 7802;
  g[3][0] =  635, g[3][1] = 6828;
*/
  r[0] = 2365;
  r[1] = 6648;
  r[2] = 1837;
  r[3] = 2186;
  for (i = 0; i < m; i++) {
    g[i][0] = (A + (B * r[i]) % p) % p;
    g[i][1] = (B + (C * r[i]) % p) % p;
  }
  Blom_Scheme(k, m, p, r, g, K);
  printf("the g matrix is:\n");
  print_matrix(m, k + 1, g);
  printf("the K matrix is:\n");
  print_matrix(m, m, K);
  a[0][0] = 1;
  a[0][1] = r[2];
  a[1][1] = 1;
  a[1][2] = r[2];
  a[2][0] = 1;
  a[2][1] = r[3];
  a[3][1] = 1;
  a[3][2] = r[3];
  b[0] = g[2][0];
  b[1] = g[2][1];
  b[2] = g[3][0];
  b[3] = g[3][1];
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      s = 0;
      for (k = 0; k < m; k++)
        s += (a[k][i] * a[k][j]) % p;
      s %= p;
      h[i][j] = s;
    }
  }
  printf("the A matrix:\n");
  print_matrix(m, n, a);
  printf("A_t * A:\n");
  print_matrix(n, n, h);
  for (i = 0; i < n; i++) {
    s = 0;
    for (j = 0; j < m; j++)
      s += (a[j][i] * b[j]) % p;
    s %= p;
    c[i] = s;
  }
  gaussian_elimination(n, p, c, x, h);
  printf("a, b, c in program:\n");
  printf("%ld\n", A);
  printf("%ld\n", B);
  printf("%ld\n", C);
  printf("a, b, c calculated:\n");
  for (i = 0; i < n; i++)
    printf("%ld\n", x[i]);
  h[0][0] = (x[0] + (x[1] * r[0]) % p) % p;
  h[0][1] = (x[1] + (x[2] * r[0]) % p) % p;
  printf("g_u given and calculated:\n");
  printf("%4ld %4ld\n", g[0][0], h[0][0]);
  printf("%4ld %4ld\n", g[0][1], h[0][1]);
  for (i = 0; i < m; i++) {
    free(a[i]);
    free(g[i]);
    free(h[i]);
    free(K[i]);
  }
  free(a);
  free(g);
  free(h);
  free(K);
  return 0;
}